ATTPCROOT  0.3.0-alpha
A ROOT-based framework for analyzing data from active target detectors
AtTPCXSReader.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- AtTPCXSReader implementation file -----
3 // ----- Created 03/07/18 by H. Alvarez -----
4 // -------------------------------------------------------------------------
5 #include "AtTPCXSReader.h"
6 
7 #include "AtVertexPropagator.h"
8 
9 #include <FairIon.h>
10 #include <FairParticle.h>
11 #include <FairPrimaryGenerator.h>
12 #include <FairRunSim.h>
13 
14 #include <TDatabasePDG.h>
15 #include <TH2.h>
16 #include <TMath.h>
17 #include <TParticle.h>
18 #include <TParticlePDG.h>
19 #include <TRandom.h>
20 #include <TVector3.h>
21 
22 #include <algorithm>
23 #include <cmath>
24 #include <cstdio>
25 #include <cstdlib>
26 #include <fstream>
27 #include <iostream>
28 
29 using std::cout;
30 using std::endl;
31 
32 Int_t AtTPCXSReader::fgNIon = 0;
33 
34 AtTPCXSReader::AtTPCXSReader() : fMult(0), fPx(0.), fPy(0.), fPz(0.), fVx(0.), fVy(0.), fVz(0.), fIon(0), fQ(0)
35 {
36  // cout << "-W- AtTPCXSReader: "
37  // << " Please do not use the default constructor! " << endl;
38 }
39 
40 AtTPCXSReader::AtTPCXSReader(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,
42  std::vector<Double_t> *pz, std::vector<Double_t> *mass)
43  : fMult(mult), fPx(0.), fPy(0.), fPz(0.), fVx(0.), fVy(0.), fVz(0.), fIon(0), fPType(0.), fQ(0)
44 {
45 
46  fgNIon++;
47 
48  fIon.reserve(fMult);
49 
50  SetXSFileName();
51 
52  TString dir = getenv("VMCWORKDIR");
53  TString XSFileName = dir + "/AtGenerators/" + fXSFileName;
54  std::cout << " AtTPCXSReader: Opening input file " << XSFileName << std::endl;
55  auto fInputXSFile = std::ifstream(XSFileName);
56  // std::ifstream* fInputXSFile = new
57  // std::ifstream("/home/ayyadlim/fair_install/AtTPCROOTv2_HAP/AtGenerators/xs_22Mgp_fusionEvaporation.txt");
58  if (!fInputXSFile.is_open())
59  Fatal("AtTPCXSReader", "Cannot open input file.");
60 
61  std::cout << "AtTPCXSReader: opening PACE cross sections..." << std::endl;
62 
63  Double_t ene[31];
64  Double_t xs[31][18];
65 
66  // fixed format
67  for (Int_t energies = 0; energies < 31; energies++) {
68  fInputXSFile >> ene[energies];
69  // std::cout << ene[energies] << " ";
70  for (Int_t xsvalues = 0; xsvalues < 18; xsvalues++) {
71  fInputXSFile >> xs[energies][xsvalues];
72  // std::cout << xs[energies][xsvalues]<< " ";
73  }
74  // std::cout << std::endl;
75  }
76 
77  fh_pdf = new TH2F("pdf", "pdf", 31, 0, 31, 18, 0, 180); // NOLINT (ROOT will cleanup)
78  for (Int_t energies = 0; energies < 31; energies++) {
79  for (Int_t xsvalues = 0; xsvalues < 18; xsvalues++) {
80  fh_pdf->SetBinContent(energies + 1, xsvalues + 1, xs[energies][xsvalues]);
81  }
82  }
83  // fh_pdf->Write();
84 
85  auto *kProton = new TParticle(); // NOLINT Probably actually a problem though
86  kProton->SetPdgCode(2212);
87 
88  auto *kNeutron = new TParticle(); // NOLINT Probably actually a problem though
89  kNeutron->SetPdgCode(2112);
90 
91  char buffer[30];
92  for (Int_t i = 0; i < fMult; i++) {
93  fPx.push_back(Double_t(a->at(i)) * px->at(i));
94  fPy.push_back(Double_t(a->at(i)) * py->at(i));
95  fPz.push_back(Double_t(a->at(i)) * pz->at(i));
96  Masses.push_back(mass->at(i) * 1000.0);
97  fWm.push_back(mass->at(i) * 1000.0);
98  FairIon *IonBuff;
99  FairParticle *ParticleBuff;
100  sprintf(buffer, "Product_Ion%d", i);
101 
102  if (a->at(i) != 1) {
103  IonBuff = new FairIon(buffer, z->at(i), a->at(i), q->at(i), 0.0, // NOLINT Probably actually a problem though
104  mass->at(i));
105  ParticleBuff = // NOLINT Probably actually a problem though
106  new FairParticle("dummyPart", 1, 1, 1.0, 0, 0.0, 0.0);
107  fPType.emplace_back("Ion");
108  std::cout << " Adding : " << buffer << std::endl;
109  } else if (a->at(i) == 1 && z->at(i) == 1) {
110  IonBuff = new FairIon("dummyIon", 50, 50, 0, 0.0, 100); // NOLINT Probably actually a problem though
111  ParticleBuff = new FairParticle(2212, kProton); // NOLINT Probably actually a problem though
112  fPType.emplace_back("Proton");
113  } else if (a->at(i) == 1 && z->at(i) == 0) {
114  IonBuff = new FairIon("dummyIon", 50, 50, 0, 0.0, 100); // NOLINT Probably actually a problem though
115  ParticleBuff = new FairParticle(2112, kNeutron); // NOLINT Probably actually a problem though
116  fPType.emplace_back("Neutron");
117  }
118 
119  std::cout << " Z " << z->at(i) << " A " << a->at(i) << std::endl;
120  // std::cout<<buffer<<std::endl;
121  fIon.push_back(IonBuff);
122  fParticle.push_back(ParticleBuff);
123  }
124 
125  FairRunSim *run = FairRunSim::Instance();
126  if (!run) {
127  std::cout << "-E- FairIonGenerator: No FairRun instantised!" << std::endl;
128  Fatal("FairIonGenerator", "No FairRun instantised!");
129  }
130 
131  for (Int_t i = 0; i < fMult; i++) {
132  if (fPType.at(i) == "Ion") {
133  std::cout << " In position " << i << " adding an : " << fPType.at(i) << std::endl;
134  run->AddNewIon(fIon.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  } else if (fPType.at(i) == "Proton") {
138  std::cout << " In position " << i << " adding an : " << fPType.at(i) << std::endl;
139  // run->AddNewParticle(fParticle.at(i));
140  std::cout << " fIon name :" << fIon.at(i)->GetName() << std::endl;
141  std::cout << " fParticle name :" << fParticle.at(i)->GetName() << std::endl;
142  std::cout << fParticle.at(i)->GetName() << std::endl;
143  } else if (fPType.at(i) == "Neutron") {
144  std::cout << " In position " << i << " adding an : " << fPType.at(i) << std::endl;
145  // run->AddNewParticle(fParticle.at(i));
146  std::cout << " fIon name :" << fIon.at(i)->GetName() << std::endl;
147  std::cout << " fParticle name :" << fParticle.at(i)->GetName() << std::endl;
148  std::cout << fParticle.at(i)->GetName() << std::endl;
149  }
150  }
151 }
152 
153 Bool_t AtTPCXSReader::ReadEvent(FairPrimaryGenerator *primGen)
154 {
155  const Double_t rad2deg = 0.0174532925;
156 
157  std::vector<Double_t> Ang; // Lab Angle of the products
158  std::vector<Double_t> Ene; // Lab Energy of the products
159  Ang.reserve(2);
160  Ene.reserve(2);
161  fPx.clear();
162  fPy.clear();
163  fPx.clear();
164  fPx.resize(fMult);
165  fPy.resize(fMult);
166  fPx.resize(fMult);
167 
168  fBeamEnergy = AtVertexPropagator::Instance()->GetEnergy();
169 
170  // Requires a non zero vertex energy and pre-generated Beam event (not punch thorugh)
171  if (AtVertexPropagator::Instance()->GetEnergy() > 0 && AtVertexPropagator::Instance()->GetDecayEvtCnt() % 2 != 0) {
172  // proton parameters come from the XS PDF
173  Double_t energyFromPDF, thetaFromPDF;
174  fh_pdf->GetRandom2(energyFromPDF, thetaFromPDF);
175 
176  Ang.push_back(thetaFromPDF * TMath::Pi() / 180); // set angle PROTON (in rad)
177  Ene.push_back(energyFromPDF); // set energy PROTON
178 
179  fPxBeam = AtVertexPropagator::Instance()->GetPx();
180  fPyBeam = AtVertexPropagator::Instance()->GetPy();
181  fPzBeam = AtVertexPropagator::Instance()->GetPz();
182 
183  // Double_t eb = fBeamEnergy + fWm.at(0); // total (beam) projectile energy = projectile kinetic e + mass
184  Double_t pb2 = fBeamEnergy * fBeamEnergy + 2.0 * fBeamEnergy * fWm.at(0); //(beam) projectile momentum squared
185  // Double_t pb = TMath::Sqrt(pb2); //(beam)projectile momentum
186  // Double_t beta = pb/(eb+fWm.at(1)); // ??check beta of the projectile+target compound check??
187  // Double_t gamma = 1.0/sqrt(1.0-beta*beta);
188  Double_t e = fBeamEnergy + fWm.at(0) + fWm.at(1); // total energy (beam+target)
189  Double_t e_cm2 = e * e - pb2; // cm energy (beam+target) squared
190  Double_t e_cm = TMath::Sqrt(e_cm2); // cm energy (beam+target)
191  Double_t t_cm = e_cm - fWm.at(2) - fWm.at(3); // kinetic energy available for final products
192 
193  // HERE REMAINS THE CALCULAtION OF THE ANGLE OF THE SCAtTER AS A FUNCTION OF THE
194  // ANGLE AND KINETIC ENERGY OF THE RECOIL (proton)
195  // Double_t p_c
196  // tan theta_scatter = p_1*cos(theta1)*sin(theta1)/(p_A - p_1(cos(theta1))*cos(theta1));
197  // with p_1(cos(theta1))=mA*E_1
198 
199  Double_t t_scatter = 0;
200  if (t_cm - energyFromPDF > 0)
201  t_scatter = t_cm - energyFromPDF;
202  else
203  std::cout << "Kinetic Energy of scatter particle negative!" << std::endl;
204 
205  Double_t theta_scatter = 0.05; // TMath::Atan(*TMath::Sin(angleFromPDF)/) (in rad)
206 
207  Ang.push_back(theta_scatter); // set angle ION DUMMY FOR THE MOMENT!!!!!!!!!!!!!!! CHECK AND SOLVE
208  Ene.push_back(t_scatter); // set energy ION
209 
210  AtVertexPropagator::Instance()->SetTrackEnergy(0, Ene.at(0));
211  AtVertexPropagator::Instance()->SetTrackAngle(0, Ang.at(0) * 180.0 / TMath::Pi());
212 
213  AtVertexPropagator::Instance()->SetTrackEnergy(1, Ene.at(1));
214  AtVertexPropagator::Instance()->SetTrackAngle(1, Ang.at(1) * 180.0 / TMath::Pi());
215 
216  fPx.at(0) = 0.0;
217  fPy.at(0) = 0.0;
218  fPz.at(0) = 0.0;
219  fPx.at(1) = 0.0;
220  fPy.at(1) = 0.0;
221  fPz.at(1) = 0.0;
222 
223  Double_t phi1 = 0., phi2 = 0.;
224  phi1 = 2 * TMath::Pi() * gRandom->Uniform(); // flat probability in phi
225  phi2 = phi1 + TMath::Pi();
226 
227  // To MeV for Euler Transformation
228  TVector3 BeamPos(AtVertexPropagator::Instance()->GetPx() * 1000, AtVertexPropagator::Instance()->GetPy() * 1000,
229  AtVertexPropagator::Instance()->GetPz() * 1000);
230 
231  TVector3 direction1 = TVector3(sin(Ang.at(0)) * cos(phi1), sin(Ang.at(0)) * sin(phi1),
232  cos(Ang.at(0))); // recoil
233 
234  TVector3 direction2 = TVector3(sin(Ang.at(1)) * cos(phi2), sin(Ang.at(1)) * sin(phi2),
235  cos(Ang.at(1))); // scatter
236 
237  Double_t p2_recoil = Ene.at(0) * Ene.at(0) + 2.0 * Ene.at(0) * fWm.at(3);
238  Double_t p2_scatter = Ene.at(1) * Ene.at(1) + 2.0 * Ene.at(1) * fWm.at(2);
239  if (p2_recoil < 0 || p2_scatter < 0)
240  std::cout << "Particle momentum negative!" << std::endl;
241 
242  Double_t p_recoil = TMath::Sqrt(p2_recoil);
243  Double_t p_scatter = TMath::Sqrt(p2_scatter);
244 
245  fPx.at(2) = p_scatter * direction2.X() / 1000.0; // To GeV for FairRoot
246  fPy.at(2) = p_scatter * direction2.Y() / 1000.0;
247  fPz.at(2) = p_scatter * direction2.Z() / 1000.0;
248 
249  fPx.at(3) = p_recoil * direction1.X() / 1000.0;
250  fPy.at(3) = p_recoil * direction1.Y() / 1000.0;
251  fPz.at(3) = p_recoil * direction1.Z() / 1000.0;
252 
253  // Particle transport begins here
254  for (Int_t i = 0; i < fMult; i++) {
255  TParticlePDG *thisPart;
256  if (fPType.at(i) == "Ion")
257  thisPart = TDatabasePDG::Instance()->GetParticle(fIon.at(i)->GetName());
258  else if (fPType.at(i) == "Proton")
259  thisPart = TDatabasePDG::Instance()->GetParticle(fParticle.at(i)->GetName());
260  else if (fPType.at(i) == "Neutron")
261  thisPart = TDatabasePDG::Instance()->GetParticle(fParticle.at(i)->GetName());
262 
263  if (!thisPart) {
264  if (fPType.at(i) == "Ion")
265  std::cout << "-W- FairIonGenerator: Ion " << fIon.at(i)->GetName() << " not found in database!"
266  << std::endl;
267  else if (fPType.at(i) == "Proton")
268  std::cout << "-W- FairIonGenerator: Particle " << fParticle.at(i)->GetName() << " not found in database!"
269  << std::endl;
270  else if (fPType.at(i) == "Neutron")
271  std::cout << "-W- FairIonGenerator: Particle " << fParticle.at(i)->GetName() << " not found in database!"
272  << std::endl;
273  return kFALSE;
274  }
275 
276  int pdgType = thisPart->PdgCode();
277 
278  // Propagate the vertex of the previous event
279  fVx = AtVertexPropagator::Instance()->GetVx();
280  fVy = AtVertexPropagator::Instance()->GetVy();
281  fVz = AtVertexPropagator::Instance()->GetVz();
282 
283  // TODO: Dirty way to propagate only the products (0 and 1 are beam and target respectively)
284  if (i > 1 && AtVertexPropagator::Instance()->GetDecayEvtCnt() && pdgType != 1000500500 &&
285  fPType.at(i) == "Ion") {
286  std::cout << "-I- FairIonGenerator: Generating ions of type " << fIon.at(i)->GetName() << " (PDG code "
287  << pdgType << ")" << std::endl;
288  std::cout << " Momentum (" << fPx.at(i) << ", " << fPy.at(i) << ", " << fPz.at(i)
289  << ") Gev from vertex (" << fVx << ", " << fVy << ", " << fVz << ") cm" << std::endl;
290  primGen->AddTrack(pdgType, fPx.at(i), fPy.at(i), fPz.at(i), fVx, fVy, fVz);
291  } else if (i > 1 && AtVertexPropagator::Instance()->GetDecayEvtCnt() && pdgType == 2212 &&
292  fPType.at(i) == "Proton") {
293  std::cout << "-I- FairIonGenerator: Generating ions of type " << fParticle.at(i)->GetName() << " (PDG code "
294  << pdgType << ")" << std::endl;
295  std::cout << " Momentum (" << fPx.at(i) << ", " << fPy.at(i) << ", " << fPz.at(i)
296  << ") Gev from vertex (" << fVx << ", " << fVy << ", " << fVz << ") cm" << std::endl;
297  primGen->AddTrack(pdgType, fPx.at(i), fPy.at(i), fPz.at(i), fVx, fVy, fVz);
298  } else if (i > 1 && AtVertexPropagator::Instance()->GetDecayEvtCnt() && pdgType == 2112 &&
299  fPType.at(i) == "Neutron") {
300  std::cout << "-I- FairIonGenerator: Generating ions of type " << fParticle.at(i)->GetName() << " (PDG code "
301  << pdgType << ")" << std::endl;
302  std::cout << " Momentum (" << fPx.at(i) << ", " << fPy.at(i) << ", " << fPz.at(i)
303  << ") Gev from vertex (" << fVx << ", " << fVy << ", " << fVz << ") cm" << std::endl;
304  primGen->AddTrack(pdgType, fPx.at(i), fPy.at(i), fPz.at(i), fVx, fVy, fVz);
305  }
306  }
307  } // if residual energy > 0
308 
309  AtVertexPropagator::Instance()
310  ->IncDecayEvtCnt(); // TODO: Okay someone should put a more suitable name but we are on a hurry...
311 
312  return kTRUE;
313 }
314 
315 ClassImp(AtTPCXSReader)
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
AtTPCXSReader::SetXSFileName
void SetXSFileName(TString name="xs_22Mgp_fusionEvaporation.txt")
Definition: AtTPCXSReader.h:49
AtVertexPropagator::GetPz
Double_t GetPz()
Definition: AtVertexPropagator.cxx:228
AtVertexPropagator::GetPy
Double_t GetPy()
Definition: AtVertexPropagator.cxx:224
ClassImp
ClassImp(AtFindVertex)
AtVertexPropagator::GetVy
Double_t GetVy()
Definition: AtVertexPropagator.cxx:200
AtTPCXSReader.h
AtTPCXSReader::ReadEvent
virtual Bool_t ReadEvent(FairPrimaryGenerator *primGen)
Definition: AtTPCXSReader.cxx:153
AtVertexPropagator::GetEnergy
Double_t GetEnergy()
Definition: AtVertexPropagator.cxx:232
AtVertexPropagator.h
AtVertexPropagator::SetTrackEnergy
void SetTrackEnergy(int trackID, double energy)
Definition: AtVertexPropagator.cxx:92
AtTPCXSReader
Definition: AtTPCXSReader.h:24
AtVertexPropagator::SetTrackAngle
void SetTrackAngle(int trackID, double angle)
Definition: AtVertexPropagator.cxx:96
AtTPCXSReader::AtTPCXSReader
AtTPCXSReader()
Definition: AtTPCXSReader.cxx:34
AtVertexPropagator::Instance
static AtVertexPropagator * Instance()
Definition: AtVertexPropagator.cxx:18
AtVertexPropagator::IncDecayEvtCnt
void IncDecayEvtCnt()
Definition: AtVertexPropagator.cxx:321
Generated on Sat Nov 11 2023 15:30:38 for ATTPCROOT by   doxygen 1.8.18