CbmRoot
CbmMuch.cxx
Go to the documentation of this file.
1 //----------------------------------------------------------------------------------------
2 //-------------- CbmMuch -----------
3 //--------------- Modified 18/10/2017 by Omveer Singh -----------
4 //----------------------------------------------------------------------------------------
5 
6 
7 #include "CbmMuch.h"
8 
9 #include "CbmGeoMuch.h"
10 #include "CbmGeoMuchPar.h"
11 #include "CbmGeometryUtils.h"
12 #include "CbmMuchPoint.h"
13 
14 #include "CbmStack.h"
15 #include "FairGeoInterface.h"
16 #include "FairGeoLoader.h"
17 #include "FairGeoNode.h"
18 #include "FairGeoRootBuilder.h"
19 #include "FairRootManager.h"
20 #include "FairRun.h"
21 #include "FairRuntimeDb.h"
22 #include "FairVolume.h"
23 
24 #include "TClonesArray.h"
25 #include "TGeoMCGeometry.h"
26 #include "TKey.h"
27 #include "TObjArray.h"
28 #include "TParticle.h"
29 #include "TVirtualMC.h"
30 
31 #include "FairGeoMedia.h"
32 #include "FairGeoMedium.h"
33 #include "TGeoBBox.h"
34 #include "TGeoBoolNode.h"
35 #include "TGeoCompositeShape.h"
36 #include "TGeoCone.h"
37 #include "TGeoManager.h"
38 #include "TGeoMatrix.h"
39 #include "TGeoTube.h"
40 #include "TGeoVolume.h"
41 
42 #include "CbmMuchAddress.h"
43 #include "CbmMuchGeoScheme.h"
44 #include "CbmMuchLayer.h"
45 #include "CbmMuchLayerSide.h"
46 #include "CbmMuchModule.h"
47 #include "CbmMuchModuleGemRadial.h"
48 #include "CbmMuchStation.h"
49 #include "TGeoMatrix.h"
50 
51 #include <cassert>
52 #include <fstream>
53 #include <iostream>
54 
55 #include "TGeoArb8.h"
56 #include <TFile.h>
57 
58 using std::cout;
59 using std::endl;
60 
61 ClassImp(CbmMuch)
62 
63  // ----- Default constructor -------------------------------------------
64  CbmMuch::CbmMuch()
65  : FairDetector()
66  , fTrackID(-1)
67  , fVolumeID(-1)
68  , fFlagID(0)
69  , fPosIn()
70  , fPosOut()
71  , fMomIn()
72  , fMomOut()
73  , fTime(0.)
74  , fLength(0.)
75  , fELoss(0.)
76  , fPosIndex(0)
77  , fMuchCollection(new TClonesArray("CbmMuchPoint"))
78  , kGeoSaved(kFALSE)
79  , flGeoPar(new TList())
80  , fPar(NULL)
81  , fVolumeName("")
82  , fCombiTrans() {
83  ResetParameters();
84  flGeoPar->SetName(GetName());
85  fVerboseLevel = 1;
86 }
87 // -------------------------------------------------------------------------
88 
89 
90 // ----- Standard constructor ------------------------------------------
91 CbmMuch::CbmMuch(const char* name, Bool_t active)
92  : FairDetector(name, active)
93  , fTrackID(-1)
94  , fVolumeID(-1)
95  , fFlagID(0)
96  , fPosIn()
97  , fPosOut()
98  , fMomIn()
99  , fMomOut()
100  , fTime(0.)
101  , fLength(0.)
102  , fELoss(0.)
103  , fPosIndex(0)
104  , fMuchCollection(new TClonesArray("CbmMuchPoint"))
105  , kGeoSaved(kFALSE)
106  , flGeoPar(new TList())
107  , fPar(NULL)
108  , fVolumeName("")
109  , fCombiTrans() {
110 
111  ResetParameters();
112  flGeoPar->SetName(GetName());
113  fVerboseLevel = 1;
114 }
115 
116 
117 // ----- Destructor ----------------------------------------------------
118 CbmMuch::~CbmMuch() {
119 
120  if (flGeoPar) delete flGeoPar;
121  if (fMuchCollection) {
122  fMuchCollection->Delete();
123  delete fMuchCollection;
124  }
125 }
126 
127 
128 // ----- Public method ProcessHits --------------------------------------
129 
130 
131 Bool_t CbmMuch::ProcessHits(FairVolume* vol) {
132  // cout<<" called process Hit****************** "<<endl;
133  if (gMC->IsTrackEntering()) {
134  fELoss = 0.;
135  fTime = gMC->TrackTime() * 1.0e09;
136  fLength = gMC->TrackLength();
137  gMC->TrackPosition(fPosIn);
138  gMC->TrackMomentum(fMomIn);
139  }
140 
141  // Sum energy loss for all steps in the active volume
142  fELoss += gMC->Edep();
143 
144  // Set additional parameters at exit of active volume. Create CbmMuchPoint.
145  if (gMC->IsTrackExiting() || gMC->IsTrackStop()
146  || gMC->IsTrackDisappeared()) {
147  fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
148  fVolumeID = vol->getMCid();
149  Int_t fDetectorId = GetDetId(vol);
150  // cout<<" check det ID 2 "<<fDetectorId<<endl;
151  if (fVerboseLevel > 2) {
152  printf(" TrackId: %i", fTrackID);
153  printf(" System: %i", CbmMuchAddress::GetSystemIndex(fDetectorId));
154  printf(" Station: %i", CbmMuchAddress::GetStationIndex(fDetectorId));
155  printf(" Layer: %i", CbmMuchAddress::GetLayerIndex(fDetectorId));
156  printf(" Side: %i", CbmMuchAddress::GetLayerSideIndex(fDetectorId));
157  printf(" Module: %i", CbmMuchAddress::GetModuleIndex(fDetectorId));
158  printf(" Vol %i \n", fVolumeID);
159  }
160  Int_t iStation = CbmMuchAddress::GetStationIndex(fDetectorId);
161  gMC->TrackPosition(fPosOut);
162  gMC->TrackMomentum(fMomOut);
163 
164  assert(iStation >= 0 && iStation < fPar->GetNStations());
165  CbmMuchStation* station =
166  (CbmMuchStation*) fPar->GetStations()->At(iStation);
167  //cout<<" track # "<<fTrackID<<" Rmin "<<station->GetRmin()<<" Rmax "<<station->GetRmax()<<" in perp "<<fPosIn.Perp()<<" out perp "<<fPosOut.Perp()<<" eloss "<<fELoss<<endl;
168  if (fPosIn.Perp() > station->GetRmax()) { return kFALSE; }
169  if (fPosOut.Perp() > station->GetRmax()) { return kFALSE; }
170  if (fPosIn.Perp() < station->GetRmin()) { return kFALSE; }
171  if (fPosOut.Perp() < station->GetRmin()) { return kFALSE; }
172 
173 
174  if (fELoss == 0.) return kFALSE;
175  AddHit(fTrackID,
176  fDetectorId,
177  TVector3(fPosIn.X(), fPosIn.Y(), fPosIn.Z()),
178  TVector3(fPosOut.X(), fPosOut.Y(), fPosOut.Z()),
179  TVector3(fMomIn.Px(), fMomIn.Py(), fMomIn.Pz()),
180  TVector3(fMomOut.Px(), fMomOut.Py(), fMomOut.Pz()),
181  fTime,
182  fLength,
183  fELoss);
184 
185  //if (fPosOut.Z()>250) printf("%f\n",fPosOut.Z());
186 
187  // Increment number of MuchPoints for this track
188  CbmStack* stack = (CbmStack*) gMC->GetStack();
189  stack->AddPoint(ECbmModuleId::kMuch);
190 
191  ResetParameters();
192  }
193  return kTRUE;
194 }
195 
196 
197 //-------------------------------------------------------------------------
198 
199 Int_t CbmMuch::GetDetId(FairVolume* vol) {
200  TString name = vol->GetName();
201  Char_t cStationNr[3] = {name[11], name[12], ' '};
202  Int_t iStation = atoi(cStationNr) - 1;
203  Int_t iLayer = TString(name[18]).Atoi() - 1;
204  Int_t iSide = (name[19] == 'b') ? 1 : 0;
205  Char_t cModuleNr[4] = {name[26], name[27], name[28], ' '};
206  Int_t iModule = atoi(cModuleNr) - 1;
207  if (iSide != 1 && iSide != 0) printf("side = %i", iSide);
208  Int_t detectorId =
209  CbmMuchAddress::GetAddress(iStation, iLayer, iSide, iModule);
210  assert(CbmMuchAddress::GetStationIndex(detectorId) == iStation);
211  assert(CbmMuchAddress::GetLayerIndex(detectorId) == iLayer);
212  assert(CbmMuchAddress::GetLayerSideIndex(detectorId) == iSide);
213  assert(CbmMuchAddress::GetModuleIndex(detectorId) == iModule);
214  assert(detectorId > 0);
215 
216  return detectorId;
217 }
218 
219 // -------------------------------------------------------------------------
220 void CbmMuch::BeginEvent() {}
221 
222 
223 // -------------------------------------------------------------------------
224 void CbmMuch::EndOfEvent() {
225  if (fVerboseLevel) Print("");
226  fMuchCollection->Delete();
227  ResetParameters();
228 }
229 
230 // -------------------------------------------------------------------------
231 void CbmMuch::Register() {
232  FairRootManager::Instance()->Register(
233  "MuchPoint", GetName(), fMuchCollection, kTRUE);
234 }
235 
236 // -------------------------------------------------------------------------
237 TClonesArray* CbmMuch::GetCollection(Int_t iColl) const {
238  if (iColl == 0)
239  return fMuchCollection;
240  else
241  return NULL;
242 }
243 // -------------------------------------------------------------------------
244 
245 
246 // -------------------------------------------------------------------------
247 void CbmMuch::Print(Option_t*) const {
248  Int_t nHits = fMuchCollection->GetEntriesFast();
249  LOG(info) << fName << ": " << nHits << " points registered in this event.";
250 }
251 
252 // -------------------------------------------------------------------------
253 void CbmMuch::Reset() {
254  fMuchCollection->Delete();
255  ResetParameters();
256 }
257 
258 // -------------------------------------------------------------------------
259 void CbmMuch::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset) {
260  Int_t nEntries = cl1->GetEntriesFast();
261  LOG(info) << fName << ": " << nEntries << " entries to add.";
262  TClonesArray& clref = *cl2;
263  CbmMuchPoint* oldpoint = NULL;
264  for (Int_t i = 0; i < nEntries; i++) {
265  oldpoint = (CbmMuchPoint*) cl1->At(i);
266  Int_t index = oldpoint->GetTrackID() + offset;
267  oldpoint->SetTrackID(index);
268  new (clref[fPosIndex]) CbmMuchPoint(*oldpoint);
269  fPosIndex++;
270  }
271  LOG(info) << fName << ": " << cl2->GetEntriesFast() << " merged entries.";
272 }
273 // -------------------------------------------------------------------------
274 
275 
276 // ----- Private method AddHit --------------------------------------------
277 CbmMuchPoint* CbmMuch::AddHit(Int_t trackID,
278  Int_t detID,
279  TVector3 posIn,
280  TVector3 posOut,
281  TVector3 momIn,
282  TVector3 momOut,
283  Double_t time,
284  Double_t length,
285  Double_t eLoss) {
286 
287  TClonesArray& clref = *fMuchCollection;
288  Int_t size = clref.GetEntriesFast();
289  if (fVerboseLevel > 1)
290  LOG(info) << fName << ": Adding Point at (" << posIn.X() << ", "
291  << posIn.Y() << ", " << posIn.Z() << ") cm, detector " << detID
292  << ", track " << trackID << ", energy loss " << eLoss * 1e06
293  << " keV";
294 
295  return new (clref[size]) CbmMuchPoint(
296  trackID, detID, posIn, posOut, momIn, momOut, time, length, eLoss);
297 }
298 
299 // ----- ConstructGeometry -----------------------------------------------
300 void CbmMuch::ConstructGeometry() {
301 
302  TString fileName = GetGeometryFileName();
303 
304  // --- Only ROOT geometries are supported
305  if (!fileName.EndsWith(".root")) {
306  LOG(fatal) << GetName() << ": Geometry format of file " << fileName.Data()
307  << " not supported.";
308  }
309 
310  if (fileName.Contains("mcbm")) {
311 
312  fFlagID = 1;
313  LOG(info) << "mcbm geometry found ";
314  }
315  ConstructRootGeometry();
316 }
317 // -------------------------------------------------------------------------
318 
319 
320 void CbmMuch::ConstructRootGeometry(TGeoMatrix*) {
321  FairRun* fRun = FairRun::Instance();
322  if (!fRun) {
323  Fatal("CreateGeometry", "No FairRun found");
324  return;
325  }
326  FairRuntimeDb* rtdb = FairRuntimeDb::instance();
327  fPar = (CbmGeoMuchPar*) (rtdb->getContainer("CbmGeoMuchPar"));
328  TObjArray* stations = fPar->GetStations();
329 
330  //-------------------------------------------------------
331 
332  CbmMuchGeoScheme* fGeoScheme = CbmMuchGeoScheme::Instance();
333 
334  fGeoScheme->Init(stations, fFlagID);
335 
336  TGeoMatrix* tempMatrix {nullptr};
337  if (Cbm::GeometryUtils::IsNewGeometryFile(
338  fgeoName, fVolumeName, &tempMatrix)) {
339  LOG(info) << "Importing MUCH geometry from ROOT file " << fgeoName.Data();
340  Cbm::GeometryUtils::ImportRootGeometry(fgeoName, this);
341  } else {
342  LOG(info) << "Constructing MUCH geometry from ROOT file "
343  << fgeoName.Data();
344  FairModule::ConstructRootGeometry();
345  }
346 
347 
348  TObjArray* fSensNodes = fPar->GetGeoSensitiveNodes();
349  TObjArray* fPassNodes = fPar->GetGeoPassiveNodes();
350  TGeoNode* ncave = gGeoManager->GetTopNode();
351 
352  fGeoScheme->ExtractGeoParameter(ncave, fVolumeName.Data());
353 
354  TString objName = fVolumeName + "_0";
355  TGeoNode* nmuch = (TGeoNode*) ncave->GetNodes()->FindObject(objName);
356  fPassNodes->Add(nmuch);
357 
358 
359  TObjArray* nodes = nmuch->GetNodes();
360 
361  for (Int_t i = 0; i < nodes->GetEntriesFast(); i++) {
362  TGeoNode* node = (TGeoNode*) nodes->At(i);
363  TString nodeName = node->GetName();
364 
365 
366  TObjArray* nodes1 = node->GetNodes();
367 
368  for (Int_t j = 0; j < nodes1->GetEntriesFast(); j++) {
369  TGeoNode* node1 = (TGeoNode*) nodes1->At(j);
370  TString node1Name = node1->GetName();
371 
372  if (node1Name.Contains("absblock")) fPassNodes->Add(node);
373  if (node1Name.Contains("muchstation")) {
374 
375  TObjArray* layers = node1->GetNodes();
376  for (Int_t l = 0; l < layers->GetEntriesFast(); l++) {
377  TGeoNode* layer = (TGeoNode*) layers->At(l);
378 
379  if (!TString(layer->GetName()).Contains("layer")) continue;
380  TObjArray* layerNodes = layer->GetNodes();
381  for (Int_t m = 0; m < layerNodes->GetEntriesFast(); m++) {
382  TGeoNode* layerNode = (TGeoNode*) layerNodes->At(m);
383  TString layerNodeName = layerNode->GetName();
384 
385  if (layerNodeName.Contains("active")) fSensNodes->Add(layerNode);
386 
387  if (layerNodeName.Contains("support")) fPassNodes->Add(layerNode);
388 
389  if (layerNodeName.Contains("cool")) fPassNodes->Add(layerNode);
390  }
391  }
392  }
393  }
394  }
395 
396  fPar->setChanged();
397  fPar->setInputVersion(fRun->GetRunId(), 1);
398 }
399 
400 
401 // ----- CheckIfSensitive -------------------------------------------------
402 Bool_t CbmMuch::CheckIfSensitive(std::string name) {
403  TString tsname = name;
404 
405 
406  if (tsname.Contains("active")) {
407  LOG(debug1) << "CbmMuch::CheckIfSensitive: Register active volume: "
408  << tsname;
409  return kTRUE;
410  }
411  return kFALSE;
412 }
CbmMuchLayerSide.h
CbmMuch::CopyClones
virtual void CopyClones(TClonesArray *cl1, TClonesArray *cl2, Int_t offset)
Definition: CbmMuch.cxx:259
CbmMuchModule.h
CbmMuchGeoScheme
Definition: CbmMuchGeoScheme.h:43
CbmMuchPoint
Definition: CbmMuchPoint.h:21
CbmMuch::fMomIn
TLorentzVector fMomIn
position
Definition: CbmMuch.h:117
CbmMuchStation
Definition: CbmMuchStation.h:22
CbmMuch::fVolumeID
Int_t fVolumeID
track index
Definition: CbmMuch.h:114
CbmMuch::fVolumeName
TString fVolumeName
parameter container
Definition: CbmMuch.h:127
CbmMuch::CbmMuch
CbmMuch()
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmMuch::fPosOut
TLorentzVector fPosOut
Definition: CbmMuch.h:116
CbmMuch::Reset
virtual void Reset()
Definition: CbmMuch.cxx:253
CbmMuch::fTrackID
Int_t fTrackID
Definition: CbmMuch.h:113
CbmGeoMuchPar::GetStations
TObjArray * GetStations()
Definition: CbmGeoMuchPar.h:40
CbmMuch::fMomOut
TLorentzVector fMomOut
Definition: CbmMuch.h:117
CbmMuchGeoScheme::Init
void Init(TObjArray *stations, Int_t flag)
Definition: CbmMuchGeoScheme.cxx:121
CbmGeoMuchPar
Definition: CbmGeoMuchPar.h:25
CbmMuchAddress::GetSystemIndex
static Int_t GetSystemIndex(Int_t address)
Definition: CbmMuchAddress.h:100
CbmMuchGeoScheme::ExtractGeoParameter
void ExtractGeoParameter(TGeoNode *muchNode, const char *volumeName)
Definition: CbmMuchGeoScheme.cxx:463
CbmMuch::flGeoPar
TList * flGeoPar
Definition: CbmMuch.h:125
CbmMuchModuleGemRadial.h
CbmStack.h
CbmGeoMuchPar::GetGeoSensitiveNodes
TObjArray * GetGeoSensitiveNodes()
Definition: CbmGeoMuchPar.h:38
CbmMuchAddress::GetModuleIndex
static Int_t GetModuleIndex(Int_t address)
Definition: CbmMuchAddress.h:112
CbmStack::AddPoint
void AddPoint(ECbmModuleId iDet)
Definition: CbmStack.cxx:383
CbmMuch::~CbmMuch
virtual ~CbmMuch()
Definition: CbmMuch.cxx:118
CbmMuch::Print
virtual void Print(Option_t *) const
Definition: CbmMuch.cxx:247
CbmMuch::GetCollection
virtual TClonesArray * GetCollection(Int_t iColl) const
Definition: CbmMuch.cxx:237
CbmMuchPoint.h
CbmGeoMuchPar.h
CbmMuch.h
CbmMuch::ConstructGeometry
virtual void ConstructGeometry()
Definition: CbmMuch.cxx:300
Cbm::GeometryUtils::IsNewGeometryFile
Bool_t IsNewGeometryFile(TString &filename)
Definition: CbmGeometryUtils.cxx:133
CbmMuchStation::GetRmax
Double_t GetRmax() const
Definition: CbmMuchStation.h:52
CbmMuch
Definition: CbmMuch.h:27
CbmMuchGeoScheme::Instance
static CbmMuchGeoScheme * Instance()
Definition: CbmMuchGeoScheme.cxx:113
CbmMuch::BeginEvent
virtual void BeginEvent()
Definition: CbmMuch.cxx:220
CbmMuch::fLength
Double32_t fLength
time
Definition: CbmMuch.h:119
CbmMuch::fPosIndex
Int_t fPosIndex
energy loss
Definition: CbmMuch.h:122
CbmMuch::fFlagID
Int_t fFlagID
volume id
Definition: CbmMuch.h:115
CbmMuchLayer.h
CbmMuchAddress::GetLayerIndex
static Int_t GetLayerIndex(Int_t address)
Definition: CbmMuchAddress.h:106
CbmMuchAddress::GetStationIndex
static Int_t GetStationIndex(Int_t address)
Definition: CbmMuchAddress.h:103
CbmMuch::ConstructRootGeometry
virtual void ConstructRootGeometry(TGeoMatrix *shift=NULL)
Definition: CbmMuch.cxx:320
CbmMuch::GetDetId
Int_t GetDetId(FairVolume *vol)
Definition: CbmMuch.cxx:199
CbmMuch::AddHit
CbmMuchPoint * AddHit(Int_t trackID, Int_t detID, TVector3 posIn, TVector3 posOut, TVector3 momIn, TVector3 momOut, Double_t time, Double_t length, Double_t eLoss)
Definition: CbmMuch.cxx:277
CbmMuch::Register
virtual void Register()
Definition: CbmMuch.cxx:231
CbmMuchAddress.h
CbmGeoMuch.h
m
__m128 m
Definition: L1/vectors/P4_F32vec4.h:26
ClassImp
ClassImp(CbmMuch) CbmMuch
Definition: CbmMuch.cxx:61
CbmMuch::fPar
CbmGeoMuchPar * fPar
Definition: CbmMuch.h:126
CbmGeoMuchPar::GetGeoPassiveNodes
TObjArray * GetGeoPassiveNodes()
Definition: CbmGeoMuchPar.h:39
Cbm::GeometryUtils::ImportRootGeometry
void ImportRootGeometry(TString &filename, FairModule *mod, TGeoMatrix *mat)
Definition: CbmGeometryUtils.cxx:140
CbmMuch::ProcessHits
virtual Bool_t ProcessHits(FairVolume *vol=0)
Definition: CbmMuch.cxx:131
ECbmModuleId::kMuch
@ kMuch
Muon detection system.
CbmMuchGeoScheme.h
CbmMuch::ResetParameters
void ResetParameters()
Definition: CbmMuch.h:158
CbmStack
Definition: CbmStack.h:45
CbmMuchAddress::GetAddress
static UInt_t GetAddress(Int_t station=0, Int_t layer=0, Int_t side=0, Int_t module=0, Int_t sector=0, Int_t channel=0)
Definition: CbmMuchAddress.cxx:43
CbmMuchStation::GetRmin
Double_t GetRmin() const
Definition: CbmMuchStation.h:51
CbmMuch::fMuchCollection
TClonesArray * fMuchCollection
Definition: CbmMuch.h:123
CbmMuch::fPosIn
TLorentzVector fPosIn
Definition: CbmMuch.h:116
CbmMuchStation.h
CbmMuchAddress::GetLayerSideIndex
static Int_t GetLayerSideIndex(Int_t address)
Definition: CbmMuchAddress.h:109
CbmMuch::CheckIfSensitive
Bool_t CheckIfSensitive(std::string name)
Definition: CbmMuch.cxx:402
CbmGeometryUtils.h
CbmMuch::EndOfEvent
virtual void EndOfEvent()
Definition: CbmMuch.cxx:224
CbmMuch::fELoss
Double32_t fELoss
length
Definition: CbmMuch.h:120
CbmMuch::fTime
Double32_t fTime
momentum
Definition: CbmMuch.h:118