ATTPCROOT  0.3.0-alpha
A ROOT-based framework for analyzing data from active target detectors
AtSiArray.cxx
Go to the documentation of this file.
1 #include "AtSiArray.h"
2 
3 #include "AtDetectorList.h"
4 #include "AtSiPoint.h"
5 #include "AtStack.h"
6 
7 #include <FairDetector.h>
8 #include <FairLogger.h>
9 #include <FairRootManager.h>
10 #include <FairVolume.h>
11 
12 #include <TClonesArray.h>
13 #include <TGeoManager.h>
14 #include <TLorentzVector.h>
15 #include <TVector3.h>
16 #include <TVirtualMC.h>
17 #include <TVirtualMCStack.h>
18 
19 #include <iostream>
20 
21 using std::cout;
22 using std::endl;
23 
24 AtSiArray::AtSiArray()
25  : FairDetector("AtSiArray", kTRUE, kAtSiArray), fTrackID(-1), fVolumeID(-1), fPos(), fMom(), fTime(-1.),
26  fLength(-1.), fELoss(-1), fPosIndex(-1), fAtSiArrayPointCollection(new TClonesArray("AtSiPoint")), fELossAcc(-1)
27 {
28  // LOG(INFO)<<" AtSiArray detector initialized ";
29 }
30 
31 AtSiArray::AtSiArray(const char *name, Bool_t active)
32  : FairDetector(name, active, kAtSiArray), fTrackID(-1), fVolumeID(-1), fPos(), fMom(), fTime(-1.), fLength(-1.),
33  fELoss(-1), fPosIndex(-1), fAtSiArrayPointCollection(new TClonesArray("AtSiPoint")), fELossAcc(-1)
34 {
35  // LOG(INFO)<<" AtSiArray detector initialized ";
36 }
37 
38 AtSiArray::~AtSiArray()
39 {
40  if (fAtSiArrayPointCollection) {
41  fAtSiArrayPointCollection->Delete();
42  delete fAtSiArrayPointCollection;
43  }
44 }
45 
46 void AtSiArray::Initialize()
47 {
48  FairDetector::Initialize();
49  // FairRuntimeDb *rtdb = FairRun::Instance()->GetRuntimeDb();
50  // auto *par = (AtSiArrayGeoPar *)(rtdb->getContainer("AtSiArrayGeoPar"));
51 }
52 
53 Bool_t AtSiArray::ProcessHits(FairVolume *vol)
54 {
57  auto *stack = dynamic_cast<AtStack *>(TVirtualMC::GetMC()->GetStack());
58  std::pair<Int_t, Int_t> AZ;
59  AZ = DecodePdG(gMC->TrackPid());
60  fVolName = gMC->CurrentVolName();
61  Int_t VolumeID = 0;
62 
63  LOG(debug) << "In AtSiArray::ProcessHits";
64  // Set parameters at entrance of volume. Reset ELoss.
65  if (TVirtualMC::GetMC()->IsTrackEntering()) {
66 
67  fELoss = 0.;
68  fTime = TVirtualMC::GetMC()->TrackTime() * 1.0e09;
69  fLength = TVirtualMC::GetMC()->TrackLength();
70  TVirtualMC::GetMC()->TrackPosition(fPosIn);
71  TVirtualMC::GetMC()->TrackMomentum(fMomIn);
72 
73  std::cout << " AtSiArray: Track is entering "
74  << "\n";
75  LOG(INFO) << " HELIOS: First hit in Volume " << fVolName;
76  LOG(INFO) << " Particle : " << gMC->ParticleName(gMC->TrackPid());
77  LOG(INFO) << " PID PdG : " << gMC->TrackPid();
78  LOG(INFO) << " Atomic Mass : " << AZ.first;
79  LOG(INFO) << " Atomic Number : " << AZ.second;
80  LOG(INFO) << " Volume ID " << gMC->CurrentVolID(VolumeID);
81  LOG(INFO) << " Track ID : " << fTrackID;
82  LOG(INFO) << " Position In : " << fPosIn.X() << " " << fPosIn.Y() << " " << fPosIn.Z() << std::endl;
83  LOG(INFO) << " Position Out : " << fPosOut.X() << " " << fPosOut.Y() << " " << fPosOut.Z() << std::endl;
84  LOG(INFO) << " Momentum In: " << fMomIn.X() << " " << fMomIn.Y() << " " << fMomIn.Z() << std::endl;
85  // LOG(INFO)<<" Total relativistic energy " <<gMC->Etot()<< FairLogger::endl;
86  // LOG(INFO)<<" Mass of the Tracked particle (gAVTP) :
87  // "<<AtVertexPropagator::Instance()->GetBeamMass()<<std::endl; LOG(INFO)<<" Mass of the Tracked particle (gMC) :
88  // "<<gMC->TrackMass()<<std::endl; LOG(INFO)<<" Initial energy of the current particle in this volume :
89  // "<<((gMC->Etot() - gMC->TrackMass()) * 1000.)<<FairLogger::endl;// Relativistic Mass
90  }
91 
92  // Sum energy loss for all steps in the active volume
93  fELoss = TVirtualMC::GetMC()->Edep();
94  // std::cout<<" Energy loss : "<<fELoss*1000<<"\n";
95  gMC->TrackPosition(fPosIn);
96  gMC->TrackMomentum(fMomIn);
97  fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
98  fVolumeID = vol->getMCid();
99  fTime = gMC->TrackTime() * 1.0e09;
100  fLength = gMC->TrackLength();
101 
102  // Create AtSiArrayPoint at exit of active volume
103  if (TVirtualMC::GetMC()->IsTrackExiting() || TVirtualMC::GetMC()->IsTrackStop() ||
104  TVirtualMC::GetMC()->IsTrackDisappeared()) {
105  fTrackID = TVirtualMC::GetMC()->GetStack()->GetCurrentTrackNumber();
106  fVolumeID = vol->getMCid();
107  TVirtualMC::GetMC()->TrackPosition(fPosOut);
108  TVirtualMC::GetMC()->TrackMomentum(fMomOut);
109 
110  /*if (fELoss == 0.) {
111  return kFALSE;
112  }*/
113 
114  if (gMC->IsTrackExiting()) {
115 
116  const Double_t *oldpos = nullptr;
117  const Double_t *olddirection = nullptr;
118  Double_t newpos[3];
119  Double_t newdirection[3];
120  Double_t safety = 0;
121 
122  gGeoManager->FindNode(fPosOut.X(), fPosOut.Y(), fPosOut.Z());
123  oldpos = gGeoManager->GetCurrentPoint();
124  olddirection = gGeoManager->GetCurrentDirection();
125 
126  for (Int_t i = 0; i < 3; i++) {
127  newdirection[i] = -1 * olddirection[i];
128  }
129 
130  gGeoManager->SetCurrentDirection(newdirection);
131  // TGeoNode *bla = gGeoManager->FindNextBoundary(2);
132  safety = gGeoManager->GetSafeDistance();
133 
134  gGeoManager->SetCurrentDirection(-newdirection[0], -newdirection[1], -newdirection[2]);
135 
136  for (Int_t i = 0; i < 3; i++) {
137  newpos[i] = oldpos[i] - (3 * safety * olddirection[i]);
138  }
139 
140  fPosOut.SetX(newpos[0]);
141  fPosOut.SetY(newpos[1]);
142  fPosOut.SetZ(newpos[2]);
143 
144  std::cout << " AtSiArray: Track is exiting "
145  << "\n";
146  LOG(INFO) << " HELIOS: First hit in Volume " << fVolName;
147  LOG(INFO) << " Particle : " << gMC->ParticleName(gMC->TrackPid());
148  LOG(INFO) << " PID PdG : " << gMC->TrackPid();
149  LOG(INFO) << " Atomic Mass : " << AZ.first;
150  LOG(INFO) << " Atomic Number : " << AZ.second;
151  LOG(INFO) << " Volume ID " << gMC->CurrentVolID(VolumeID);
152  LOG(INFO) << " Track ID : " << fTrackID;
153  LOG(INFO) << " Position In : " << fPosIn.X() << " " << fPosIn.Y() << " " << fPosIn.Z() << std::endl;
154  LOG(INFO) << " Position Out : " << fPosOut.X() << " " << fPosOut.Y() << " " << fPosOut.Z() << std::endl;
155  LOG(INFO) << " Momentum In: " << fMomIn.X() << " " << fMomIn.Y() << " " << fMomIn.Z() << std::endl;
156  LOG(INFO) << " Momentum Out: " << fMomOut.X() << " " << fMomOut.Y() << " " << fMomOut.Z() << std::endl;
157  }
158 
159  /* AddHit(fTrackID,
160  fVolumeID,
161  fVolName,
162  fDetCopyID,
163  TVector3(fPosIn.X(), fPosIn.Y(), fPosIn.Z()),
164  TVector3(fPosOut.X(), fPosOut.Y(), fPosOut.Z()),
165  TVector3(fMomIn.Px(), fMomIn.Py(), fMomIn.Pz()),
166  TVector3(fMomOut.Px(), fMomOut.Py(), fMomOut.Pz()),
167  fTime,
168  fLength,
169  fELoss,
170  0.0,
171  0.0,
172  AZ.first,
173  AZ.second);
174 
175 
176 
177 
178  AtStack* stack = static_cast<AtStack*>(TVirtualMC::GetMC()->GetStack());
179  stack->AddPoint(kAtSiArray);*/
180  }
181 
182  AddHit(fTrackID, fVolumeID, fVolName, fDetCopyID, TVector3(fPosIn.X(), fPosIn.Y(), fPosIn.Z()),
183  TVector3(fPosOut.X(), fPosOut.Y(), fPosOut.Z()), TVector3(fMomIn.Px(), fMomIn.Py(), fMomIn.Pz()),
184  TVector3(fMomOut.Px(), fMomOut.Py(), fMomOut.Pz()), fTime, fLength, fELoss, 0.0, 0.0, AZ.first, AZ.second);
185 
186  stack->AddPoint(kAtSiArray);
187 
188  return kTRUE;
189 }
190 
191 void AtSiArray::EndOfEvent()
192 {
193 
194  fAtSiArrayPointCollection->Clear();
195 }
196 
197 void AtSiArray::Register()
198 {
199 
206  FairRootManager::Instance()->Register("AtSiArrayPoint", "AtSiArray", fAtSiArrayPointCollection, kTRUE);
207 }
208 
209 TClonesArray *AtSiArray::GetCollection(Int_t iColl) const
210 {
211  if (iColl == 0) {
212  return fAtSiArrayPointCollection;
213  } else {
214  return nullptr;
215  }
216 }
217 
218 void AtSiArray::Reset()
219 {
220  fAtSiArrayPointCollection->Clear();
221 }
222 
223 void AtSiArray::Print(Option_t *option) const
224 {
225  Int_t nHits = fAtSiArrayPointCollection->GetEntriesFast();
226  LOG(INFO) << "Si Array: " << nHits << " points registered in this event";
227 }
228 
229 void AtSiArray::ConstructGeometry()
230 {
231  TString fileName = GetGeometryFileName();
232  if (fileName.EndsWith(".geo")) {
233  LOG(INFO) << "Constructing Si Array geometry from ASCII file " << fileName;
234  // ConstructASCIIGeometry();
235  } else if (fileName.EndsWith(".root")) {
236  LOG(INFO) << "Constructing Si Array geometry from ROOT file " << fileName;
237  ConstructRootGeometry();
238  } else {
239  std::cout << "Geometry format not supported." << std::endl;
240  }
241 }
242 
243 Bool_t AtSiArray::CheckIfSensitive(std::string name)
244 {
245 
246  TString tsname = name;
247  if (tsname.Contains("silicon")) {
248  LOG(INFO) << " Si Array geometry: Sensitive volume found: " << tsname;
249  return kTRUE;
250  }
251  return kFALSE;
252 }
253 
254 AtSiPoint *AtSiArray::AddHit(Int_t trackID, Int_t detID, TVector3 pos, TVector3 mom, Double_t time, Double_t length,
255  Double_t eLoss)
256 {
257  TClonesArray &clref = *fAtSiArrayPointCollection;
258  Int_t size = clref.GetEntriesFast();
259  return new (clref[size]) AtSiPoint(trackID, detID, pos, mom, time, length, eLoss);
260 }
261 
262 // ----- Private method AddHit --------------------------------------------
263 AtSiPoint *AtSiArray::AddHit(Int_t trackID, Int_t detID, TString VolName, Int_t detCopyID, TVector3 posIn,
264  TVector3 posOut, TVector3 momIn, TVector3 momOut, Double_t time, Double_t length,
265  Double_t eLoss, Double_t EIni, Double_t AIni, Int_t A, Int_t Z)
266 {
267  TClonesArray &clref = *fAtSiArrayPointCollection;
268  Int_t size = clref.GetEntriesFast();
269  // std::cout<< "AtTPC: Adding Point at (" << posIn.X() << ", " << posIn.Y() << ", " << posIn.Z() << ") cm, detector
270  // " << detID << ", track " << trackID
271  //<< ", energy loss " << eLoss << " MeV" <<" with accumulated Energy Loss : "<<fELossAcc<<" MeV "<< std::endl;
272  if (fVerboseLevel > 1)
273  LOG(INFO) << "Si Array: Adding Point at (" << posIn.X() << ", " << posIn.Y() << ", " << posIn.Z()
274  << ") cm, detector " << detID << ", track " << trackID << ", energy loss " << eLoss * 1e06 << " keV";
275 
276  return new (clref[size]) AtSiPoint(trackID, detID, posIn, posOut, momIn, momOut, time, length, eLoss, VolName,
277  detCopyID, EIni, AIni, A, Z);
278 }
279 
280 std::pair<Int_t, Int_t> AtSiArray::DecodePdG(Int_t PdG_Code)
281 {
282  Int_t A = PdG_Code / 10 % 1000;
283  Int_t Z = PdG_Code / 10000 % 1000;
284 
285  std::pair<Int_t, Int_t> nucleus;
286 
287  if (PdG_Code == 2212) {
288  nucleus.first = 1;
289  nucleus.second = 1;
290  } else if (PdG_Code == 2112) {
291  nucleus.first = 1;
292  nucleus.second = 0;
293  } else {
294  nucleus.first = A;
295  nucleus.second = Z;
296  }
297 
298  return nucleus;
299 }
300 
301 ClassImp(AtSiArray);
AtSiArray::CheckIfSensitive
Bool_t CheckIfSensitive(std::string name)
Definition: AtSiArray.cxx:243
AtDetectorList.h
AtSiArray::AtSiArray
AtSiArray()
Definition: AtSiArray.cxx:24
AtSiArray
Definition: AtSiArray.h:22
AtSiArray.h
AtSiArray::Reset
virtual void Reset()
Definition: AtSiArray.cxx:218
AtSiArray::EndOfEvent
virtual void EndOfEvent()
Definition: AtSiArray.cxx:191
ClassImp
ClassImp(AtSiArray)
AtSiArray::ProcessHits
virtual Bool_t ProcessHits(FairVolume *v=0)
Definition: AtSiArray.cxx:53
AtSiArray::ConstructGeometry
void ConstructGeometry()
Definition: AtSiArray.cxx:229
AtSiPoint
Definition: AtSiPoint.h:14
AtSiArray::Initialize
virtual void Initialize()
Definition: AtSiArray.cxx:46
AtSiArray::AddHit
AtSiPoint * AddHit(Int_t trackID, Int_t detID, TVector3 pos, TVector3 mom, Double_t time, Double_t length, Double_t eLoss)
Definition: AtSiArray.cxx:254
AtSiArray::Print
virtual void Print(Option_t *option="") const
Definition: AtSiArray.cxx:223
AtStack.h
AtSiArray::~AtSiArray
virtual ~AtSiArray()
Definition: AtSiArray.cxx:38
AtSiArray::Register
virtual void Register()
Definition: AtSiArray.cxx:197
AtSiArray::GetCollection
virtual TClonesArray * GetCollection(Int_t iColl) const
Definition: AtSiArray.cxx:209
kAtSiArray
@ kAtSiArray
Definition: AtDetectorList.h:28
AtStack
Definition: AtStack.h:54
AtSiPoint.h
AtSiArray::DecodePdG
std::pair< Int_t, Int_t > DecodePdG(Int_t PdG_Code)
Definition: AtSiArray.cxx:280
Generated on Sat Nov 11 2023 15:30:38 for ATTPCROOT by   doxygen 1.8.18