ATTPCROOT  0.3.0-alpha
A ROOT-based framework for analyzing data from active target detectors
AtTPC_Background.cxx
Go to the documentation of this file.
1 #include "AtTPC_Background.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 <TMath.h>
11 #include <TParticle.h>
12 #include <TRandom.h>
13 
14 #include <algorithm>
15 #include <cmath>
16 #include <cstdio>
17 #include <iostream>
18 #include <utility> // for move
19 
20 constexpr float amu = 931.494;
21 
22 Int_t AtTPC_Background::fgNIon = 0;
23 
24 AtTPC_Background::AtTPC_Background() : fMult(0), fPx(0.), fPy(0.), fPz(0.), fVx(0.), fVy(0.), fVz(0.), fIon(0)
25 {
26  // cout << "-W- AtTPCIonGenerator: "
27  // << " Please do not use the default constructor! " << endl;
28 }
29 
30 // ----- Default constructor ------------------------------------------
31 AtTPC_Background::AtTPC_Background(const char *name, std::vector<Int_t> *z, std::vector<Int_t> *a,
32  std::vector<Int_t> *q, Int_t mult, std::vector<Double_t> *px,
33  std::vector<Double_t> *py, std::vector<Double_t> *pz, std::vector<Double_t> *mass,
34  std::vector<Double_t> *Ex)
35  : fPx(0.), fPy(0.), fPz(0.), fMult(mult), fVx(0.), fVy(0.), fVz(0.), fIon(0)
36 {
37 
38  fgNIon++;
39 
40  fIon.reserve(fMult);
41 
42  char buffer[30];
43 
44  for (Int_t i = 0; i < fMult; i++) {
45 
46  fPx.push_back(px->at(i));
47  fPy.push_back(py->at(i));
48  fPz.push_back(pz->at(i));
49  Masses.push_back(mass->at(i));
50  fExEnergy.push_back(Ex->at(i));
51  // fWm.push_back( mass->at(i));
52 
53  std::unique_ptr<FairIon> IonBuff = nullptr;
54  std::unique_ptr<FairParticle> ParticleBuff = nullptr;
55  sprintf(buffer, "Product_Ion%d", i);
56  if (a->at(i) != 1) {
57  IonBuff = std::make_unique<FairIon>(buffer, z->at(i), a->at(i), q->at(i), 0.0, mass->at(i) * amu / 1000.0);
58  ParticleBuff = std::make_unique<FairParticle>("dummyPart", 1, 1, 1.0, 0, 0.0, 0.0);
59  fPType.emplace_back("Ion");
60  std::cout << " Adding : " << buffer << std::endl;
61 
62  } else if (a->at(i) == 1 && z->at(i) == 1) {
63 
64  IonBuff = std::make_unique<FairIon>(buffer, z->at(i), a->at(i), q->at(i), 0.0, mass->at(i) * amu / 1000.0);
65  auto *kProton = new TParticle(); // NOLINT
66  kProton->SetPdgCode(2212);
67  ParticleBuff = std::make_unique<FairParticle>(2212, kProton);
68  fPType.emplace_back("Proton");
69  }
70 
71  std::cout << " Z " << z->at(i) << " A " << a->at(i) << std::endl;
72  // std::cout<<buffer<<std::endl;
73  fIon.push_back(std::move(IonBuff));
74  fParticle.push_back(std::move(ParticleBuff));
75  }
76 
77  FairRunSim *run = FairRunSim::Instance();
78  if (!run) {
79  std::cout << "-E- FairIonGenerator: No FairRun instantised!" << std::endl;
80  Fatal("FairIonGenerator", "No FairRun instantised!");
81  }
82 
83  for (Int_t i = 0; i < fMult; i++) {
84 
85  if (fPType.at(i) == "Ion") {
86  std::cout << " In position " << i << " adding an : " << fPType.at(i) << std::endl;
87  run->AddNewIon(fIon.at(i).get()); // NOLINT
88  std::cout << " fIon name :" << fIon.at(i)->GetName() << std::endl;
89  std::cout << " fParticle name :" << fParticle.at(i)->GetName() << std::endl;
90 
91  } else if (fPType.at(i) == "Proton") {
92  std::cout << " In position " << i << " adding an : " << fPType.at(i) << std::endl;
93  run->AddNewParticle(fParticle.at(i).get()); // NOLINT
94  std::cout << " fIon name :" << fIon.at(i)->GetName() << std::endl;
95  std::cout << " fParticle name :" << fParticle.at(i)->GetName() << std::endl;
96  std::cout << fParticle.at(i)->GetName() << std::endl;
97  }
98  }
99 }
100 
101 Double_t AtTPC_Background::omega(Double_t x, Double_t y, Double_t z)
102 {
103  return sqrt(x * x + y * y + z * z - 2 * x * y - 2 * y * z - 2 * x * z);
104 }
105 
106 // Rotation of a 3D vector around an arbitrary axis
107 // Rodriges Formula
108 std::vector<Double_t>
109 AtTPC_Background::TRANSF(std::vector<Double_t> *from, std::vector<Double_t> *to, std::vector<Double_t> *vin)
110 {
111 
112  static std::vector<double> vout(3);
113  double n[3];
114  double normn, normf, normt;
115  double alpha, a, b;
116 
117  normf = sqrt(pow(from->at(0), 2) + pow(from->at(1), 2) + pow(from->at(2), 2));
118  normt = sqrt(pow(to->at(0), 2) + pow(to->at(1), 2) + pow(to->at(2), 2));
119 
120  alpha =
121  acos(((from->at(0)) * (to->at(0)) + (from->at(1)) * (to->at(1)) + (from->at(2)) * (to->at(2))) / (normf * normt));
122  a = sin(alpha);
123  b = 1.0 - cos(alpha);
124 
125  if (fabs(alpha) < 0.000001) {
126  vout.at(0) = vin->at(0);
127  vout.at(1) = vin->at(1);
128  vout.at(2) = vin->at(2);
129 
130  } else {
131  n[0] = ((from->at(1)) * (to->at(2)) - (from->at(2)) * (to->at(1)));
132  n[1] = ((from->at(2)) * (to->at(0)) - (from->at(0)) * (to->at(2)));
133  n[2] = ((from->at(0)) * (to->at(1)) - (from->at(1)) * (to->at(0)));
134 
135  // std::cout<<from->at(0)<<" "<<from->at(1)<<" "<<from->at(2)<<std::endl;
136  // std::cout<<to->at(0)<<" "<<to->at(1)<<" "<<to->at(2)<<std::endl;
137  // std::cout<<vin->at(0)<<" "<<vin->at(1)<<" "<<vin->at(2)<<std::endl;
138  normn = sqrt(pow(n[0], 2) + pow(n[1], 2) + pow(n[2], 2));
139  n[0] = n[0] / normn;
140  n[1] = n[1] / normn;
141  n[2] = n[2] / normn;
142 
143  vout.at(0) = (1 - b * (pow(n[2], 2) + pow(n[1], 2))) * (vin->at(0)) +
144  (-a * n[2] + b * n[0] * n[1]) * (vin->at(1)) + (a * n[1] + b * n[0] * n[2]) * (vin->at(2));
145  vout.at(1) = (a * n[2] + b * n[0] * n[1]) * (vin->at(0)) +
146  (1 - b * (pow(n[2], 2) + pow(n[0], 2))) * (vin->at(1)) +
147  (-a * n[0] + b * n[1] * n[2]) * (vin->at(2));
148  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)) +
149  (1 - b * (pow(n[1], 2) + pow(n[0], 2))) * (vin->at(2));
150  }
151 
152  return vout;
153 }
154 
155 Double_t *AtTPC_Background::TwoB(Double_t m1b, Double_t m2b, Double_t m3b, Double_t m4b, Double_t Kb, Double_t thetacm)
156 {
157 
158  static Double_t kinrec[2];
159 
160  double Et1 = Kb + m1b;
161  double Et2 = m2b;
162 
163  double s = pow(m1b, 2) + pow(m2b, 2) + 2 * m2b * Et1;
164  double t = 0.;
165  double u = 0.;
166 
167  double a = 4. * m2b * s;
168  double b = (pow(m3b, 2) - pow(m4b, 2) + s) * (pow(m1b, 2) - pow(m2b, 2) - s);
169  double c =
170  AtTPC_Background::omega(s, pow(m1b, 2), pow(m2b, 2)) * AtTPC_Background::omega(s, pow(m3b, 2), pow(m4b, 2));
171 
172  double Et3 = (c * cos((180. - thetacm) * 3.1415926535 / 180) - b) / a; // estamos viendo el recoil en cm (pi-theta)
173  // double Et4 = (Et1 + m2b - Et3);
174 
175  double K3 = Et3 - m3b;
176  // double K4 = Et4 - m4b;
177 
178  //------------------Mandestam variables
179  // t = pow(m2b, 2) + pow(m4b, 2) - 2 * m2b * Et4;
180  u = pow(m2b, 2) + pow(m3b, 2) - 2 * m2b * Et3;
181 
182  double theta_lab = acos(
183  ((s - pow(m1b, 2) - pow(m2b, 2)) * (pow(m2b, 2) + pow(m3b, 2) - u) +
184  2. * pow(m2b, 2) * (pow(m2b, 2) + pow(m4b, 2) - s - u)) /
185  (AtTPC_Background::omega(s, pow(m1b, 2), pow(m2b, 2)) * AtTPC_Background::omega(u, pow(m2b, 2), pow(m3b, 2))));
186 
187  kinrec[0] = K3;
188  kinrec[1] = theta_lab;
189 
190  // std::cout<<K3 <<" "<<K4<<" "<<theta_lab<<std::endl;
191 
192  return kinrec;
193 }
194 
195 std::vector<Double_t> AtTPC_Background::BreakUp(std::vector<Double_t> *Pdeuteron)
196 {
197 
198  static std::vector<double> Pproton(3);
199 
200  double mp = 1.0078250322 * 931.494; // proton
201  double mn = 1.0086649158 * 931.494; // neutron
202  double md = mp + mn + 1.0 * (gRandom->Uniform()); // Deuteron unbound around 1 MeV excitation energy
203  double pdL[3] = {Pdeuteron->at(0), Pdeuteron->at(1), Pdeuteron->at(2)};
204  double EdL = sqrt(pow(pdL[0], 2) + pow(pdL[1], 2) + pow(pdL[2], 2) + pow(md, 2));
205 
206  // deuteron breaks up into pn
207  double normPdL = sqrt(pow(pdL[0], 2) + pow(pdL[1], 2) + pow(pdL[2], 2));
208  double betad = normPdL / EdL;
209  double gammad = 1.0 / sqrt(1.0 - pow(betad, 2));
210  double S_pn = pow(md, 2);
211  double Pcpn = 0.5 * AtTPC_Background::omega(S_pn, pow(mp, 2), pow(mn, 2)) / sqrt(S_pn);
212 
213  //-----------generate isotropically theta and phi of particles p and n
214  double ran1 = (gRandom->Uniform());
215  double ran2 = (gRandom->Uniform());
216  double thetapn = acos(2 * ran1 - 1.);
217  double phipn = 2 * TMath::Pi() * ran2;
218 
219  // generate the p and n in the 2H rest frame
220  double pprest[3], pnrest[3];
221  pprest[2] = Pcpn * cos(thetapn);
222  pprest[0] = Pcpn * sin(thetapn) * cos(phipn);
223  pprest[1] = Pcpn * sin(thetapn) * sin(phipn);
224  pnrest[2] = -1 * pprest[2];
225  pnrest[0] = -1 * pprest[0];
226  pnrest[1] = -1 * pprest[1];
227  double Eprest = sqrt(pow(Pcpn, 2) + pow(mp, 2));
228  double Enrest = sqrt(pow(Pcpn, 2) + pow(mn, 2));
229 
230  // boost to 2He frame
231  double ppL[3], pnL[3];
232  ppL[0] = pprest[0];
233  ppL[1] = pprest[1];
234  ppL[2] = gammad * (pprest[2] + betad * Eprest);
235  // double EpL = sqrt(pow(mp, 2) + pow(ppL[0], 2) + pow(ppL[1], 2) + pow(ppL[2], 2));
236 
237  pnL[0] = pnrest[0];
238  pnL[1] = pnrest[1];
239  pnL[2] = gammad * (pnrest[2] + betad * Enrest);
240  // double EnL = sqrt(pow(mn, 2) + pow(pnL[0], 2) + pow(pnL[1], 2) + pow(pnL[2], 2));
241 
242  // rotate to the 2H direction
243  std::vector<double> fvfrom(3);
244  std::vector<double> fvto(3);
245  std::vector<double> fvin(3);
246  std::vector<double> fvout(3);
247  fvfrom.at(0) = 0;
248  fvfrom.at(1) = 0;
249  fvfrom.at(2) = 1;
250  fvto.at(0) = pdL[0];
251  fvto.at(1) = pdL[1];
252  fvto.at(2) = pdL[2];
253  fvin.at(0) = ppL[0];
254  fvin.at(1) = ppL[1];
255  fvin.at(2) = ppL[2];
256  fvout = AtTPC_Background::TRANSF(&fvfrom, &fvto, &fvin);
257  Pproton.at(0) = fvout.at(0);
258  Pproton.at(1) = fvout.at(1);
259  Pproton.at(2) = fvout.at(2);
260 
261  // std::cout<<"El breakup***************** "<<betad<<" "<<S_pn<<" "<<Pcpn<<std::endl;
262 
263  return Pproton;
264 }
265 
266 // ----- Public method ReadEvent --------------------------------------
267 Bool_t AtTPC_Background::ReadEvent(FairPrimaryGenerator *primGen)
268 {
269 
270  fIsDecay = kFALSE;
271 
272  fBeamEnergy = AtVertexPropagator::Instance()->GetEnergy();
273  std::cout << " -I- AtTPC_Background Residual energy : " << AtVertexPropagator::Instance()->GetEnergy() << std::endl;
274 
275  fPxBeam = AtVertexPropagator::Instance()->GetPx();
276  fPyBeam = AtVertexPropagator::Instance()->GetPy();
277  fPzBeam = AtVertexPropagator::Instance()->GetPz();
278 
279  // fPxBeam = fPx.at(0) ;
280  // fPyBeam = fPy.at(0) ;
281  // fPzBeam = fPz.at(0) ;
282 
283  if (fBeamEnergy == 0) {
284  std::cout << "-I- AtTP_Background : No solution!" << std::endl;
285  AtVertexPropagator::Instance()->SetValidKine(kFALSE);
286  }
287 
288  if (!AtVertexPropagator::Instance()->GetValidKine()) {
289 
290  fPx.at(2) = 0.; // To GeV for FairRoot
291  fPy.at(2) = 0.;
292  fPz.at(2) = 0.;
293 
294  fPx.at(3) = 0.;
295  fPy.at(3) = 0.;
296  fPz.at(3) = 0.;
297 
298  fPx.at(4) = 0.;
299  fPy.at(4) = 0.;
300  fPz.at(4) = 0.;
301 
302  fPx.at(5) = 0.;
303  fPy.at(5) = 0.;
304  fPz.at(5) = 0.;
305  }
306 
307  m1 = Masses.at(0) * amu + fExEnergy.at(0);
308  m2 = Masses.at(1) * amu + fExEnergy.at(1);
309  m3 = Masses.at(2) * amu; // recoil 1
310  m4 = Masses.at(3) * amu + fExEnergy.at(3); // ejectile 1
311  m7 = Masses.at(4) * amu; // recoil 2
312  m8 = Masses.at(5) * amu + fExEnergy.at(5); // ejectile 2
313  K1 = sqrt(pow(fPx.at(0), 2) + pow(fPy.at(0), 2) + pow(fPz.at(0), 2) + pow(m1, 2)) - m1;
314  // K1 = sqrt(pow(fPxBeam*1000.0,2) + pow(fPyBeam*1000.0,2) + pow( fPzBeam*1000.0,2) + pow(m1,2)) - m1;
315 
316  Double_t fThetaCmsMin = 0.;
317  Double_t fThetaCmsMax = 8;
318 
319  Double_t costhetamin = TMath::Cos(fThetaCmsMin * TMath::DegToRad());
320  Double_t costhetamax = TMath::Cos(fThetaCmsMax * TMath::DegToRad());
321 
322  /*
323  //proton 1 from (d,p)
325  Double_t thetacmsInput = TMath::ACos( (costhetamax - costhetamin )*gRandom->Uniform() + costhetamin
326 )*TMath::RadToDeg(); Double_t* kin2B1 = AtTPC_Background::TwoB(m1, m2, m3, m4, K1, thetacmsInput); Double_t phi1 =
327 2*TMath::Pi() * gRandom->Uniform(); //flat probability in phi Double_t krec = *(kin2B1+0); Double_t angrec =
328 *(kin2B1+1); Prec = sqrt( pow(krec,2) + 2*krec*m3); fPx.at(2) = (Prec*sin(angrec)*cos(phi1) )/1000.0; // To GeV for
329 FairRoot fPy.at(2) = (Prec*sin(angrec)*sin(phi1) )/1000.0; // To GeV for FairRoot fPz.at(2) = (Prec*cos(angrec)
330 )/1000.0; // To GeV for FairRoot
331  //std::cout<<"Kin 1 "<< angrec<<" "<<krec<<std::endl;
332 
333 
334 
335  //proton 2 from (d,p)
337  thetacmsInput = TMath::ACos( (costhetamax - costhetamin )*gRandom->Uniform() + costhetamin )*TMath::RadToDeg();
338  Double_t* kin2B2 = AtTPC_Background::TwoB(m1, m2, m3, m4, K1, thetacmsInput);
339  Double_t phi2 = 2*TMath::Pi() * gRandom->Uniform(); //flat probability in phi
340  krec = *(kin2B2+0);
341  angrec = *(kin2B2+1);
342  Prec = sqrt( pow(krec,2) + 2*krec*m3);
343  fPx.at(3) = (Prec*sin(angrec)*cos(phi2) )/1000.0; // To GeV for FairRoot
344 fPy.at(3) = (Prec*sin(angrec)*sin(phi2) )/1000.0; // To GeV for FairRoot
345 fPz.at(3) = (Prec*cos(angrec) )/1000.0; // To GeV for FairRoot
346  //std::cout<<"Kin 2 "<< angrec<<" "<<krec<<std::endl;
347  */
348 
349  // proton 1 from breakup
351  Double_t thetacmsInput =
352  TMath::ACos((costhetamax - costhetamin) * gRandom->Uniform() + costhetamin) * TMath::RadToDeg();
353  Double_t *kin2B1 = AtTPC_Background::TwoB(m1, m2, m7, m8, K1, thetacmsInput);
354  Double_t phi1 = 2 * TMath::Pi() * gRandom->Uniform(); // flat probability in phi
355  Double_t krec = *(kin2B1 + 0);
356  Double_t angrec = *(kin2B1 + 1);
357  Prec = sqrt(pow(krec, 2) + 2 * krec * m2);
358  std::vector<double> Pdeut(3);
359  Pdeut.at(0) = (Prec * sin(angrec) * cos(phi1));
360  Pdeut.at(1) = (Prec * sin(angrec) * sin(phi1));
361  Pdeut.at(2) = (Prec * cos(angrec));
362  std::vector<double> Pprot1 = AtTPC_Background::BreakUp(&Pdeut);
363  fPx.at(2) = (Pprot1.at(0)) / 1000.0; // To GeV for FairRoot
364  fPy.at(2) = (Pprot1.at(1)) / 1000.0; // To GeV for FairRoot
365  fPz.at(2) = (Pprot1.at(2)) / 1000.0; // To GeV for FairRoot
366  // std::cout<<"Kin 1 "<< acos(Pprot1.at(2)/( sqrt( pow(Pprot1.at(0),2) + pow(Pprot1.at(1),2) + pow(Pprot1.at(2),2))
367  // ))*180.0/3.1415<<" "<<sqrt( pow(Pprot1.at(0),2) + pow(Pprot1.at(1),2) + pow(Pprot1.at(2),2) + pow(m3,2) )
368  // -m3<<std::endl;
369 
370  // proton 2 from breakup
372  thetacmsInput = TMath::ACos((costhetamax - costhetamin) * gRandom->Uniform() + costhetamin) * TMath::RadToDeg();
373  Double_t *kin2B2 = AtTPC_Background::TwoB(m1, m2, m7, m8, K1, thetacmsInput);
374  Double_t phi2 = 2 * TMath::Pi() * gRandom->Uniform(); // flat probability in phi
375  krec = *(kin2B2 + 0);
376  angrec = *(kin2B2 + 1);
377  Prec = sqrt(pow(krec, 2) + 2 * krec * m2);
378  Pdeut.at(0) = (Prec * sin(angrec) * cos(phi2));
379  Pdeut.at(1) = (Prec * sin(angrec) * sin(phi2));
380  Pdeut.at(2) = (Prec * cos(angrec));
381  std::vector<double> Pprot2 = AtTPC_Background::BreakUp(&Pdeut);
382  fPx.at(3) = (Pprot2.at(0)) / 1000.0; // To GeV for FairRoot
383  fPy.at(3) = (Pprot2.at(1)) / 1000.0; // To GeV for FairRoot
384  fPz.at(3) = (Pprot2.at(2)) / 1000.0; // To GeV for FairRoot
385  // std::cout<<"Kin 2 "<< acos(Pprot2.at(2)/( sqrt( pow(Pprot2.at(0),2) + pow(Pprot2.at(1),2) + pow(Pprot2.at(2),2))
386  // ))*180.0/3.1415<<" "<<sqrt( pow(Pprot2.at(0),2) + pow(Pprot2.at(1),2) + pow(Pprot2.at(2),2) + pow(m3,2) )
387  // -m3<<std::endl;
388 
389  // std::cout<<"Kin 2 "<< angrec<<" "<<krec<<std::endl;
390 
391  // proton 3 from breakup
393  thetacmsInput = TMath::ACos((costhetamax - costhetamin) * gRandom->Uniform() + costhetamin) * TMath::RadToDeg();
394  Double_t *kin2B3 = AtTPC_Background::TwoB(m1, m2, m7, m8, K1, thetacmsInput);
395  Double_t phi3 = 2 * TMath::Pi() * gRandom->Uniform(); // flat probability in phi
396  krec = *(kin2B3 + 0);
397  angrec = *(kin2B3 + 1);
398  Prec = sqrt(pow(krec, 2) + 2 * krec * m2);
399  Pdeut.at(0) = (Prec * sin(angrec) * cos(phi3));
400  Pdeut.at(1) = (Prec * sin(angrec) * sin(phi3));
401  Pdeut.at(2) = (Prec * cos(angrec));
402  std::vector<double> Pprot3 = AtTPC_Background::BreakUp(&Pdeut);
403  fPx.at(4) = (Pprot3.at(0)) / 1000.0; // To GeV for FairRoot
404  fPy.at(4) = (Pprot3.at(1)) / 1000.0; // To GeV for FairRoot
405  fPz.at(4) = (Pprot3.at(2)) / 1000.0; // To GeV for FairRoot
406  // std::cout<<"Kin 3 "<< acos(Pprot3.at(2)/( sqrt( pow(Pprot3.at(0),2) + pow(Pprot3.at(1),2) + pow(Pprot3.at(2),2))
407  // ))*180.0/3.1415<<" "<<sqrt( pow(Pprot3.at(0),2) + pow(Pprot3.at(1),2) + pow(Pprot3.at(2),2) + pow(m3,2) )
408  // -m3<<std::endl;
409 
410  // proton 4 from breakup
412  thetacmsInput = TMath::ACos((costhetamax - costhetamin) * gRandom->Uniform() + costhetamin) * TMath::RadToDeg();
413  Double_t *kin2B4 = AtTPC_Background::TwoB(m1, m2, m7, m8, K1, thetacmsInput);
414  Double_t phi4 = 2 * TMath::Pi() * gRandom->Uniform(); // flat probability in phi
415  krec = *(kin2B4 + 0);
416  angrec = *(kin2B4 + 1);
417  Prec = sqrt(pow(krec, 2) + 2 * krec * m2);
418  Pdeut.at(0) = (Prec * sin(angrec) * cos(phi4));
419  Pdeut.at(1) = (Prec * sin(angrec) * sin(phi4));
420  Pdeut.at(2) = (Prec * cos(angrec));
421  std::vector<double> Pprot4 = AtTPC_Background::BreakUp(&Pdeut);
422  fPx.at(5) = (Pprot4.at(0)) / 1000.0; // To GeV for FairRoot
423  fPy.at(5) = (Pprot4.at(1)) / 1000.0; // To GeV for FairRoot
424  fPz.at(5) = (Pprot4.at(2)) / 1000.0; // To GeV for FairRoot
425  // std::cout<<"Kin 4 "<< acos(Pprot4.at(2)/( sqrt( pow(Pprot4.at(0),2) + pow(Pprot4.at(1),2) + pow(Pprot4.at(2),2))
426  // ))*180.0/3.1415<<" "<<sqrt( pow(Pprot4.at(0),2) + pow(Pprot4.at(1),2) + pow(Pprot4.at(2),2) + pow(m3,2) )
427  // -m3<<std::endl;
428 
429  do {
430  // random_z = 100.0*(gRandom->Uniform()); //cm
431  random_r = 1.0 * (gRandom->Gaus(0, 1)); // cm
432  random_phi = 2.0 * TMath::Pi() * (gRandom->Uniform()); // rad
433 
434  } while (fabs(random_r) > 4.7); // cut at 2 sigma
435 
436  for (Int_t i = 0; i < fMult; i++) {
437 
438  int pdgType = 2212;
439 
440  fVx = random_r * cos(random_phi);
441  fVy = random_r * sin(random_phi);
442  fVz = 100.0 * (gRandom->Uniform()); // cm
443 
444  if (i > 1 && AtVertexPropagator::Instance()->GetDecayEvtCnt() && pdgType == 2212) {
445  // TODO: Dirty way to propagate only the products (0 and 1 are beam and target respectively)
446 
447  // std::cout << "-I- FairIonGenerator: Generating ions of type "
448  //<< fIon.at(i)->GetName() << " (PDG code " << pdgType << ")" << std::endl;
449  std::cout << " Momentum (" << fPx.at(i) << ", " << fPy.at(i) << ", " << fPz.at(i) << ") Gev from vertex ("
450  << fVx << ", " << fVy << ", " << fVz << ") cm" << std::endl;
451 
452  primGen->AddTrack(pdgType, fPx.at(i), fPy.at(i), fPz.at(i), fVx, fVy, fVz);
453  }
454  }
455 
456  AtVertexPropagator::Instance()
457  ->IncDecayEvtCnt(); // TODO: Okay someone should put a more suitable name but we are on a hurry...
458 
459  return kTRUE;
460 }
461 
462 ClassImp(AtTPC_Background)
AtTPC_Background::BreakUp
virtual std::vector< Double_t > BreakUp(std::vector< Double_t > *Pdeuteron)
Definition: AtTPC_Background.cxx:195
AtTPC_Background::AtTPC_Background
AtTPC_Background()
Definition: AtTPC_Background.cxx:24
AtTPC_Background::ReadEvent
virtual Bool_t ReadEvent(FairPrimaryGenerator *primGen)
Definition: AtTPC_Background.cxx:267
AtTPC_Background::TRANSF
virtual std::vector< Double_t > TRANSF(std::vector< Double_t > *from, std::vector< Double_t > *to, std::vector< Double_t > *vin)
Definition: AtTPC_Background.cxx:109
AtVertexPropagator::SetValidKine
void SetValidKine(Bool_t val)
Definition: AtVertexPropagator.cxx:325
AtVertexPropagator::GetPx
Double_t GetPx()
Definition: AtVertexPropagator.cxx:220
AtTPC_Background::omega
virtual Double_t omega(Double_t x, Double_t y, Double_t z)
Definition: AtTPC_Background.cxx:101
AtVertexPropagator::GetPz
Double_t GetPz()
Definition: AtVertexPropagator.cxx:228
AtVertexPropagator::GetPy
Double_t GetPy()
Definition: AtVertexPropagator.cxx:224
ClassImp
ClassImp(AtFindVertex)
AtVertexPropagator::GetDecayEvtCnt
Int_t GetDecayEvtCnt()
Definition: AtVertexPropagator.cxx:192
amu
constexpr float amu
Definition: AtTPC_Background.cxx:20
y
const double * y
Definition: lmcurve.cxx:20
AtVertexPropagator::GetEnergy
Double_t GetEnergy()
Definition: AtVertexPropagator.cxx:232
AtVertexPropagator.h
AtVertexPropagator::Instance
static AtVertexPropagator * Instance()
Definition: AtVertexPropagator.cxx:18
AtTPC_Background.h
AtVertexPropagator::IncDecayEvtCnt
void IncDecayEvtCnt()
Definition: AtVertexPropagator.cxx:321
c
constexpr auto c
Definition: AtRadialChargeModel.cxx:20
AtTPC_Background
Definition: AtTPC_Background.h:25
AtTPC_Background::TwoB
virtual Double_t * TwoB(Double_t m1b, Double_t m2b, Double_t m3b, Double_t m4b, Double_t Kb, Double_t thetacm)
Definition: AtTPC_Background.cxx:155
Generated on Sat Nov 11 2023 15:30:37 for ATTPCROOT by   doxygen 1.8.18