ATTPCROOT  0.3.0-alpha
A ROOT-based framework for analyzing data from active target detectors
AtTPC_d2He.cxx
Go to the documentation of this file.
1 #include "AtTPC_d2He.h"
2 
3 #include "AtVertexPropagator.h"
4 
5 #include <FairIon.h>
6 #include <FairParticle.h>
7 #include <FairPrimaryGenerator.h>
8 #include <FairRunSim.h>
9 
10 #include <TDatabasePDG.h>
11 #include <TMath.h>
12 #include <TMathBase.h>
13 #include <TParticle.h>
14 #include <TParticlePDG.h>
15 #include <TRandom.h>
16 #include <TVector3.h>
17 
18 #include <algorithm>
19 #include <cmath>
20 #include <cstdio>
21 #include <iostream>
22 
23 constexpr auto cRED = "\033[1;31m";
24 constexpr auto cYELLOW = "\033[1;33m";
25 constexpr auto cNORMAL = "\033[0m";
26 constexpr auto cGREEN = "\033[1;32m";
27 constexpr auto cBLUE = "\033[1;34m";
28 
29 constexpr float amu = 931.494;
30 
31 Int_t AtTPC_d2He::fgNIon = 0;
32 
33 AtTPC_d2He::AtTPC_d2He() : fMult(0), fPx(0.), fPy(0.), fPz(0.), fVx(0.), fVy(0.), fVz(0.), fIon(0)
34 {
35  // cout << "-W- AtTPCIonGenerator: "
36  // << " Please do not use the default constructor! " << endl;
37 }
38 
39 inline Double_t sign(Double_t num)
40 {
41  if (num > 0)
42  return 1;
43  return (num == 0) ? 1 : -1;
44 }
45 
46 // ----- Default constructor ------------------------------------------
47 AtTPC_d2He::AtTPC_d2He(const char *name, std::vector<Int_t> *z, std::vector<Int_t> *a, std::vector<Int_t> *q,
48  Int_t mult, std::vector<Double_t> *px, std::vector<Double_t> *py, std::vector<Double_t> *pz,
49  std::vector<Double_t> *mass, std::vector<Double_t> *Ex, std::vector<Double_t> *cross1,
50  std::vector<Double_t> *cross2, std::vector<Double_t> *cross3, Int_t N_data)
51  : fPx(0.), fPy(0.), fPz(0.), fMult(mult), fVx(0.), fVy(0.), fVz(0.), fIon(0)
52 {
53 
54  fgNIon++;
55 
56  fIon.reserve(fMult);
57 
58  char buffer[30];
59  auto *kProton = new TParticle(); // NOLINT but probably a problem
60  kProton->SetPdgCode(2212);
61 
62  auto *kNeutron = new TParticle(); // NOLINT but probably a problem
63  kNeutron->SetPdgCode(2112);
64 
65  // Read the cross section table
66  for (Int_t i = 0; i < N_data; i++) {
67  inp1.push_back(cross1->at(i));
68  inp2.push_back(cross2->at(i));
69  inp3.push_back(cross3->at(i));
70  fCStot += inp3.at(i);
71  // std::cout<<"===================================================================="<<std::endl;
72  // std::cout<<inp1.at(i)<<" "<<inp2.at(i)<<" "<<inp3.at(i)<<std::endl;
73  // std::cout<<"===================================================================="<<std::endl;
74  }
75 
76  fN = N_data;
77 
78  for (Int_t i = 0; i < fMult; i++) {
79 
80  fPx.push_back(px->at(i));
81  fPy.push_back(py->at(i));
82  fPz.push_back(pz->at(i));
83  Masses.push_back(mass->at(i));
84  fExEnergy.push_back(Ex->at(i));
85  // fWm.push_back( mass->at(i));
86 
87  FairIon *IonBuff;
88  FairParticle *ParticleBuff;
89  sprintf(buffer, "Product_Ion%d", i);
90  if (a->at(i) != 1) {
91  IonBuff = new FairIon(buffer, z->at(i), a->at(i), q->at(i), 0.0, // NOLINT but probably a problem
92  mass->at(i) * amu / 1000.0);
93  ParticleBuff = new FairParticle("dummyPart", 1, 1, 1.0, 0, 0.0, 0.0); // NOLINT but probably a problem
94  fPType.emplace_back("Ion");
95  // std::cout<<" Adding : "<<buffer<<std::endl;
96 
97  } else if (a->at(i) == 1 && z->at(i) == 1) {
98 
99  IonBuff = new FairIon(buffer, z->at(i), a->at(i), q->at(i), 0.0, // NOLINT but probably a problem
100  mass->at(i) * amu / 1000.0);
101  ParticleBuff = new FairParticle(2212, kProton); // NOLINT but probably a problem
102  fPType.emplace_back("Proton");
103 
104  } else if (a->at(i) == 1 && z->at(i) == 0) {
105 
106  IonBuff = new FairIon(buffer, z->at(i), a->at(i), q->at(i), 0.0, // NOLINT but probably a problem
107  mass->at(i) * amu / 1000.0);
108  ParticleBuff = new FairParticle(2112, kNeutron); // NOLINT but probably a problem
109  fPType.emplace_back("Neutron");
110  }
111 
112  // std::cout<<" Z "<<z->at(i)<<" A "<<a->at(i)<<std::endl;
113  // std::cout<<buffer<<std::endl;
114  fIon.push_back(IonBuff);
115  fParticle.push_back(ParticleBuff);
116  }
117 
118  FairRunSim *run = FairRunSim::Instance();
119  if (!run) {
120  std::cout << "-E- FairIonGenerator: No FairRun instantised!" << std::endl;
121  Fatal("FairIonGenerator", "No FairRun instantised!");
122  }
123 
124  for (Int_t i = 0; i < fMult; i++) {
125 
126  if (fPType.at(i) == "Ion") {
127  // std::cout<<" In position "<<i<<" adding an : "<<fPType.at(i)<<std::endl;
128  run->AddNewIon(fIon.at(i)); // NOLINT
129  // std::cout<<" fIon name :"<<fIon.at(i)->GetName()<<std::endl;
130  // std::cout<<" fParticle name :"<<fParticle.at(i)->GetName()<<std::endl;
131 
132  } else if (fPType.at(i) == "Proton") {
133  // std::cout<<" In position "<<i<<" adding an : "<<fPType.at(i)<<std::endl;
134  // run->AddNewParticle(fParticle.at(i));
135  // std::cout<<" fIon name :"<<fIon.at(i)->GetName()<<std::endl;
136  // std::cout<<" fParticle name :"<<fParticle.at(i)->GetName()<<std::endl;
137  // std::cout<<fParticle.at(i)->GetName()<<std::endl;
138 
139  } else if (fPType.at(i) == "Neutron") {
140 
141  // std::cout<<" In position "<<i<<" adding an : "<<fPType.at(i)<<std::endl;
142  // run->AddNewParticle(fParticle.at(i));
143  // std::cout<<" fIon name :"<<fIon.at(i)->GetName()<<std::endl;
144  // std::cout<<" fParticle name :"<<fParticle.at(i)->GetName()<<std::endl;
145  // std::cout<<fParticle.at(i)->GetName()<<std::endl;
146  }
147  }
148 }
149 
150 // Rotation of a 3D vector around an arbitrary axis
151 // Rodriges Formula
152 std::vector<Double_t>
153 AtTPC_d2He::TRANSF(std::vector<Double_t> *from, std::vector<Double_t> *to, std::vector<Double_t> *vin)
154 {
155 
156  static std::vector<double> vout(3);
157  double n[3];
158  double normn, normf, normt;
159  double alpha, a, b;
160 
161  normf = sqrt(pow(from->at(0), 2) + pow(from->at(1), 2) + pow(from->at(2), 2));
162  normt = sqrt(pow(to->at(0), 2) + pow(to->at(1), 2) + pow(to->at(2), 2));
163 
164  alpha =
165  acos(((from->at(0)) * (to->at(0)) + (from->at(1)) * (to->at(1)) + (from->at(2)) * (to->at(2))) / (normf * normt));
166  a = sin(alpha);
167  b = 1.0 - cos(alpha);
168 
169  if (fabs(alpha) < 0.000001) {
170  vout.at(0) = vin->at(0);
171  vout.at(1) = vin->at(1);
172  vout.at(2) = vin->at(2);
173 
174  } else {
175  n[0] = ((from->at(1)) * (to->at(2)) - (from->at(2)) * (to->at(1)));
176  n[1] = ((from->at(2)) * (to->at(0)) - (from->at(0)) * (to->at(2)));
177  n[2] = ((from->at(0)) * (to->at(1)) - (from->at(1)) * (to->at(0)));
178 
179  // std::cout<<from->at(0)<<" "<<from->at(1)<<" "<<from->at(2)<<std::endl;
180  // std::cout<<to->at(0)<<" "<<to->at(1)<<" "<<to->at(2)<<std::endl;
181  // std::cout<<vin->at(0)<<" "<<vin->at(1)<<" "<<vin->at(2)<<std::endl;
182  normn = sqrt(pow(n[0], 2) + pow(n[1], 2) + pow(n[2], 2));
183  n[0] = n[0] / normn;
184  n[1] = n[1] / normn;
185  n[2] = n[2] / normn;
186 
187  vout.at(0) = (1 - b * (pow(n[2], 2) + pow(n[1], 2))) * (vin->at(0)) +
188  (-a * n[2] + b * n[0] * n[1]) * (vin->at(1)) + (a * n[1] + b * n[0] * n[2]) * (vin->at(2));
189  vout.at(1) = (a * n[2] + b * n[0] * n[1]) * (vin->at(0)) +
190  (1 - b * (pow(n[2], 2) + pow(n[0], 2))) * (vin->at(1)) +
191  (-a * n[0] + b * n[1] * n[2]) * (vin->at(2));
192  vout.at(2) = (-a * n[1] + b * n[0] * n[2]) * (vin->at(0)) + (a * n[0] + b * n[1] * n[2]) * (vin->at(1)) +
193  (1 - b * (pow(n[1], 2) + pow(n[0], 2))) * (vin->at(2));
194  }
195 
196  return vout;
197 }
198 
199 Double_t AtTPC_d2He::omega(Double_t x, Double_t y, Double_t z)
200 {
201  return sqrt(x * x + y * y + z * z - 2 * x * y - 2 * y * z - 2 * x * z);
202 }
203 
204 // ----- Public method ReadEvent --------------------------------------
205 Bool_t AtTPC_d2He::ReadEvent(FairPrimaryGenerator *primGen)
206 {
207 
208  std::vector<Double_t> Ang; // Lab Angle of the products
209  std::vector<Double_t> Ene; // Lab Energy of the products
210  Ang.resize(4);
211  Ene.resize(4);
212 
213  fIsDecay = kFALSE;
214 
215  fBeamEnergy = AtVertexPropagator::Instance()->GetEnergy();
216  std::cout << " -I- AtTPC_d2He Residual energy : " << AtVertexPropagator::Instance()->GetEnergy() << std::endl;
217 
218  fPxBeam = AtVertexPropagator::Instance()->GetPx();
219  fPyBeam = AtVertexPropagator::Instance()->GetPy();
220  fPzBeam = AtVertexPropagator::Instance()->GetPz();
221 
222  // fPxBeam = fPx.at(0) ;
223  // fPyBeam = fPy.at(0) ;
224  // fPzBeam = fPz.at(0) ;
225 
226  if (fBeamEnergy == 0) {
227  std::cout << "-I- AtTP_d2He : No solution!" << std::endl;
228  AtVertexPropagator::Instance()->SetValidKine(kFALSE);
229  }
230 
231  if (!AtVertexPropagator::Instance()->GetValidKine()) {
232 
233  fPx.at(2) = 0.; // To GeV for FairRoot
234  fPy.at(2) = 0.;
235  fPz.at(2) = 0.;
236 
237  fPx.at(3) = 0.;
238  fPy.at(3) = 0.;
239  fPz.at(3) = 0.;
240 
241  fPx.at(4) = 0.;
242  fPy.at(4) = 0.;
243  fPz.at(4) = 0.;
244 
245  fPx.at(5) = 0.;
246  fPy.at(5) = 0.;
247  fPz.at(5) = 0.;
248 
249  Ene.at(0) = 0.0;
250  Ang.at(0) = 0.0;
251  Ene.at(1) = 0.0;
252  Ang.at(1) = 0.0;
253  Ene.at(2) = 0.0;
254  Ang.at(2) = 0.0;
255  Ene.at(3) = 0.0;
256  Ang.at(3) = 0.0;
257  }
258 
259  else {
260  // MC to distribute the events with the cross section
261  // fN==1 used for the efficiency map (1 simu per theta-epp bin)
262  if (fN == 1) { // Depp=0.25 and Dtheta_cm=0.5 in ACCBA file
263  if (inp2.at(0) == 0)
264  theta_cm = (inp2.at(0) + 0.25 * gRandom->Uniform()) *
265  TMath::DegToRad(); // carefull in the analysis, these bins have 2 times more stat
266  else
267  theta_cm = (inp2.at(0) - 0.25 + 0.5 * gRandom->Uniform()) * TMath::DegToRad();
268  phi_cm = 2 * TMath::Pi() * (gRandom->Uniform());
269  epsilon = inp1.at(0) - 0.125 + 0.25 * gRandom->Uniform();
270  } else {
271  do {
272  ran_theta = fN * gRandom->Uniform();
273  ranX = fCStot * gRandom->Uniform();
274  } while (ranX > inp3.at(ran_theta));
275 
276  // Depp=0.25 and Dtheta_cm=0.5 in ACCBA file, here Depp=0.5 and Dtheta_cm=1 introduce smearing
277  theta_cm = TMath::Abs(inp2.at(ran_theta) - 0.5 + gRandom->Uniform()) * TMath::DegToRad();
278  phi_cm = 2 * TMath::Pi() * (gRandom->Uniform());
279  epsilon = TMath::Abs(inp1.at(ran_theta) - 0.25 + 0.5 * gRandom->Uniform());
280  }
281  // std::cout<<"===================================================================="<<std::endl;
282  // std::cout<<theta_cm*TMath::RadToDeg()<<" "<<phi_cm<<" "<<epsilon<<std::endl;
283  // std::cout<<fCStot<<std::endl;
284  // //std::cout<<fPxBeam<<" "<<fPyBeam<<" "<<fPzBeam<<std::endl;
285  // std::cout<<"===================================================================="<<std::endl;
286 
287  // dirty way to include more than one excited state
288  /*test_var = gRandom->Uniform();
289  if(test_var>= 0 && test_var<0.25) Ex_ejectile = 0.0;
290  if(test_var>= 0.25 && test_var<0.50) Ex_ejectile = 5.0;
291  if(test_var>= 0.50 && test_var<0.75) Ex_ejectile = 10.0;
292  if(test_var>= 0.75 && test_var<1.0) Ex_ejectile = 20.0;
293  */
294  Ex_ejectile = fExEnergy.at(2);
295 
296  m1 = Masses.at(0) * amu + fExEnergy.at(0);
297  m2 = Masses.at(1) * amu + fExEnergy.at(1);
298  m3 = Masses.at(2) * amu + Ex_ejectile; // ejectile
299  m4 = Masses.at(3) * amu + epsilon; // 2he
300  m7 = Masses.at(4) * amu;
301  m8 = Masses.at(5) * amu;
302  // K1 = sqrt(pow(fPx.at(0),2) + pow( fPy.at(0),2) + pow( fPz.at(0),2) + pow(m1,2)) - m1;
303  K1 = sqrt(pow(fPxBeam * 1000.0, 2) + pow(fPyBeam * 1000.0, 2) + pow(fPzBeam * 1000.0, 2) + pow(m1, 2)) - m1;
304 
305  p1L[0] = fPxBeam * 1000.0;
306  p1L[1] = fPyBeam * 1000.0;
307  p1L[2] = fPzBeam * 1000.0;
308  // p1L[0] = fPx.at(0);
309  // p1L[1] = fPy.at(0);
310  // p1L[2] = fPz.at(0);
311 
312  E1L = K1 + m1;
313 
314  // get COM parameters
315  beta_cm = p1L[2] / (E1L + m2);
316  gamma_cm = 1.0 / sqrt(1.0 - pow(beta_cm, 2));
317  S = 2. * E1L * m2 + pow(m1, 2) + pow(m2, 2);
318  Pcm = 0.5 * (AtTPC_d2He::omega(S, pow(m3, 2), pow(m4, 2))) / sqrt(S);
319 
320  // std::cout<<beta_cm<<" "<<gamma_cm<<" "<<p1L[2]<<" "<<K1<<" "<<m1<<" "<<(AtTPC_d2He::omega(S, pow(m3,2),
321  // pow(m4,2)))<<" "<<sqrt(S)<<std::endl;
322 
323  // generate cm angles and momenta (for now isotropic) and corresponding momenta
324  p4C[2] = Pcm * cos(TMath::Pi() - theta_cm); // Pi -thethacm because inv kinematics
325  p4C[0] = Pcm * sin(TMath::Pi() - theta_cm) * cos(phi_cm);
326  p4C[1] = Pcm * sin(TMath::Pi() - theta_cm) * sin(phi_cm);
327  E4C = sqrt(pow(Pcm, 2) + pow(m4, 2));
328 
329  for (int i = 0; i < 3; i++) {
330  p3C[i] = -1. * p4C[i];
331  }
332  E3C = sqrt(pow(Pcm, 2) + pow(m3, 2));
333 
334  // transformation to the lab frame
335  p3L[0] = p3C[0];
336  p3L[1] = p3C[1];
337  p3L[2] = gamma_cm * (p3C[2] + beta_cm * E3C);
338  E3L = sqrt(pow(p3L[0], 2) + pow(p3L[1], 2) + pow(p3L[2], 2) + pow(m3, 2));
339 
340  p4L[0] = p4C[0];
341  p4L[1] = p4C[1];
342  p4L[2] = gamma_cm * (p4C[2] + beta_cm * E4C);
343  E4L = sqrt(pow(p4L[0], 2) + pow(p4L[1], 2) + pow(p4L[2], 2) + pow(m4, 2));
344 
345  // std::cout<<E3L-m3<<" "<<E4L-m4<<" "<<p3L[2]<<" "<<p4L[2]<<std::endl;
346 
347  // rotate back to the beam axis
348  fvto.clear();
349  fvto.resize(3);
350  fvin.clear();
351  fvin.resize(3);
352  fvout.clear();
353  fvout.resize(3);
354  fvfrom.clear();
355  fvfrom.resize(3);
356  fvfrom.at(0) = 0;
357  fvfrom.at(1) = 0;
358  fvfrom.at(2) = 1;
359  fvto.at(0) = p1L[0];
360  fvto.at(1) = p1L[1];
361  fvto.at(2) = p1L[2];
362  fvin.at(0) = p3L[0];
363  fvin.at(1) = p3L[1];
364  fvin.at(2) = p3L[2];
365  fvout = AtTPC_d2He::TRANSF(&fvfrom, &fvto, &fvin);
366  p3L[0] = fvout.at(0);
367  p3L[1] = fvout.at(1);
368  p3L[2] = fvout.at(2);
369  E3L = sqrt(pow(p3L[0], 2) + pow(p3L[1], 2) + pow(p3L[2], 2) + pow(m3, 2));
370 
371  fvin.at(0) = p4L[0];
372  fvin.at(1) = p4L[1];
373  fvin.at(2) = p4L[2];
374  fvout = AtTPC_d2He::TRANSF(&fvfrom, &fvto, &fvin);
375  p4L[0] = fvout.at(0);
376  p4L[1] = fvout.at(1);
377  p4L[2] = fvout.at(2);
378  E4L = sqrt(pow(p4L[0], 2) + pow(p4L[1], 2) + pow(p4L[2], 2) + pow(m4, 2));
379 
380  // Particle 4 is the 2He which is unbound and will thus decay into 2 protons
381  normP4L = sqrt(pow(p4L[0], 2) + pow(p4L[1], 2) + pow(p4L[2], 2));
382  beta4 = normP4L / E4L;
383  gamma4 = 1.0 / sqrt(1.0 - pow(beta4, 2));
384  S_78 = pow(m4, 2);
385  Pc78 = 0.5 * AtTPC_d2He::omega(S_78, pow(m7, 2), pow(m8, 2)) / sqrt(S_78);
386 
387  //-----------generate isotropically theta and phi of particles 7 and 8
388  ran1 = (gRandom->Uniform());
389  ran2 = (gRandom->Uniform());
390  theta78 = acos(2 * ran1 - 1.);
391  phi78 = 2 * TMath::Pi() * ran2;
392 
393  // generate the protons in the 2He rest frame
394  p7rest[2] = Pc78 * cos(theta78);
395  p7rest[0] = Pc78 * sin(theta78) * cos(phi78);
396  p7rest[1] = Pc78 * sin(theta78) * sin(phi78);
397  p8rest[2] = -1 * p7rest[2];
398  p8rest[0] = -1 * p7rest[0];
399  p8rest[1] = -1 * p7rest[1];
400  E7rest = sqrt(pow(Pc78, 2) + pow(m7, 2));
401  E8rest = sqrt(pow(Pc78, 2) + pow(m8, 2));
402 
403  // boost to 2He frame
404  p7L[0] = p7rest[0];
405  p7L[1] = p7rest[1];
406  p7L[2] = gamma4 * (p7rest[2] + beta4 * E7rest);
407  E7L = sqrt(pow(m7, 2) + pow(p7L[0], 2) + pow(p7L[1], 2) + pow(p7L[2], 2));
408 
409  p8L[0] = p8rest[0];
410  p8L[1] = p8rest[1];
411  p8L[2] = gamma4 * (p8rest[2] + beta4 * E8rest);
412  E8L = sqrt(pow(m8, 2) + pow(p8L[0], 2) + pow(p8L[1], 2) + pow(p8L[2], 2));
413 
414  // rotate to the 2He direction
415  fvto.at(0) = p4L[0];
416  fvto.at(1) = p4L[1];
417  fvto.at(2) = p4L[2];
418  fvin.at(0) = p7L[0];
419  fvin.at(1) = p7L[1];
420  fvin.at(2) = p7L[2];
421  fvout = AtTPC_d2He::TRANSF(&fvfrom, &fvto, &fvin);
422  p7L[0] = fvout.at(0);
423  p7L[1] = fvout.at(1);
424  p7L[2] = fvout.at(2);
425 
426  fvin.at(0) = p8L[0];
427  fvin.at(1) = p8L[1];
428  fvin.at(2) = p8L[2];
429  fvout = AtTPC_d2He::TRANSF(&fvfrom, &fvto, &fvin);
430  p8L[0] = fvout.at(0);
431  p8L[1] = fvout.at(1);
432  p8L[2] = fvout.at(2);
433 
434  E7L = sqrt(pow(m7, 2) + pow(p7L[0], 2) + pow(p7L[1], 2) + pow(p7L[2], 2));
435  E8L = sqrt(pow(m8, 2) + pow(p8L[0], 2) + pow(p8L[1], 2) + pow(p8L[2], 2));
436 
437  fPx.at(2) = p3L[0] / 1000.0; // To GeV for FairRoot
438  fPy.at(2) = p3L[1] / 1000.0;
439  fPz.at(2) = p3L[2] / 1000.0;
440 
441  fPx.at(3) = p4L[0] / 1000.0; // To GeV for FairRoot
442  fPy.at(3) = p4L[1] / 1000.0;
443  fPz.at(3) = p4L[2] / 1000.0;
444 
445  fPx.at(4) = p7L[0] / 1000.0; // To GeV for FairRoot
446  fPy.at(4) = p7L[1] / 1000.0;
447  fPz.at(4) = p7L[2] / 1000.0;
448 
449  fPx.at(5) = p8L[0] / 1000.0; // To GeV for FairRoot
450  fPy.at(5) = p8L[1] / 1000.0;
451  fPz.at(5) = p8L[2] / 1000.0;
452 
453  Ene.at(0) = E3L - m3; // beam like particle
454  Ang.at(0) = acos((p1L[0] * p3L[0] + p1L[1] * p3L[1] + p1L[2] * p3L[2]) /
455  (sqrt(pow(p1L[0], 2) + pow(p1L[1], 2) + pow(p1L[2], 2)) *
456  sqrt(pow(p3L[0], 2) + pow(p3L[1], 2) + pow(p3L[2], 2)))) *
457  TMath::RadToDeg();
458  Ene.at(1) = E4L - m4; // 2He
459  Ang.at(1) = acos((p1L[0] * p4L[0] + p1L[1] * p4L[1] + p1L[2] * p4L[2]) /
460  (sqrt(pow(p1L[0], 2) + pow(p1L[1], 2) + pow(p1L[2], 2)) *
461  sqrt(pow(p4L[0], 2) + pow(p4L[1], 2) + pow(p4L[2], 2)))) *
462  TMath::RadToDeg();
463  Ene.at(2) = E7L - m7; // proton 1
464  Ang.at(2) = acos((p1L[0] * p7L[0] + p1L[1] * p7L[1] + p1L[2] * p7L[2]) /
465  (sqrt(pow(p1L[0], 2) + pow(p1L[1], 2) + pow(p1L[2], 2)) *
466  sqrt(pow(p7L[0], 2) + pow(p7L[1], 2) + pow(p7L[2], 2)))) *
467  TMath::RadToDeg();
468  Ene.at(3) = E8L - m8;
469  Ang.at(3) = acos((p1L[0] * p8L[0] + p1L[1] * p8L[1] + p1L[2] * p8L[2]) /
470  (sqrt(pow(p1L[0], 2) + pow(p1L[1], 2) + pow(p1L[2], 2)) *
471  sqrt(pow(p8L[0], 2) + pow(p8L[1], 2) + pow(p8L[2], 2)))) *
472  TMath::RadToDeg();
473 
474  if (std::isnan(Ene.at(0)) || std::isnan(Ene.at(1)) || std::isnan(Ene.at(2)) || std::isnan(Ene.at(3))) {
475 
476  fPx.at(2) = 0.; // To GeV for FairRoot
477  fPy.at(2) = 0.;
478  fPz.at(2) = 0.;
479 
480  fPx.at(3) = 0.;
481  fPy.at(3) = 0.;
482  fPz.at(3) = 0.;
483 
484  fPx.at(4) = 0.;
485  fPy.at(4) = 0.;
486  fPz.at(4) = 0.;
487 
488  fPx.at(5) = 0.;
489  fPy.at(5) = 0.;
490  fPz.at(5) = 0.;
491  }
492 
493  } // if solution is valid
494 
495  /*
496  Double_t phi7 = atan2(p7L[1], p7L[0]) * TMath::RadToDeg();
497  if (phi7 < 0)
498  phi7 = (phi7 + 360.0);
499 
500  Double_t phi8 = atan2(p8L[1], p8L[0]) * TMath::RadToDeg();
501  if (phi8 < 0)
502  phi8 = (phi8 + 360.0);
503  */
504  // std::cout << " -I- ===== AtTPC_d2He - Kinematics ====== "<<Ex_ejectile<<std::endl;
505  // std::cout << " Scattered energy:" << Ene.at(0) << " MeV" << std::endl;
506  // std::cout << " Scattered angle:" << Ang.at(0) << " deg" << std::endl;
507  // std::cout << " proton1 energy:" << Ene.at(2) << " MeV" << std::endl;
508  // std::cout << " proton1 angle:" << Ang.at(2) << " deg" << std::endl;
509  // std::cout << " proton1 angle phi:" << phi78 << " deg" << std::endl;
510  // std::cout << " proton2 energy:" << Ene.at(3) << " MeV" << std::endl;
511  // std::cout << " proton2 angle:" << Ang.at(3) << " deg" << std::endl;
512  // // std::cout << " proton2 angle phi:" << phi78 << " deg" << std::endl;
513  // std::cout << " 2He kinetic energy:" << Ene.at(1) << " MeV" << std::endl;
514  // std::cout << " 2He lab angle:" << Ang.at(1) << " deg" << std::endl;
515 
516  AtVertexPropagator::Instance()->SetTrackEnergy(0, Ene.at(0));
517  AtVertexPropagator::Instance()->SetTrackAngle(0, Ang.at(0));
518  AtVertexPropagator::Instance()->SetTrackEnergy(1, Ene.at(2));
519  AtVertexPropagator::Instance()->SetTrackAngle(1, Ang.at(2));
520  AtVertexPropagator::Instance()->SetTrackEnergy(2, Ene.at(3));
521  AtVertexPropagator::Instance()->SetTrackAngle(2, Ang.at(3));
522 
523  TVector3 ScatP(fPx.at(2), fPy.at(2), fPz.at(2));
524  AtVertexPropagator::Instance()->SetScatterP(ScatP);
525  AtVertexPropagator::Instance()->SetScatterEx(Ex_ejectile);
526 
527  /*
528  do{
529  random_z = 100.0*(gRandom->Uniform()); //cm
530 
531  random_r = 1.0*(gRandom->Gaus(0,1)); //cm
532  random_phi = 2.0*TMath::Pi()*(gRandom->Uniform()); //rad
533 
534  }while( fabs(random_r) > 4.7 ); //cut at 2 sigma
535  */
536  TVector3 d2HeVtx = AtVertexPropagator::Instance()->Getd2HeVtx();
537  fVx = d2HeVtx.X();
538  fVy = d2HeVtx.Y();
539  fVz = d2HeVtx.Z();
540  std::cout << cYELLOW << "vertex in AtTPC_d2He " << fVx << " " << fVy << " " << fVz << cNORMAL << std::endl;
541 
542  for (Int_t i = 0; i < fMult; i++) {
543  TParticlePDG *thisPart;
544 
545  if (fPType.at(i) == "Ion")
546  thisPart = TDatabasePDG::Instance()->GetParticle(fIon.at(i)->GetName());
547  else if (fPType.at(i) == "Proton")
548  thisPart = TDatabasePDG::Instance()->GetParticle(fParticle.at(i)->GetName());
549  else if (fPType.at(i) == "Neutron")
550  thisPart = TDatabasePDG::Instance()->GetParticle(fParticle.at(i)->GetName());
551 
552  if (!thisPart) {
553  if (fPType.at(i) == "Ion")
554  std::cout << "-W- FairIonGenerator: Ion " << fIon.at(i)->GetName() << " not found in database!"
555  << std::endl;
556  else if (fPType.at(i) == "Proton")
557  std::cout << "-W- FairIonGenerator: Particle " << fParticle.at(i)->GetName() << " not found in database!"
558  << std::endl;
559  else if (fPType.at(i) == "Neutron")
560  std::cout << "-W- FairIonGenerator: Particle " << fParticle.at(i)->GetName() << " not found in database!"
561  << std::endl;
562  return kFALSE;
563  }
564 
565  int pdgType = thisPart->PdgCode();
566 
567  // Propagate the vertex of the previous event
568 
569  // fVx = AtVertexPropagator::Instance()->GetVx();
570  // fVy = AtVertexPropagator::Instance()->GetVy();
571  // fVz = AtVertexPropagator::Instance()->GetVz();
572 
573  /*
574  fVx = random_r*cos(random_phi);
575  fVy = random_r*sin(random_phi);
576  fVz = random_z;
577  */
578  // cout<<AtVertexPropagator::Instance()->GetVx(); <<" "<<AtVertexPropagator::Instance()->GetVy();<<"
579  // "<<AtVertexPropagator::Instance()->GetVz();<<" "<<endl;
580 
581  // TVector3 d2HeVtx(fVx,fVy,fVz);
582  // AtVertexPropagator::Instance()->Setd2HeVtx(d2HeVtx);
583 
584  if (i > 1 && i != 3 && AtVertexPropagator::Instance()->GetDecayEvtCnt() && pdgType != 1000500500 &&
585  fPType.at(i) == "Ion") {
586  // TODO: Dirty way to propagate only the products (0 and 1 are beam and target respectively)
587  // i=3 is excluded because corresponds to 2He
588  /* std::cout << "-I- FairIonGenerator: Generating ions of type "
589  << fIon.at(i)->GetName() << " (PDG code " << pdgType << ")" << std::endl;
590  std::cout << " Momentum (" << fPx.at(i) << ", " << fPy.at(i) << ", " << fPz.at(i)
591  << ") Gev from vertex (" << fVx << ", " << fVy
592  << ", " << fVz << ") cm" << std::endl;*/
593 
594  primGen->AddTrack(pdgType, fPx.at(i), fPy.at(i), fPz.at(i), fVx, fVy, fVz);
595 
596  } else if (i > 1 && i != 3 && AtVertexPropagator::Instance()->GetDecayEvtCnt() && pdgType == 2212 &&
597  fPType.at(i) == "Proton") {
598 
599  /* std::cout << "-I- FairIonGenerator: Generating ions of type "
600  << fParticle.at(i)->GetName() << " (PDG code " << pdgType << ")" << std::endl;
601  std::cout << " Momentum (" << fPx.at(i) << ", " << fPy.at(i) << ", " << fPz.at(i)
602  << ") Gev from vertex (" << fVx << ", " << fVy
603  << ", " << fVz << ") cm" << std::endl;*/
604 
605  // primGen->AddTrack(pdgType, fPx.at(i), fPy.at(i), fPz.at(i), fVx, fVy, fVz);
606  primGen->AddTrack(pdgType, fPx.at(i), fPy.at(i), fPz.at(i), fVx, fVy, fVz);
607 
608  } else if (i > 1 && i != 3 && AtVertexPropagator::Instance()->GetDecayEvtCnt() && pdgType == 2112 &&
609  fPType.at(i) == "Neutron") {
610 
611  /* std::cout << "-I- FairIonGenerator: Generating ions of type "
612  << fParticle.at(i)->GetName() << " (PDG code " << pdgType << ")" << std::endl;
613  std::cout << " Momentum (" << fPx.at(i) << ", " << fPy.at(i) << ", " << fPz.at(i)
614  << ") Gev from vertex (" << fVx << ", " << fVy
615  << ", " << fVz << ") cm" << std::endl;
616  */
617  primGen->AddTrack(pdgType, fPx.at(i), fPy.at(i), fPz.at(i), fVx, fVy, fVz);
618  }
619  }
620 
621  AtVertexPropagator::Instance()
622  ->IncDecayEvtCnt(); // TODO: Okay someone should put a more suitable name but we are on a hurry...
623  AtVertexPropagator::Instance()->Getd2HeEvt();
624 
625  return kTRUE;
626 }
627 
628 ClassImp(AtTPC_d2He)
AtTPC_d2He
Definition: AtTPC_d2He.h:24
AtVertexPropagator::SetValidKine
void SetValidKine(Bool_t val)
Definition: AtVertexPropagator.cxx:325
amu
constexpr float amu
Definition: AtTPC_d2He.cxx:29
AtVertexPropagator::GetPx
Double_t GetPx()
Definition: AtVertexPropagator.cxx:220
cNORMAL
constexpr auto cNORMAL
Definition: AtTPC_d2He.cxx:25
cRED
constexpr auto cRED
Definition: AtTPC_d2He.cxx:23
AtVertexPropagator::GetPz
Double_t GetPz()
Definition: AtVertexPropagator.cxx:228
AtTPC_d2He::TRANSF
virtual std::vector< Double_t > TRANSF(std::vector< Double_t > *from, std::vector< Double_t > *to, std::vector< Double_t > *vin)
Definition: AtTPC_d2He.cxx:153
cYELLOW
constexpr auto cYELLOW
Definition: AtTPC_d2He.cxx:24
AtVertexPropagator::GetPy
Double_t GetPy()
Definition: AtVertexPropagator.cxx:224
ClassImp
ClassImp(AtFindVertex)
AtVertexPropagator::SetScatterP
void SetScatterP(TVector3 val)
Definition: AtVertexPropagator.cxx:161
AtVertexPropagator::Getd2HeEvt
Bool_t Getd2HeEvt()
Definition: AtVertexPropagator.cxx:284
AtTPC_d2He::ReadEvent
virtual Bool_t ReadEvent(FairPrimaryGenerator *primGen)
Definition: AtTPC_d2He.cxx:205
AtTPC_d2He::omega
virtual Double_t omega(Double_t x, Double_t y, Double_t z)
Definition: AtTPC_d2He.cxx:199
cBLUE
constexpr auto cBLUE
Definition: AtTPC_d2He.cxx:27
y
const double * y
Definition: lmcurve.cxx:20
AtVertexPropagator::Getd2HeVtx
TVector3 Getd2HeVtx()
Definition: AtVertexPropagator.cxx:308
AtVertexPropagator::GetEnergy
Double_t GetEnergy()
Definition: AtVertexPropagator.cxx:232
AtVertexPropagator.h
AtVertexPropagator::SetTrackEnergy
void SetTrackEnergy(int trackID, double energy)
Definition: AtVertexPropagator.cxx:92
AtTPC_d2He.h
AtTPC_d2He::AtTPC_d2He
AtTPC_d2He()
Definition: AtTPC_d2He.cxx:33
AtVertexPropagator::SetTrackAngle
void SetTrackAngle(int trackID, double angle)
Definition: AtVertexPropagator.cxx:96
sign
Double_t sign(Double_t num)
Definition: AtTPC_d2He.cxx:39
AtVertexPropagator::SetScatterEx
void SetScatterEx(Double_t val)
Definition: AtVertexPropagator.cxx:165
AtVertexPropagator::Instance
static AtVertexPropagator * Instance()
Definition: AtVertexPropagator.cxx:18
cGREEN
constexpr auto cGREEN
Definition: AtTPC_d2He.cxx:26
AtVertexPropagator::IncDecayEvtCnt
void IncDecayEvtCnt()
Definition: AtVertexPropagator.cxx:321
Generated on Sat Nov 11 2023 15:30:38 for ATTPCROOT by   doxygen 1.8.18