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