BmnRoot
Loading...
Searching...
No Matches
BmnFD.cxx
Go to the documentation of this file.
1#include <iostream>
2
3#include "TClonesArray.h"
4#include "TGeoMCGeometry.h"
5#include "TGeoManager.h"
6#include "TLorentzVector.h"
7#include "TParticle.h"
8#include "TVirtualMC.h"
9#include "TGeoArb8.h"
10
11#include "FairGeoInterface.h"
12#include "FairGeoLoader.h"
13#include "FairGeoNode.h"
14#include "BmnFDGeo.h"
15#include "FairGeoRootBuilder.h"
16#include "CbmStack.h"
17#include "BmnFD.h"
18#include "BmnFDPoint.h"
19
20#include "TVector3.h"
21#include "FairMCPoint.h"
22#include "FairRootManager.h"
23#include "FairVolume.h"
24// add on for debug
25//#include "FairGeoG3Builder.h"
26#include "FairRuntimeDb.h"
27#include "TObjArray.h"
28#include "FairRun.h"
29
30#include "TParticlePDG.h"
31
32// ----- Default constructor -------------------------------------------
33BmnFD::BmnFD() {
34 fFDCollection = new TClonesArray("BmnFDPoint");
35 volDetector = 0;
36 fPosIndex = 0;
37 // fpreflag = 0;
38 //fpostflag = 0;
39 // fEventID=-1;
40 fVerboseLevel = 1;
41
42}
43// -------------------------------------------------------------------------
44
45// ----- Standard constructor ------------------------------------------
46BmnFD::BmnFD(const char* name, Bool_t active)
47 : FairDetector(name, active) {
48 fFDCollection = new TClonesArray("BmnFDPoint");
49 fPosIndex = 0;
50 volDetector = 0;
51 //fpreflag = 0;
52 //fpostflag = 0;
53 //fEventID=-1;
54 fVerboseLevel = 1;
55}
56// -------------------------------------------------------------------------
57
58
59
60// ----- Destructor ----------------------------------------------------
61BmnFD::~BmnFD() {
62 if (fFDCollection) {
63 fFDCollection->Delete();
64 delete fFDCollection;
65 }
66
67}
68// -------------------------------------------------------------------------
69
70// ----- Public method Intialize ---------------------------------------
71void BmnFD::Initialize() {
72 // Init function
73
74 FairDetector::Initialize();
75 //FairRun* sim = FairRun::Instance();
76 //FairRuntimeDb* rtdb=sim->GetRuntimeDb();
77}
78// -------------------------------------------------------------------------
79void BmnFD::BeginEvent(){
80 // Begin of the event
81
82}
83
84
85
86// ----- Public method ProcessHits --------------------------------------
87Bool_t BmnFD::ProcessHits(FairVolume* vol) {
88
89 // if (TMath::Abs(gMC->TrackCharge()) <= 0) return kFALSE;
90
91 static const Double_t dP = 1.032 ;
92 static const Double_t BirkC1 = 0.013/dP;
93 static const Double_t BirkC2 = 9.6e-6/(dP * dP);
94
95 static Double_t lightYield;
96 static Double_t timeIn;
97 static Double_t timeOut;
98 static Double_t lengthtrack;
99 //static Double_t de_dx;
100
101
102 Int_t ivol = vol->getMCid();
103 //TLorentzVector tPos1, tMom1;
104 //TLorentzVector tPos, tMom;
105
106
107
108 //#define EDEBUG
109#ifdef EDEBUG
110 static Int_t lEDEBUGcounter=0;
111 if (lEDEBUGcounter<1)
112 std::cout << "EDEBUG-- BmnBd1::ProcessHits: entered" << gMC->CurrentVolPath() << endl;
113#endif
114
115 if(gMC->IsTrackEntering()){
116
117 ResetParameters();
118
119 fELoss = 0.;
120 fIsPrimary = 0;
121 fCharge = -1;
122 fPdgId = 0;
123 lightYield = 0.;
124 lengthtrack = 0.;
125
126 TLorentzVector PosIn;
127 gMC->TrackPosition(PosIn);
128 fPosIn.SetXYZ(PosIn.X(), PosIn.Y(), PosIn.Z());
129
130 TLorentzVector MomIn;
131 gMC->TrackMomentum(MomIn);
132 fMomIn.SetXYZ(MomIn.Px(), MomIn.Py(), MomIn.Pz());
133
134 timeIn = gMC->TrackTime() * 1.0e09;
135
136
137#ifdef EDEBUG
138 //gMC->TrackPosition(tPos1);
139 //gMC->TrackMomentum(tMom1);
140#endif
141
142 TParticle* part = 0;
143 part = gMC->GetStack()->GetCurrentTrack();
144
145 if (part) {
146 fIsPrimary = (Int_t)part->IsPrimary();
147 fCharge = (Int_t)part->GetPDG()->Charge() / 3.;
148 fPdgId = (Int_t)part->GetPdgCode();
149 }
150
151 }
152
153 if (gMC->TrackCharge() != 0) {
154 Double_t BirkC1Mod = 0;
155
156 if (TMath::Abs(gMC->TrackCharge())>=2)
157
158 BirkC1Mod = BirkC1 * 7.2/12.6;
159 else
160 BirkC1Mod = BirkC1;
161 if (gMC->TrackStep() > 0){
162
163 Double_t dedxcm = gMC -> Edep() * 1000./gMC -> TrackStep(); //[MeV/cm]
164 Double_t curLightYield = gMC -> Edep() * 1000./(1. + BirkC1Mod * dedxcm + BirkC2 * dedxcm * dedxcm);
165 curLightYield /= 1000.; //[GeV]
166 lightYield += curLightYield;
167 }
168
169 }
170
171 Double_t eLoss = gMC->Edep();
172 if (eLoss > 0)
173 fELoss += eLoss;
174
175 Double_t slengthtrack = gMC-> TrackStep();
176 lengthtrack += slengthtrack;
177
178
179 if ((gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()) && (fELoss > 0))
180 {
181 TLorentzVector PosOut;
182 gMC->TrackPosition(PosOut);
183 fPosOut.SetXYZ(PosOut.X(), PosOut.Y(), PosOut.Z());
184
185 TLorentzVector MomOut;
186 gMC->TrackMomentum(MomOut);
187 fMomOut.SetXYZ(MomOut.Px(), MomOut.Py(), MomOut.Pz());
188
189 timeOut = gMC->TrackTime() * 1.0e09;
190
191#ifndef EDEBUG
192
193 if (fELoss == 0. ) return kFALSE;
194
195#endif
196
197 fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
198
199 Double_t time = gMC->TrackTime() * 1.0e09;
200 Double_t length = gMC->TrackLength();
201
202 Int_t copyNo;
203 /*Int_t ivol1 = */gMC->CurrentVolID(copyNo);
204
205 Int_t iCell ;
206 gMC->CurrentVolOffID(1, iCell);
207
208
209#ifdef EDEBUG
210 if (lEDEBUGcounter<100) {
211 std::cout << "EDEBUG-- BmnFD::ProcessHits: TrackID:" << fTrackID <<
212 // " ELoss: " << fELoss <<
213 // " particle: " << (part->GetName()) <<
214 " " << gMC->CurrentVolPath() << " " << tPos.Z() <<
215 // " " << (gMC->GetStack()->GetCurrentTrack()->GetMother(1)) <<
216 " " << ivol << "=="<< gMC->CurrentVolID(copyNo) << ","<< copyNo <<
217 " " << gMC->CurrentVolOffID(1,iCell) << " " << iCell <<
218 " " << gMC->CurrentVolOffName(1) << " " << gMC->CurrentVolOffName(0) <<
219 // " " << vol->getRealName() << " " << gMC->CurrentVolPath() <<
220 // " ivol,iCell,copyNo= " << ivol << ","<< iCell << ","<< copyNo <<
221 // " " << vol->getRealName() << " "<< gMC->CurrentVolName() << " "<< gMC->CurrentVolPath() <<
222 // " " << ivol << ","<< vol->getVolumeId() << " : "<< gMC->CurrentVolID(copyNo) << ","<< copyNo <<
223 // " "<< gMC->CurrentVolOffName(2) << " "<< gMC->CurrentVolOffName(3) <<
224 std::endl;
225 lEDEBUGcounter++;
226 }
227
228 AddHit(fTrackID, ivol, copyNo, fPosIn, fPosOut,
229 fMomIn, fMomOut,
230 time, length, fELoss, fIsPrimary, fCharge, fPdgId, lightYield, timeIn, timeOut, lengthtrack);
231#else
232
233 AddHit(fTrackID, ivol, copyNo, fPosIn, fPosOut,
234 fMomIn, fMomOut,
235 time, length, fELoss, fIsPrimary, fCharge, fPdgId, lightYield, timeIn, timeOut, lengthtrack);
236#endif
237
239 // Int_t nBdPoints = (points & (1<<30)) >> 30;
240 // nBdPoints ++;
241 // if (nBdPoints > 1) nBdPoints = 1;
242 // points = ( points & ( ~ (1<<30) ) ) | (nBdPoints << 30);
245
246 ((CbmStack*)gMC->GetStack())->AddPoint(kFD);
247
248 }
249
250 return kTRUE;
251
252#undef EDEBUG
253}
254
255// ----- Public method EndOfEvent -----------------------------------------
256void BmnFD::EndOfEvent() {
257 if (fVerboseLevel) Print("");
258 Reset();
259}
260
261// ----- Public method Register -------------------------------------------
262void BmnFD::Register() {
263 FairRootManager::Instance()->Register("FDPoint","FD", fFDCollection, kTRUE);
264}
265
266// ----- Public method GetCollection --------------------------------------
267TClonesArray* BmnFD::GetCollection(Int_t iColl) const {
268 if (iColl == 0) return fFDCollection;
269
270 return NULL;
271}
272
273// ----- Public method Print ----------------------------------------------
274void BmnFD::Print(Option_t*) const {
275 Int_t nHits = fFDCollection->GetEntriesFast();
276 cout << "-I- BmnFD: " << nHits << " points registered in this event."
277 << endl;
278
279 if (fVerboseLevel>1)
280 for (Int_t i=0; i<nHits; i++) (*fFDCollection)[i]->Print();
281}
282
283
284// ----- Public method Reset ----------------------------------------------
285void BmnFD::Reset() {
286 fFDCollection->Delete();
287
288 fPosIndex = 0;
289}
290
291// guarda in FairRootManager::CopyClones
292// ----- Public method CopyClones -----------------------------------------
293void BmnFD::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset ) {
294 Int_t nEntries = cl1->GetEntriesFast();
295 //cout << "-I- BmnFD: " << nEntries << " entries to add." << endl;
296 TClonesArray& clref = *cl2;
297 BmnFDPoint* oldpoint = NULL;
298 for (Int_t i=0; i<nEntries; i++) {
299 oldpoint = (BmnFDPoint*) cl1->At(i);
300 Int_t index = oldpoint->GetTrackID() + offset;
301 oldpoint->SetTrackID(index);
302 new (clref[fPosIndex]) BmnFDPoint(*oldpoint);
303 fPosIndex++;
304 }
305 cout << " -I- BmnFD: " << cl2->GetEntriesFast() << " merged entries."
306 << endl;
307}
308
309//--------------------------------------------------------------------------------------------------------------------------------------
310void BmnFD::ConstructGeometry()
311{
312 TString fileName = GetGeometryFileName();
313 if(fileName.EndsWith(".root"))
314 {
315 LOG(info) << "Constructing FD geometry from ROOT file " << fileName.Data();
316 ConstructRootGeometry();
317 }
318 else if ( fileName.EndsWith(".geo") )
319 {
320 LOG(info) << "Constructing FD geometry from ASCII file " << fileName.Data();
321 ConstructAsciiGeometry();
322 }
323 else LOG(fatal) << "Geometry format of FD file " << fileName.Data() << " not supported.";
324}
325//--------------------------------------------------------------------------------------------------------------------------------------
326void BmnFD::ConstructAsciiGeometry()
327{
328 FairGeoLoader* geoLoad = FairGeoLoader::Instance();
329 FairGeoInterface* geoFace = geoLoad->getGeoInterface();
330 BmnFDGeo* FDGeo = new BmnFDGeo();
331 FDGeo->setGeomFile(GetGeometryFileName());
332 geoFace->addGeoModule(FDGeo);
333
334 Bool_t rc = geoFace->readSet(FDGeo);
335 if (rc) FDGeo->create(geoLoad->getGeoBuilder());
336 TList* volList = FDGeo->getListOfVolumes();
337
338 // store geo parameter
339 FairRun *fRun = FairRun::Instance();
340 FairRuntimeDb *rtdb = FairRun::Instance()->GetRuntimeDb();
341 BmnFDGeoPar* par = (BmnFDGeoPar*)(rtdb->getContainer("BmnFDGeoPar"));
342 TObjArray *fSensNodes = par->GetGeoSensitiveNodes();
343 TObjArray *fPassNodes = par->GetGeoPassiveNodes();
344
345 TListIter iter(volList);
346 FairGeoNode *node = nullptr;
347 FairGeoVolume *aVol = nullptr;
348
349 while( (node = (FairGeoNode*)iter.Next()) )
350 {
351 aVol = dynamic_cast<FairGeoVolume*>(node);
352
353 if(node->isSensitive()) fSensNodes->AddLast(aVol);
354 else fPassNodes->AddLast(aVol);
355 }
356
357 par->setChanged();
358 par->setInputVersion(fRun->GetRunId(),1);
359 ProcessNodes( volList );
360}
361
362//--------------------------------------------------------------------------------------------------------------------------------------
363
364Bool_t BmnFD::CheckIfSensitive(std::string name)
365{
366 TString tsname = name;
367 if (tsname.Contains("Active")) return kTRUE;
368
369 return kFALSE;
370}
371//-------------------------------------------------------------------------
372
373// ----- Private method AddHit --------------------------------------------
374BmnFDPoint* BmnFD::AddHit( Int_t trackID, Int_t detID, Int_t copyNo, TVector3 posIn, TVector3 posOut,
375 TVector3 momIn, TVector3 momOut, Double_t time, Double_t length, Double_t eLoss, Bool_t isPrimary,
376 Double_t charge, Int_t pdgId, Double_t lightYield, Double_t timeIn, Double_t timeOut, Double_t lengthtrack) {
377
378 TClonesArray& clref = *fFDCollection;
379 Int_t size = clref.GetEntriesFast();
380 return new(clref[size]) BmnFDPoint(trackID, detID, copyNo, posIn, posOut, momIn, momOut, time, length, eLoss, isPrimary,
381 charge, pdgId, lightYield, timeIn, timeOut, lengthtrack);
382}
i
int i
Definition P4_F32vec4.h:22
kFD
@ kFD
Definition BmnDetectorList.h:26
BmnFDGeoPar::GetGeoSensitiveNodes
TObjArray * GetGeoSensitiveNodes()
Definition BmnFDGeoPar.h:20
BmnFDGeoPar::GetGeoPassiveNodes
TObjArray * GetGeoPassiveNodes()
Definition BmnFDGeoPar.h:21
BmnFD::BmnFD
BmnFD()
Definition BmnFD.cxx:33
BmnFD::ConstructGeometry
virtual void ConstructGeometry()
Definition BmnFD.cxx:310
BmnFD::Register
virtual void Register()
Definition BmnFD.cxx:262
BmnFD::CheckIfSensitive
virtual Bool_t CheckIfSensitive(std::string name)
Definition BmnFD.cxx:364
BmnFD::~BmnFD
virtual ~BmnFD()
Definition BmnFD.cxx:61
BmnFD::GetCollection
virtual TClonesArray * GetCollection(Int_t iColl) const
Definition BmnFD.cxx:267
BmnFD::BeginEvent
virtual void BeginEvent()
Definition BmnFD.cxx:79
BmnFD::Print
virtual void Print(Option_t *) const
Definition BmnFD.cxx:274
BmnFD::ConstructAsciiGeometry
virtual void ConstructAsciiGeometry()
Definition BmnFD.cxx:326
BmnFD::Reset
virtual void Reset()
Definition BmnFD.cxx:285
BmnFD::Initialize
virtual void Initialize()
Definition BmnFD.cxx:71
BmnFD::EndOfEvent
virtual void EndOfEvent()
Definition BmnFD.cxx:256
BmnFD::CopyClones
virtual void CopyClones(TClonesArray *cl1, TClonesArray *cl2, Int_t offset)
Definition BmnFD.cxx:293
BmnFD::AddHit
BmnFDPoint * AddHit(Int_t trackID, Int_t detID, Int_t copyNo, TVector3 posIn, TVector3 posOut, TVector3 momIn, TVector3 momOut, Double_t tof, Double_t length, Double_t eLoss, Bool_t isPrimary, Double_t charge, Int_t pdgId, Double_t lightYield, Double_t timeIn, Double_t timeOut, Double_t lengthtrack)
Definition BmnFD.cxx:374
BmnFD::ProcessHits
virtual Bool_t ProcessHits(FairVolume *vol=0)
Definition BmnFD.cxx:87