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