13 #include <FairMCEventHeader.h>
14 #include <FairRunSim.h>
16 #include <FairGeoBuilder.h>
17 #include <FairGeoMedia.h>
19 #include "TClonesArray.h"
20 #include "TGeoManager.h"
22 #include "TGeoPhysicalNode.h"
23 #include "TGeoVolume.h"
25 #include "TObjArray.h"
26 #include "TParticle.h"
27 #include "TVirtualMC.h"
30 #include <boost/lexical_cast.hpp>
31 #include <boost/regex.hpp>
48 , fTofCollection(new TClonesArray(
"CbmTofPoint"))
52 , fbOnePointPerTrack(kFALSE)
53 , fbIsNewTrack(kFALSE)
55 , fCurrentNodePath(
"")
56 , fCurrentModuleType(0)
57 , fCurrentModuleIndex(0)
58 , fCurrentCounterIndex(0)
61 , fInactiveCounterIDs()
64 , fbProcessAnyTrack(kFALSE)
65 , fbAllCountersInactive(kFALSE) {
82 , fTofCollection(new TClonesArray(
"CbmTofPoint"))
86 , fbOnePointPerTrack(kFALSE)
87 , fbIsNewTrack(kFALSE)
89 , fCurrentNodePath(
"")
90 , fCurrentModuleType(0)
91 , fCurrentModuleIndex(0)
92 , fCurrentCounterIndex(0)
95 , fInactiveCounterIDs()
98 , fbProcessAnyTrack(kFALSE)
99 , fbAllCountersInactive(kFALSE) {
114 if ((CounterInBeam.second).second) {
delete (CounterInBeam.second).second; }
120 FairDetector::Initialize();
124 Bool_t isSimulation = kTRUE;
143 std::get<0>(*itInactiveCounter),
144 std::get<1>(*itInactiveCounter),
145 std::get<2>(*itInactiveCounter),
172 Int_t iNCells = tCurrentPoint->
GetNCells();
173 tCurrentPoint->SetTime(tCurrentPoint->GetTime() / iNCells);
174 tCurrentPoint->SetLength(tCurrentPoint->GetLength() / iNCells);
175 tCurrentPoint->SetX(tCurrentPoint->GetX() / iNCells);
176 tCurrentPoint->SetY(tCurrentPoint->GetY() / iNCells);
177 tCurrentPoint->SetZ(tCurrentPoint->GetZ() / iNCells);
178 tCurrentPoint->
SetPx(tCurrentPoint->GetPx() / iNCells);
179 tCurrentPoint->
SetPy(tCurrentPoint->GetPy() / iNCells);
180 tCurrentPoint->
SetPz(tCurrentPoint->GetPz() / iNCells);
182 if (1 == tCurrentPoint->GetNLinks()) {
184 FairLink tLinkToTrack = tCurrentPoint->GetLink(0);
185 tLinkToTrack.SetLink(
187 tCurrentPoint->SetLink(tLinkToTrack);
195 FairMCEventHeader* tEventHeader = FairRunSim::Instance()->GetMCEventHeader();
197 Double_t dTargetVertexT = tEventHeader->GetT();
199 Double_t dGlobalTargetCoordinates[3] = {
200 tEventHeader->GetX(), tEventHeader->GetY(), tEventHeader->GetZ()};
201 Double_t dGlobalTargetCoordinates1[3] = {
202 tEventHeader->GetX(), tEventHeader->GetY(), tEventHeader->GetZ() + 1.};
203 Double_t dGlobalCounterCoordinates[3] = {0., 0., 0.};
204 Double_t dLocalTargetCoordinates[3] = {0., 0., 0.};
205 Double_t dLocalTargetCoordinates1[3] = {0., 0., 0.};
206 Double_t dLocalCounterCoordinates[3] = {0., 0., 0.};
208 TGeoPhysicalNode* tCurrentNode(NULL);
209 Int_t iModuleType(0);
210 Int_t iModuleIndex(0);
211 Int_t iCounterIndex(0);
212 Int_t iUniqueCounterId(0);
217 iModuleType = std::get<0>(CounterInBeam.first);
218 iModuleIndex = std::get<1>(CounterInBeam.first);
219 iCounterIndex = std::get<2>(CounterInBeam.first);
227 tCurrentNode = (CounterInBeam.second).second;
230 tCurrentNode->GetMatrix()->MasterToLocal(dGlobalTargetCoordinates,
231 dLocalTargetCoordinates);
232 tCurrentNode->GetMatrix()->MasterToLocal(dGlobalTargetCoordinates1,
233 dLocalTargetCoordinates1);
236 dLocalCounterCoordinates[0] =
237 dLocalTargetCoordinates[0]
238 - dLocalTargetCoordinates[2]
239 * (dLocalTargetCoordinates1[0] - dLocalTargetCoordinates[0])
240 / (dLocalTargetCoordinates1[2] - dLocalTargetCoordinates[2]);
241 dLocalCounterCoordinates[1] =
242 dLocalTargetCoordinates[1]
243 - dLocalTargetCoordinates[2]
244 * (dLocalTargetCoordinates1[1] - dLocalTargetCoordinates[1])
245 / (dLocalTargetCoordinates1[2] - dLocalTargetCoordinates[2]);
249 if (tCurrentNode->GetShape()->Contains(dLocalCounterCoordinates)) {
251 tCurrentNode->GetMatrix()->LocalToMaster(dLocalCounterCoordinates,
252 dGlobalCounterCoordinates);
256 Double_t dBeamMomentumLab = FairRunSim::Instance()->GetBeamMom();
257 Double_t dCounterBeamTime = dTargetVertexT;
260 Double_t dCounterTargetDistance =
261 dGlobalCounterCoordinates[2] - dGlobalTargetCoordinates[2];
263 if (0. < dBeamMomentumLab) {
264 Double_t dBeamVelocityLab =
266 / TMath::Sqrt(TMath::Power(dBeamMomentumLab, 2.)
267 + TMath::Power(0.938271998, 2.))
269 dCounterBeamTime += dCounterTargetDistance / dBeamVelocityLab;
275 tBeamPoint->SetDetectorID(iUniqueCounterId);
276 tBeamPoint->SetEventID(
278 tBeamPoint->SetTime(dCounterBeamTime * 1.0e09);
279 tBeamPoint->SetLength(dCounterTargetDistance);
280 tBeamPoint->SetX(dGlobalCounterCoordinates[0]);
281 tBeamPoint->SetY(dGlobalCounterCoordinates[1]);
282 tBeamPoint->SetZ(dGlobalCounterCoordinates[2]);
283 tBeamPoint->
SetPz(dBeamMomentumLab);
303 Int_t iTrackID = gMC->GetStack()->GetCurrentTrackNumber();
305 Double_t dTrackEnergyDeposit = gMC->Edep();
308 Bool_t bCounterPointExists = kFALSE;
314 for (Int_t iPoint =
fTofCollection->GetEntriesFast() - 1; iPoint >= 0;
319 if (tCounterPoint->GetDetectorID() == iCounterID
320 && tCounterPoint->GetTrackID() == iTrackID) {
321 bCounterPointExists = kTRUE;
330 if (gMC->IsTrackEntering()) {
331 Double_t dTrackTime = gMC->TrackTime() * 1.0e09;
332 Double_t dTrackLength = gMC->TrackLength();
333 Double_t dTrackPositionX(0.);
334 Double_t dTrackPositionY(0.);
335 Double_t dTrackPositionZ(0.);
336 gMC->TrackPosition(dTrackPositionX, dTrackPositionY, dTrackPositionZ);
337 Double_t dTrackMomentumX(0.);
338 Double_t dTrackMomentumY(0.);
339 Double_t dTrackMomentumZ(0.);
340 Double_t dTrackEnergy(0.);
342 dTrackMomentumX, dTrackMomentumY, dTrackMomentumZ, dTrackEnergy);
344 if (bCounterPointExists) {
345 tCounterPoint->SetTime(tCounterPoint->GetTime() + dTrackTime);
346 tCounterPoint->SetLength(tCounterPoint->GetLength() + dTrackLength);
347 tCounterPoint->SetEnergyLoss(tCounterPoint->GetEnergyLoss()
348 + dTrackEnergyDeposit);
349 tCounterPoint->SetX(tCounterPoint->GetX() + dTrackPositionX);
350 tCounterPoint->SetY(tCounterPoint->GetY() + dTrackPositionY);
351 tCounterPoint->SetZ(tCounterPoint->GetZ() + dTrackPositionZ);
352 tCounterPoint->
SetPx(tCounterPoint->GetPx() + dTrackMomentumX);
353 tCounterPoint->
SetPy(tCounterPoint->GetPy() + dTrackMomentumY);
354 tCounterPoint->
SetPz(tCounterPoint->GetPz() + dTrackMomentumZ);
360 tCounterPoint->SetTrackID(iTrackID);
361 tCounterPoint->SetDetectorID(iCounterID);
362 tCounterPoint->SetEventID(
365 tCounterPoint->SetTime(dTrackTime);
366 tCounterPoint->SetLength(dTrackLength);
367 tCounterPoint->SetEnergyLoss(dTrackEnergyDeposit);
368 tCounterPoint->SetX(dTrackPositionX);
369 tCounterPoint->SetY(dTrackPositionY);
370 tCounterPoint->SetZ(dTrackPositionZ);
371 tCounterPoint->
SetPx(dTrackMomentumX);
372 tCounterPoint->
SetPy(dTrackMomentumY);
373 tCounterPoint->
SetPz(dTrackMomentumZ);
383 tCounterPoint->SetEnergyLoss(tCounterPoint->GetEnergyLoss()
384 + dTrackEnergyDeposit);
389 if (gMC->IsTrackEntering()) {
391 fTime = gMC->TrackTime() * 1.0e09;
393 gMC->TrackPosition(
fPos);
394 gMC->TrackMomentum(
fMom);
401 if (((0 == gMC->GetStack()->GetCurrentTrack()->GetPdgCode())
403 (gMC->TrackCharge() != 0))
404 && (gMC->IsTrackExiting() || gMC->IsTrackStop()
405 || gMC->IsTrackDisappeared())) {
407 fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
411 LOG(debug2) <<
"CbmTof::TID: " <<
fTrackID;
419 LOG(debug2) << Form(
" x: %6.2f",
fPos.X());
420 LOG(debug2) << Form(
" y: %6.2f",
fPos.Y());
421 LOG(debug2) << Form(
" z: %6.2f",
fPos.Z());
450 if (fVerboseLevel)
Print();
459 FairRootManager::Instance()->Register(
478 LOG(info) << fName <<
": " << nHits <<
" points registered in this event.";
480 if (fVerboseLevel > 1)
481 for (Int_t
i = 0;
i < nHits;
i++)
497 Int_t nEntries = cl1->GetEntriesFast();
498 LOG(info) <<
"CbmTof: " << nEntries <<
" entries to add.";
499 TClonesArray& clref = *cl2;
501 for (Int_t
i = 0;
i < nEntries;
i++) {
503 Int_t index = oldpoint->GetTrackID() + offset;
504 oldpoint->SetTrackID(index);
508 LOG(info) <<
"CbmTof: " << cl2->GetEntriesFast() <<
" merged entries.";
512 TString fileName = GetGeometryFileName();
513 if (fileName.EndsWith(
".root")) {
518 LOG(fatal) <<
"Geometry format " << fileName <<
" not supported.";
524 LOG(info) <<
"Importing TOF geometry from ROOT file " << fgeoName.Data();
527 LOG(info) <<
"Constructing TOF geometry from ROOT file " << fgeoName.Data();
528 FairModule::ConstructRootGeometry();
535 Int_t iCounterIndex) {
537 std::make_tuple(iModuleType, iModuleIndex, iCounterIndex));
545 Int_t iCounterIndex) {
547 std::make_tuple(iModuleType, iModuleIndex, iCounterIndex));
555 Int_t iCounterIndex) {
556 fCountersInBeam[std::make_tuple(iModuleType, iModuleIndex, iCounterIndex)];
564 (CounterInBeam.second).second =
565 new TGeoPhysicalNode(((CounterInBeam.second).first).Data());
574 TString tsname = name;
575 if (tsname.Contains(
"Cell")) {
601 Int_t size = clref.GetEntriesFast();
602 return new (clref[size])