ATTPCROOT  0.3.0-alpha
A ROOT-based framework for analyzing data from active target detectors
AtTPC2Body.cxx
Go to the documentation of this file.
1 #include "AtTPC2Body.h"
2 
3 #include "AtEulerTransformation.h"
4 #include "AtVertexPropagator.h"
5 
6 #include <FairIon.h>
7 #include <FairLogger.h>
8 #include <FairParticle.h>
9 #include <FairPrimaryGenerator.h>
10 #include <FairRunSim.h>
11 
12 #include <TDatabasePDG.h>
13 #include <TMath.h>
14 #include <TParticle.h>
15 #include <TParticlePDG.h>
16 #include <TRandom.h>
17 #include <TVector3.h>
18 
19 #include <algorithm>
20 #include <cmath>
21 #include <cstdio>
22 #include <iostream>
23 #include <memory>
24 
25 constexpr auto cRED = "\033[1;31m";
26 constexpr auto cYELLOW = "\033[1;33m";
27 constexpr auto cNORMAL = "\033[0m";
28 constexpr auto cGREEN = "\033[1;32m";
29 constexpr auto cBLUE = "\033[1;34m";
30 
31 Int_t AtTPC2Body::fgNIon = 0;
32 
33 AtTPC2Body::AtTPC2Body() : fMult(0), fPx(0.), fPy(0.), fPz(0.), fVx(0.), fVy(0.), fVz(0.), fIon(0), fQ(0)
34 {
35  // cout << "-W- AtTPCIonGenerator: "
36  // << " Please do not use the default constructor! " << endl;
37 }
38 
39 // ----- Default constructor ------------------------------------------
40 AtTPC2Body::AtTPC2Body(const char *name, std::vector<Int_t> *z, std::vector<Int_t> *a, std::vector<Int_t> *q,
41  Int_t mult, std::vector<Double_t> *px, std::vector<Double_t> *py, std::vector<Double_t> *pz,
42  std::vector<Double_t> *mass, std::vector<Double_t> *Ex, Double_t ResEner, Double_t MinCMSAng,
43  Double_t MaxCMSAng)
44  : fMult(mult), fPx(0.), fPy(0.), fPz(0.), fVx(0.), fVy(0.), fVz(0.), fIon(0), fPType(0.), fQ(0),
45  fThetaCmsMin(MinCMSAng), fThetaCmsMax(MaxCMSAng), fBeamEnergy_buff(ResEner), fNoSolution(kFALSE),
46  fIsFixedTargetPos(kFALSE), fIsFixedMomentum(kFALSE)
47 {
48  fgNIon++;
49  fIon.reserve(fMult);
50 
51  char buffer[30];
52  auto *kProton = new TParticle(); // NOLINT but probably a problem
53  kProton->SetPdgCode(2212);
54  auto *kNeutron = new TParticle(); // NOLINT but probably a problem
55  kNeutron->SetPdgCode(2112);
56 
57  for (Int_t i = 0; i < fMult; i++) {
58 
59  fPx.push_back(Double_t(a->at(i)) * px->at(i));
60  fPy.push_back(Double_t(a->at(i)) * py->at(i));
61  fPz.push_back(Double_t(a->at(i)) * pz->at(i));
62  Masses.push_back(mass->at(i) * 1000.0);
63  fExEnergy.push_back(Ex->at(i));
64  fWm.push_back(mass->at(i) * 1000.0 * 0.93149401 + Ex->at(i));
65  FairIon *IonBuff;
66  FairParticle *ParticleBuff;
67  sprintf(buffer, "Product_Ion%d", i);
68 
69  if (a->at(i) != 1) {
70 
71  IonBuff = new FairIon(buffer, z->at(i), a->at(i), q->at(i), 0.0, mass->at(i)); // NOLINT but probably a problem
72  ParticleBuff = new FairParticle("dummyPart", 1, 1, 1.0, 0, 0.0, 0.0); // NOLINT but probably a problem
73  fPType.emplace_back("Ion");
74  std::cout << " Adding : " << buffer << std::endl;
75 
76  } else if (a->at(i) == 1 && z->at(i) == 1) {
77 
78  IonBuff = new FairIon("dummyIon", 50, 50, 0, 0.0, 100); // NOLINT but probably a problem
79  ParticleBuff = new FairParticle(2212, kProton); // NOLINT but probably a problem
80  fPType.emplace_back("Proton");
81 
82  } else if (a->at(i) == 1 && z->at(i) == 0) {
83 
84  IonBuff = new FairIon("dummyIon", 50, 50, 0, 0.0, 100); // NOLINT but probably a problem
85  ParticleBuff = new FairParticle(2112, kNeutron); // NOLINT but probably a problem
86  fPType.emplace_back("Neutron");
87  }
88 
89  std::cout << " Z " << z->at(i) << " A " << a->at(i) << std::endl;
90  // std::cout<<buffer<<std::endl;
91  fIon.push_back(IonBuff);
92  fParticle.push_back(ParticleBuff);
93  }
94 
95  FairRunSim *run = FairRunSim::Instance();
96  if (!run) {
97  std::cout << "-E- FairIonGenerator: No FairRun instantised!" << std::endl;
98  Fatal("FairIonGenerator", "No FairRun instantised!");
99  }
100 
101  for (Int_t i = 0; i < fMult; i++) {
102 
103  if (fPType.at(i) == "Ion") {
104 
105  std::cout << " In position " << i << " adding an : " << fPType.at(i) << std::endl;
106  run->AddNewIon(fIon.at(i)); // NOLINT
107  std::cout << " fIon name :" << fIon.at(i)->GetName() << std::endl;
108  std::cout << " fParticle name :" << fParticle.at(i)->GetName() << std::endl;
109 
110  } else if (fPType.at(i) == "Proton") {
111 
112  std::cout << " In position " << i << " adding an : " << fPType.at(i) << std::endl;
113  // run->AddNewParticle(fParticle.at(i));
114  std::cout << " fIon name :" << fIon.at(i)->GetName() << std::endl;
115  std::cout << " fParticle name :" << fParticle.at(i)->GetName() << std::endl;
116  std::cout << fParticle.at(i)->GetName() << std::endl;
117 
118  } else if (fPType.at(i) == "Neutron") {
119 
120  std::cout << " In position " << i << " adding an : " << fPType.at(i) << std::endl;
121  // run->AddNewParticle(fParticle.at(i));
122  std::cout << " fIon name :" << fIon.at(i)->GetName() << std::endl;
123  std::cout << " fParticle name :" << fParticle.at(i)->GetName() << std::endl;
124  std::cout << fParticle.at(i)->GetName() << std::endl;
125  }
126  }
127 }
128 
129 void AtTPC2Body::SetFixedTargetPosition(double vx, double vy, double vz)
130 {
131  std::cout << " -I- AtTPC2Body : Fixed target position at " << vx << " " << vy << " " << vz << std::endl;
132  fIsFixedTargetPos = kTRUE;
133  fVx = vx;
134  fVy = vy;
135  fVz = vz;
136 }
137 
138 void AtTPC2Body::SetFixedBeamMomentum(double px, double py, double pz)
139 {
140  std::cout << " -I- AtTPC2Body : Fixed beam momentum " << px << " " << py << " " << pz << std::endl;
141  fIsFixedMomentum = kTRUE;
142  fPxBeam_buff = px;
143  fPyBeam_buff = py;
144  fPzBeam_buff = pz;
145 }
146 
147 // ----- Public method ReadEvent --------------------------------------
148 Bool_t AtTPC2Body::ReadEvent(FairPrimaryGenerator *primGen)
149 {
150 
151  std::vector<Double_t> Ang; // Lab Angle of the products
152  std::vector<Double_t> Ene; // Lab Energy of the products
153  Ang.reserve(2);
154  Ang.reserve(2);
155  fPx.clear();
156  fPy.clear();
157  fPx.clear();
158 
159  fPx.resize(fMult);
160  fPy.resize(fMult);
161  fPx.resize(fMult);
162 
163  Double_t costhetamin = TMath::Cos(fThetaCmsMin * TMath::DegToRad());
164  Double_t costhetamax = TMath::Cos(fThetaCmsMax * TMath::DegToRad());
165  // Double_t thetacmsInput = fThetaCmsMin + ((fThetaCmsMax-fThetaCmsMin)*gRandom->Uniform());
167  Double_t thetacmsInput =
168  TMath::ACos((costhetamax - costhetamin) * gRandom->Uniform() + costhetamin) * TMath::RadToDeg();
169 
170  std::cout << cBLUE << " -I- AtTPC2Body : Random CMS Theta angle in degrees : " << thetacmsInput << cNORMAL
171  << std::endl;
172  const Double_t rad2deg = 0.0174532925;
173 
174  if (fIsFixedTargetPos) {
175  fBeamEnergy = fBeamEnergy_buff;
176  AtVertexPropagator::Instance()->SetValidKine(kTRUE);
177  std::cout << cBLUE << " -I- AtTPC2Body Beam energy (Fixed Target mode) : " << fBeamEnergy << cNORMAL << std::endl;
178 
179  } else {
180  fBeamEnergy = AtVertexPropagator::Instance()->GetEnergy();
181 
182  std::cout << cBLUE << " -I- AtTPC2Body Residual energy (Active Target mode) : "
183  << AtVertexPropagator::Instance()->GetEnergy() << cNORMAL << std::endl;
184  }
185 
186  if (fBeamEnergy > 0 &&
187  (AtVertexPropagator::Instance()->GetDecayEvtCnt() % 2 != 0 ||
188  fIsFixedTargetPos)) { // Requires a non zero vertex energy and pre-generated Beam event (not punch thorugh)
189 
190  if (fIsFixedTargetPos) {
191  fPxBeam = fPxBeam_buff;
192  fPyBeam = fPyBeam_buff;
193  fPzBeam = fPzBeam_buff;
194 
195  } else {
196 
197  fPxBeam = AtVertexPropagator::Instance()->GetPx();
198  fPyBeam = AtVertexPropagator::Instance()->GetPy();
199  fPzBeam = AtVertexPropagator::Instance()->GetPz();
200  }
201 
202  Double_t eb = fBeamEnergy + fWm.at(0);
203  Double_t pb2 = fBeamEnergy * fBeamEnergy + 2.0 * fBeamEnergy * fWm.at(0);
204  Double_t pb = TMath::Sqrt(pb2);
205  Double_t beta = pb / (eb + fWm.at(1));
206  Double_t gamma = 1.0 / sqrt(1.0 - beta * beta);
207 
208  Double_t thetacms = thetacmsInput * rad2deg; // degree to radian
209 
210  Double_t thetacmr = TMath::Pi() - thetacms;
211  Double_t e = fBeamEnergy + fWm.at(0) + fWm.at(1);
212  Double_t e_cm2 = e * e - pb2;
213  Double_t e_cm = TMath::Sqrt(e_cm2);
214  Double_t t_cm = e_cm - fWm.at(2) - fWm.at(3);
215 
216  if (t_cm < 0.0) {
217  std::cout << "-I- AtTPC2Body : No solution!" << std::endl;
218  fNoSolution = kTRUE;
219  AtVertexPropagator::Instance()->SetValidKine(kFALSE);
220  // return kFALSE;
221  }
222 
223  if (AtVertexPropagator::Instance()->GetValidKine()) {
224 
225  Double_t t_cm2 = t_cm * t_cm;
226  Double_t t3_cm = (t_cm2 + 2. * fWm.at(3) * t_cm) / (t_cm + fWm.at(2) + fWm.at(3)) / 2.0;
227  Double_t t4_cm = (t_cm2 + 2. * fWm.at(2) * t_cm) / (t_cm + fWm.at(2) + fWm.at(3)) / 2.0;
228  Double_t p3_cm2 = t3_cm * t3_cm + 2.0 * t3_cm * fWm.at(2);
229  Double_t p3_cm = TMath::Sqrt(p3_cm2);
230  Double_t tg_thetalabs =
231  p3_cm * TMath::Sin(thetacms) /
232  (gamma * (p3_cm * TMath::Cos(thetacms) + beta * TMath::Sqrt(p3_cm * p3_cm + fWm.at(2) * fWm.at(2))));
233 
234  if (tg_thetalabs >= 1.0e6) {
235  Ang.push_back(TMath::Pi() / 2.0);
236  } else {
237  Ang.push_back(TMath::ATan(tg_thetalabs));
238  }
239 
240  if (Ang.at(0) < 0.0)
241  Ang.at(0) = TMath::Pi() + Ang.at(0);
242 
243  Double_t p4_cm2 = t4_cm * t4_cm + 2. * t4_cm * fWm.at(3);
244  Double_t p4_cm = TMath::Sqrt(p4_cm2);
245  Double_t tg_thetalabr =
246  p4_cm * TMath::Sin(thetacmr) /
247  (gamma * (p4_cm * TMath::Cos(thetacmr) + beta * TMath::Sqrt(p4_cm * p4_cm + fWm.at(3) * fWm.at(3))));
248  // std::cout<<" tg_thetalabr : "<<tg_thetalabr<<std::endl;
249 
250  if (tg_thetalabr > 1.0e6) {
251  Ang.push_back(TMath::Pi() / 2.0);
252  } else {
253  Ang.push_back(TMath::ATan(tg_thetalabr));
254  }
255 
256  if (Ang.at(1) < 0.0)
257  Ang.at(1) = TMath::Pi() + Ang.at(1);
258 
259  // Lorentz transformations to lab -----
260 
261  Double_t p3_cmx = p3_cm * sin(thetacms);
262  Double_t p3_cmz = p3_cm * cos(thetacms);
263  Double_t p3_labx = p3_cmx;
264  Double_t p3_labz = gamma * (p3_cmz + beta * (t3_cm + fWm.at(2)));
265  Double_t p3_lab = TMath::Sqrt(p3_labx * p3_labx + p3_labz * p3_labz);
266  Ene.push_back(TMath::Sqrt(p3_lab * p3_lab + fWm.at(2) * fWm.at(2)) - fWm.at(2));
267 
268  Double_t p4_cmx = p4_cm * sin(thetacmr);
269  Double_t p4_cmz = p4_cm * cos(thetacmr);
270  Double_t p4_labx = p4_cmx;
271  Double_t p4_labz = gamma * (p4_cmz + beta * (t4_cm + fWm.at(3)));
272  Double_t p4_lab = TMath::Sqrt(p4_labx * p4_labx + p4_labz * p4_labz);
273  Ene.push_back(TMath::Sqrt(p4_lab * p4_lab + fWm.at(3) * fWm.at(3)) - fWm.at(3));
274 
275  if (!AtVertexPropagator::Instance()->GetValidKine()) {
276  std::cout << " -I- ===== AtTPC2Body - Kinematics ====== " << std::endl;
277  std::cout << " -I- ===== No Valid Solution on Beam Event for these kinematics (probably due to threshold "
278  "energy) ====== "
279  << std::endl;
280  std::cout << " -I- ===== Setting Values to 0 ====== " << std::endl;
281  Ene.at(0) = 0.0;
282  Ang.at(0) = 0.0;
283  Ene.at(1) = 0.0;
284  Ang.at(1) = 0.0;
285 
286  } else {
287 
288  std::cout << cBLUE << " -I- ===== AtTPC2Body - Kinematics ====== " << std::endl;
289  std::cout << " Decay of scattered particle enabled : " << kIsDecay << "\n";
290  std::cout << " Scattered energy:" << Ene.at(0) << " MeV" << std::endl;
291  std::cout << " Scattered angle:" << Ang.at(0) * 180 / TMath::Pi() << " deg" << std::endl;
292  std::cout << " Recoil energy:" << Ene.at(1) << " MeV" << std::endl;
293  std::cout << " Recoiled angle:" << Ang.at(1) * 180.0 / TMath::Pi() << " deg" << cNORMAL << std::endl;
294  }
295 
296  // AtVertexPropagator::Instance()->SetRecoilE(Ene.at(1));
297  AtVertexPropagator::Instance()->SetTrackEnergy(1, Ene.at(1));
298  AtVertexPropagator::Instance()->SetTrackAngle(1, Ang.at(1) * 180.0 / TMath::Pi());
299 
300  AtVertexPropagator::Instance()->SetTrackEnergy(0, Ene.at(0));
301  AtVertexPropagator::Instance()->SetTrackAngle(0, Ang.at(0) * 180.0 / TMath::Pi());
302 
303  fPx.at(0) = 0.0;
304  fPy.at(0) = 0.0;
305  fPz.at(0) = 0.0;
306 
307  fPx.at(1) = 0.0;
308  fPy.at(1) = 0.0;
309  fPz.at(1) = 0.0;
310 
311  fPx.at(2) = p3_labx / 1000.0; // To GeV for FairRoot
312  fPy.at(2) = 0.0;
313  fPz.at(2) = p3_labz / 1000.0;
314 
315  fPx.at(3) = p4_labx / 1000.0;
316  fPy.at(3) = 0.0;
317  fPz.at(3) = p4_labz / 1000.0;
318 
319  Double_t phiBeam1 = 0., phiBeam2 = 0.;
320 
321  phiBeam1 = 2 * TMath::Pi() * gRandom->Uniform(); // flat probability in phi
322  phiBeam2 = phiBeam1 + TMath::Pi();
323 
324  // std::cout<<" Propagated Entrance Position 2 - X : "<<AtVertexPropagator::Instance()->GetVx()<<" - Y :
325  // "<<AtVertexPropagator::Instance()->GetVy()<<" - Z :
326  // "<<AtVertexPropagator::Instance()->GetVz()<<std::endl; std::cout<<" Propagated Stop Position - X :
327  // "<<AtVertexPropagator::Instance()->GetInVx()<<" - Y :
328  // "<<AtVertexPropagator::Instance()->GetInVy()<<" - Z :
329  // "<<AtVertexPropagator::Instance()->GetInVz()<<std::endl;
330 
331  /* TVector3 BeamPos( AtVertexPropagator::Instance()->GetVx() - AtVertexPropagator::Instance()->GetInVx() ,
332  AtVertexPropagator::Instance()->GetVy() - AtVertexPropagator::Instance()->GetInVy() ,
333  AtVertexPropagator::Instance()->GetVz() - AtVertexPropagator::Instance()->GetInVz() ); std::cout << " Beam
334  Theta (Pos) : "<<BeamPos.Theta()*180.0/TMath::Pi()<<std::endl; std::cout << " Beam Phi (Pos) :
335  "<<BeamPos.Phi()*180.0/TMath::Pi()<<std::endl;*/
336 
337  TVector3 BeamPos(fPxBeam * 1000, fPyBeam * 1000, fPzBeam * 1000); // To MeV for Euler Transformation
338  // TVector3 BeamPos(1.0,1.0,0.0);
339  LOG(DEBUG) << " Beam Theta (Mom) : " << BeamPos.Theta() * 180.0 / TMath::Pi();
340  LOG(DEBUG) << " Beam Phi (Mom) : " << BeamPos.Phi() * 180.0 / TMath::Pi();
341 
342  Double_t thetaLab1, phiLab1, thetaLab2, phiLab2;
343  auto EulerTransformer = std::make_unique<AtEulerTransformation>();
344  EulerTransformer->SetBeamDirectionAtVertexTheta(BeamPos.Theta());
345  EulerTransformer->SetBeamDirectionAtVertexPhi(BeamPos.Phi());
346 
347  EulerTransformer->SetThetaInBeamSystem(Ang.at(0));
348  EulerTransformer->SetPhiInBeamSystem(phiBeam1);
349  EulerTransformer->DoTheEulerTransformationBeam2Lab(); // Euler transformation for particle 1
350  // EulerTransformer->PrintResults();
351 
352  thetaLab1 = EulerTransformer->GetThetaInLabSystem();
353  phiLab1 = EulerTransformer->GetPhiInLabSystem();
354  LOG(DEBUG) << " Scattered angle Phi :" << phiBeam1 * 180.0 / TMath::Pi() << " deg";
355  LOG(DEBUG) << " Scattered angle Theta (Euler) :" << thetaLab1 * 180.0 / TMath::Pi() << " deg";
356  LOG(DEBUG) << " Scattered angle Phi (Euler) :" << phiLab1 * 180.0 / TMath::Pi() << " deg";
357 
358  /*TVector3 direction1 = TVector3(sin(thetaLab1)*cos(phiLab1),
359  sin(thetaLab1)*sin(phiLab1),
360  cos(thetaLab1));*/
361 
362  TVector3 direction1 = TVector3(sin(thetaLab1) * cos(phiLab1), sin(thetaLab1) * sin(phiLab1), cos(thetaLab1));
363 
364  EulerTransformer->SetThetaInBeamSystem(Ang.at(1));
365  EulerTransformer->SetPhiInBeamSystem(phiBeam2);
366  EulerTransformer->DoTheEulerTransformationBeam2Lab(); // Euler transformation for particle 2
367  // EulerTransformer->PrintResults();
368 
369  thetaLab2 = EulerTransformer->GetThetaInLabSystem();
370  phiLab2 = EulerTransformer->GetPhiInLabSystem();
371 
372  TVector3 direction2 = TVector3(sin(thetaLab2) * cos(phiLab2), sin(thetaLab2) * sin(phiLab2), cos(thetaLab2));
373 
374  LOG(DEBUG) << " Recoiled angle Phi :" << phiBeam2 * 180.0 / TMath::Pi() << " deg";
375  LOG(DEBUG) << " Recoiled angle Theta (Euler) :" << thetaLab2 * 180.0 / TMath::Pi() << " deg";
376  LOG(DEBUG) << " Recoiled angle Phi (Euler) :" << phiLab2 * 180.0 / TMath::Pi() << " deg";
377 
378  LOG(DEBUG) << " Phi Diference :" << (phiBeam1 * 180.0 / TMath::Pi()) - (phiBeam2 * 180.0 / TMath::Pi())
379  << " deg";
380  LOG(DEBUG) << " Phi Diference (Euler) :" << (phiLab1 * 180.0 / TMath::Pi()) - (phiLab2 * 180.0 / TMath::Pi())
381  << " deg";
382 
383  LOG(DEBUG) << " Direction 1 Theta : " << direction1.Theta() * 180.0 / TMath::Pi();
384  LOG(DEBUG) << " Direction 1 Phi : " << direction1.Phi() * 180.0 / TMath::Pi();
385  LOG(DEBUG) << " Direction 2 Theta : " << direction2.Theta() * 180.0 / TMath::Pi();
386  LOG(DEBUG) << " Direction 2 Phi : " << direction2.Phi() * 180.0 / TMath::Pi();
387 
388  fPx.at(2) = p3_lab * direction1.X() / 1000.0; // To GeV for FairRoot
389  fPy.at(2) = p3_lab * direction1.Y() / 1000.0;
390  fPz.at(2) = p3_lab * direction1.Z() / 1000.0;
391 
392  fPx.at(3) = p4_lab * direction2.X() / 1000.0;
393  fPy.at(3) = p4_lab * direction2.Y() / 1000.0;
394  fPz.at(3) = p4_lab * direction2.Z() / 1000.0;
395 
396  } else {
397 
398  fPx.at(2) = 0.0; // To GeV for FairRoot
399  fPy.at(2) = 0.0;
400  fPz.at(2) = 0.0;
401 
402  fPx.at(3) = 0.0;
403  fPy.at(3) = 0.0;
404  fPz.at(3) = 0.0;
405  }
406 
407  /* if(!AtVertexPropagator::Instance()->GetValidKine()){
408 
409  fPx.at(2) = 0.0; // To GeV for FairRoot
410  fPy.at(2) = 0.0;
411  fPz.at(2) = 0.0;
412 
413  fPx.at(3) = 0.0;
414  fPy.at(3) = 0.0;
415  fPz.at(3) = 0.0;
416 
417  }else{
418 
419  fPx.at(2) = p3_lab*direction1.X()/1000.0; // To GeV for FairRoot
420  fPy.at(2) = p3_lab*direction1.Y()/1000.0;
421  fPz.at(2) = p3_lab*direction1.Z()/1000.0;
422 
423  fPx.at(3) = p4_lab*direction2.X()/1000.0;
424  fPy.at(3) = p4_lab*direction2.Y()/1000.0;
425  fPz.at(3) = p4_lab*direction2.Z()/1000.0;
426 
427  }*/
428 
429  // Debugging purposes
430 
431  /*TVector3 debug(fPx.at(2),fPy.at(2),fPz.at(2));
432  std::cout<<" p3_lab : "<<p3_lab<<std::endl;
433  std::cout<<" Debug momentum : "<<debug.Mag()<<std::endl;
434  std::cout<<" Debug Phi : "<<debug.Phi()*180.0/TMath::Pi()<<std::endl;
435  std::cout<<" Debug Theta : "<<debug.Theta()*180.0/TMath::Pi()<<std::endl;
436  TVector3 debug2(fPx.at(3),fPy.at(3),fPz.at(3));
437  std::cout<<" p4_lab : "<<p4_lab<<std::endl;
438  std::cout<<" Debug momentum 2 : "<<debug2.Mag()<<std::endl;
439  std::cout<<" Debug Phi 2 : "<<debug2.Phi()*180.0/TMath::Pi()<<std::endl;
440  std::cout<<" Debug Theta 2 : "<<debug2.Theta()*180.0/TMath::Pi()<<std::endl;*/
441 
442  // Particle transport begins here
443 
444  for (Int_t i = 0; i < fMult; i++) {
445 
446  TParticlePDG *thisPart = nullptr;
447 
448  if (fPType.at(i) == "Ion")
449  thisPart = TDatabasePDG::Instance()->GetParticle(fIon.at(i)->GetName());
450  else if (fPType.at(i) == "Proton")
451  thisPart = TDatabasePDG::Instance()->GetParticle(fParticle.at(i)->GetName());
452  else if (fPType.at(i) == "Neutron")
453  thisPart = TDatabasePDG::Instance()->GetParticle(fParticle.at(i)->GetName());
454 
455  if (!thisPart) {
456 
457  if (fPType.at(i) == "Ion")
458  std::cout << "-W- FairIonGenerator: Ion " << fIon.at(i)->GetName() << " not found in database!"
459  << std::endl;
460  else if (fPType.at(i) == "Proton")
461  std::cout << "-W- FairIonGenerator: Particle " << fParticle.at(i)->GetName() << " not found in database!"
462  << std::endl;
463  else if (fPType.at(i) == "Neutron")
464  std::cout << "-W- FairIonGenerator: Particle " << fParticle.at(i)->GetName() << " not found in database!"
465  << std::endl;
466 
467  return kFALSE;
468  }
469 
470  int pdgType = thisPart->PdgCode();
471 
472  // Propagate the vertex of the previous event
473  if (fIsFixedTargetPos) {
474 
475  // fVx = 0.0;
476  // fVy = 0.0;
477  // fVz = 0.0;
478 
479  } else {
480  fVx = AtVertexPropagator::Instance()->GetVx();
481  fVy = AtVertexPropagator::Instance()->GetVy();
482  fVz = AtVertexPropagator::Instance()->GetVz();
483  }
484 
485  // For decay generators
486  TVector3 ScatP(fPx.at(2), fPy.at(2), fPz.at(2));
487  AtVertexPropagator::Instance()->SetScatterP(ScatP);
488  AtVertexPropagator::Instance()->SetScatterEx(fExEnergy.at(2));
489 
490  Int_t trackIdCut = 0;
491 
492  if (kIsDecay)
493  trackIdCut = 2; // Remove beam and decaying particle
494  else
495  trackIdCut = 1; // Remove beam
496 
497  if (i > trackIdCut && (AtVertexPropagator::Instance()->GetDecayEvtCnt() || fIsFixedTargetPos) &&
498  pdgType != 1000500500 && fPType.at(i) == "Ion") {
499 
500  std::cout << cBLUE << "-I- FairIonGenerator: Generating ions of type " << fIon.at(i)->GetName()
501  << " (PDG code " << pdgType << ")" << std::endl;
502  std::cout << " Momentum (" << fPx.at(i) << ", " << fPy.at(i) << ", " << fPz.at(i)
503  << ") Gev from vertex (" << fVx << ", " << fVy << ", " << fVz << ") cm" << std::endl;
504  primGen->AddTrack(pdgType, fPx.at(i), fPy.at(i), fPz.at(i), fVx, fVy, fVz);
505 
506  } else if (i > 1 && (AtVertexPropagator::Instance()->GetDecayEvtCnt() || fIsFixedTargetPos) &&
507  pdgType == 2212 && fPType.at(i) == "Proton") {
508 
509  std::cout << "-I- FairIonGenerator: Generating ions of type " << fParticle.at(i)->GetName() << " (PDG code "
510  << pdgType << ")" << std::endl;
511  std::cout << " Momentum (" << fPx.at(i) << ", " << fPy.at(i) << ", " << fPz.at(i)
512  << ") Gev from vertex (" << fVx << ", " << fVy << ", " << fVz << ") cm" << std::endl;
513  primGen->AddTrack(pdgType, fPx.at(i), fPy.at(i), fPz.at(i), fVx, fVy, fVz);
514 
515  } else if (i > 1 && (AtVertexPropagator::Instance()->GetDecayEvtCnt() || fIsFixedTargetPos) &&
516  pdgType == 2112 && fPType.at(i) == "Neutron") {
517 
518  std::cout << "-I- FairIonGenerator: Generating ions of type " << fParticle.at(i)->GetName() << " (PDG code "
519  << pdgType << ")" << std::endl;
520  std::cout << " Momentum (" << fPx.at(i) << ", " << fPy.at(i) << ", " << fPz.at(i)
521  << ") Gev from vertex (" << fVx << ", " << fVy << ", " << fVz << ") cm" << cNORMAL << std::endl;
522  primGen->AddTrack(pdgType, fPx.at(i), fPy.at(i), fPz.at(i), fVx, fVy, fVz);
523  }
524  }
525 
526  } // if residual energy > 0
527 
528  if (kIsDecay == kFALSE) // Only increases the reaction counter if decay is not expected
529  AtVertexPropagator::Instance()->IncDecayEvtCnt();
530 
531  return kTRUE;
532 }
533 
534 ClassImp(AtTPC2Body)
cNORMAL
constexpr auto cNORMAL
Definition: AtTPC2Body.cxx:27
AtVertexPropagator::SetValidKine
void SetValidKine(Bool_t val)
Definition: AtVertexPropagator.cxx:325
AtVertexPropagator::GetVx
Double_t GetVx()
Definition: AtVertexPropagator.cxx:196
AtVertexPropagator::GetVz
Double_t GetVz()
Definition: AtVertexPropagator.cxx:204
AtVertexPropagator::GetPx
Double_t GetPx()
Definition: AtVertexPropagator.cxx:220
AtTPC2Body::AtTPC2Body
AtTPC2Body()
Definition: AtTPC2Body.cxx:33
AtVertexPropagator::GetPz
Double_t GetPz()
Definition: AtVertexPropagator.cxx:228
AtVertexPropagator::GetPy
Double_t GetPy()
Definition: AtVertexPropagator.cxx:224
ClassImp
ClassImp(AtFindVertex)
cGREEN
constexpr auto cGREEN
Definition: AtTPC2Body.cxx:28
AtVertexPropagator::SetScatterP
void SetScatterP(TVector3 val)
Definition: AtVertexPropagator.cxx:161
cRED
constexpr auto cRED
Definition: AtTPC2Body.cxx:25
cBLUE
constexpr auto cBLUE
Definition: AtTPC2Body.cxx:29
AtVertexPropagator::GetVy
Double_t GetVy()
Definition: AtVertexPropagator.cxx:200
AtVertexPropagator::GetEnergy
Double_t GetEnergy()
Definition: AtVertexPropagator.cxx:232
AtVertexPropagator.h
AtTPC2Body::ReadEvent
virtual Bool_t ReadEvent(FairPrimaryGenerator *primGen)
Definition: AtTPC2Body.cxx:148
AtVertexPropagator::SetTrackEnergy
void SetTrackEnergy(int trackID, double energy)
Definition: AtVertexPropagator.cxx:92
AtTPC2Body.h
AtTPC2Body
Definition: AtTPC2Body.h:25
cYELLOW
constexpr auto cYELLOW
Definition: AtTPC2Body.cxx:26
AtVertexPropagator::SetTrackAngle
void SetTrackAngle(int trackID, double angle)
Definition: AtVertexPropagator.cxx:96
AtTPC2Body::SetFixedTargetPosition
void SetFixedTargetPosition(double vx, double vy, double vz)
Definition: AtTPC2Body.cxx:129
AtVertexPropagator::SetScatterEx
void SetScatterEx(Double_t val)
Definition: AtVertexPropagator.cxx:165
AtVertexPropagator::Instance
static AtVertexPropagator * Instance()
Definition: AtVertexPropagator.cxx:18
AtTPC2Body::SetFixedBeamMomentum
void SetFixedBeamMomentum(double px, double py, double pz)
Definition: AtTPC2Body.cxx:138
AtVertexPropagator::IncDecayEvtCnt
void IncDecayEvtCnt()
Definition: AtVertexPropagator.cxx:321
AtEulerTransformation.h
Generated on Sat Nov 11 2023 15:30:37 for ATTPCROOT by   doxygen 1.8.18