ATTPCROOT  0.3.0-alpha
A ROOT-based framework for analyzing data from active target detectors
AtTPCIonDecay.cxx
Go to the documentation of this file.
1 /*
2 
3 */
4 
5 #include "AtTPCIonDecay.h"
6 
7 #include "AtVertexPropagator.h"
8 
9 #include <FairIon.h>
10 #include <FairLogger.h>
11 #include <FairParticle.h>
12 #include <FairPrimaryGenerator.h>
13 #include <FairRunSim.h>
14 
15 #include <TDatabasePDG.h>
16 #include <TGenPhaseSpace.h>
17 #include <TLorentzVector.h>
18 #include <TMath.h>
19 #include <TParticle.h>
20 #include <TParticlePDG.h>
21 #include <TRandom.h>
22 #include <TString.h>
23 #include <TVector3.h>
24 
25 #include <algorithm>
26 #include <cmath>
27 #include <cstdio>
28 #include <iostream>
29 #include <iterator>
30 #include <memory>
31 #include <utility>
32 
33 constexpr float amu = 931.494;
34 
35 Int_t AtTPCIonDecay::fgNIon = 0;
36 
37 constexpr auto cRED = "\033[1;31m";
38 constexpr auto cYELLOW = "\033[1;33m";
39 constexpr auto cNORMAL = "\033[0m";
40 constexpr auto cGREEN = "\033[1;32m";
41 constexpr auto cBLUE = "\033[1;34m";
42 constexpr auto cORANGEWARNING = "\033[29;5;202m";
43 constexpr auto cBLINKINGRED = "\033[32;5m";
44 
45 AtTPCIonDecay::AtTPCIonDecay()
46  : fMult(0), fPx(0.), fPy(0.), fPz(0.), fVx(0.), fVy(0.), fVz(0.), fIon(0), fParticle(0), fPType(0), fNbCases(0),
47  fSepEne(0), fMasses(0), fQ(0), fPxBeam(0.), fPyBeam(0.), fPzBeam(0.)
48 {
49  // cout << "-W- AtTPCIonGenerator: "
50  // << " Please do not use the default constructor! " << endl;
51 }
52 
53 // ----- Default constructor ------------------------------------------
54 AtTPCIonDecay::AtTPCIonDecay(std::vector<std::vector<Int_t>> *z, std::vector<std::vector<Int_t>> *a,
55  std::vector<std::vector<Int_t>> *q, std::vector<std::vector<Double_t>> *mass, Int_t ZB,
56  Int_t AB, Double_t BMass, Double_t TMass, Double_t ExEnergy, std::vector<Double_t> *SepEne)
57  : fMult(0), fNbCases(a->size()), fPx(0.), fPy(0.), fPz(0.), fVx(0.), fVy(0.), fVz(0.), fIon(0), fParticle(0),
58  fPType(0), fSepEne(0), fMasses(0), fQ(0), fPxBeam(0.), fPyBeam(0.), fPzBeam(0.)
59 {
60 
61  char buffer[40];
62 
63  fgNIon++;
64 
65  for (Int_t i = 0; i < fNbCases; i++)
66  fMult.push_back(a->at(i).size());
67  fParticle.resize(fNbCases);
68  fIon.resize(fNbCases);
69  fPType.resize(fNbCases);
70 
71  fZBeam = ZB;
72  fABeam = AB;
73  fBeamMass = BMass * amu / 1000.0;
74  fTargetMass = TMass * amu / 1000.0;
75  fSepEne = SepEne[0];
76  fMasses = mass[0];
77  fIsSequentialDecay = kFALSE;
78  fExEnergy = ExEnergy;
79 
80  FairRunSim *run = FairRunSim::Instance();
81  if (!run) {
82  std::cout << "-E- FairIonGenerator: No FairRun instantised!" << std::endl;
83  Fatal("FairIonGenerator", "No FairRun instantised!");
84  return;
85  }
86 
87  for (Int_t k = 0; k < fNbCases; k++) {
88  for (Int_t i = 0; i < fMult.at(k); i++) {
89 
90  std::unique_ptr<FairIon> IonBuff = nullptr;
91  std::unique_ptr<FairParticle> ParticleBuff = nullptr;
92  sprintf(buffer, "Product_Ion_dec%d_%d", k, i);
93  if (a->at(k).at(i) != 1) {
94  IonBuff = std::make_unique<FairIon>(buffer, z->at(k).at(i), a->at(k).at(i), q->at(k).at(i), 0.0,
95  mass->at(k).at(i) * amu / 1000.0);
96  ParticleBuff = std::make_unique<FairParticle>("dummyPart", 1, 1, 1.0, 0, 0.0, 0.0);
97  fPType.at(k).push_back("Ion");
98  run->AddNewIon(IonBuff.get());
99 
100  } else if (a->at(k).at(i) == 1 && z->at(k).at(i) == 1) {
101  IonBuff = std::make_unique<FairIon>(buffer, z->at(k).at(i), a->at(k).at(i), q->at(k).at(i), 0.0,
102  mass->at(k).at(i) * amu / 1000.0);
103  auto *kProton = new TParticle(); // NOLINT
104  kProton->SetPdgCode(2212);
105  ParticleBuff = std::make_unique<FairParticle>(2212, kProton);
106  fPType.at(k).push_back("Proton");
107 
108  } else if (a->at(k).at(i) == 1 && z->at(k).at(i) == 0) {
109  IonBuff = std::make_unique<FairIon>(buffer, z->at(k).at(i), a->at(k).at(i), q->at(k).at(i), 0.0,
110  mass->at(k).at(i) * amu / 1000.0);
111  auto *kNeutron = new TParticle(); // NOLINT
112  kNeutron->SetPdgCode(2112);
113  ParticleBuff = std::make_unique<FairParticle>(2112, kNeutron);
114  fPType.at(k).push_back("Neutron");
115  }
116  fIon.at(k).push_back(std::move(IonBuff));
117  fParticle.at(k).push_back(std::move(ParticleBuff));
118  } // for mult
119  } // for case
120 }
121 
122 // ----- Public method ReadEvent --------------------------------------
123 Bool_t AtTPCIonDecay::ReadEvent(FairPrimaryGenerator *primGen)
124 {
125 
126  Double_t ExEject = AtVertexPropagator::Instance()->GetScatterEx() / 1000.0; // in GeV
127  Bool_t IsGoodCase = kFALSE;
128  fIsDecay = kFALSE;
129  Double_t excitationEnergy = 0.0;
130 
131  LOG(INFO) << cBLUE << " AtTPCIonDecay - Decay energy - Excitation energy from reaction : " << ExEject
132  << " and from task : " << fExEnergy
133  << ". Beam energy : " << AtVertexPropagator::Instance()->GetEnergy() / 1000.0 << " GeV . Is Sequential? "
134  << fIsSequentialDecay << cNORMAL << "\n";
135  LOG(INFO) << cORANGEWARNING
136  << " AtTPCIonDecay - Warning: Temporary warning message to control the flow of generators.Please, check "
137  "that if the decay comes from beam fusion, the energy from reaction is 0"
138  << cNORMAL << "\n";
139 
140  if (ExEject > 0.0 && !fIsSequentialDecay) {
141  LOG(INFO) << cBLINKINGRED
142  << " AtTPCIonDecay - Warning, Inconsistent variables: Recoil excitation energy from Vertex propagator "
143  "greater than 0 but sequential decay not enabled! Continue at your own risk!"
144  << cNORMAL << "\n";
145 
146  } else if (fIsSequentialDecay && fExEnergy > 0.0) {
147 
148  LOG(INFO) << cBLINKINGRED
149  << " AtTPCIonDecay - Warning, Inconsistent variables: Sequential decay should take the Ex energy from "
150  "the reaction generator! Continue at your own risk!"
151  << cNORMAL << "\n";
152 
153  } else if (ExEject > 0.0 && fExEnergy > 0.0) {
154  LOG(INFO)
155  << cBLINKINGRED
156  << " AtTPCIonDecay - Warning, Inconsistent variables: Both, excitation energy from Vertex propagator and "
157  "excitation energy from task (introduced through the macro) are positive! Continue at your own risk!"
158  << cNORMAL << "\n";
159  }
160 
161  std::vector<Int_t> GoodCases;
162 
163  // Test for decay from reaction
164  for (Int_t i = 0; i < fNbCases; i++) {
165  // if(ExEject*1000.0>fSepEne.at(i)) {
166  if (ExEject > -1) { // Forcing all good cases
167  GoodCases.push_back(i);
168  IsGoodCase = kTRUE;
169  }
170  }
171  if (IsGoodCase) {
172  int RandVar = (int)(GoodCases.size()) * gRandom->Uniform();
173  auto it = GoodCases.begin();
174  std::advance(it, RandVar);
175  Int_t Case = *it;
176  // LOG(INFO)<<"iterator "<<" "<<Case<<" "<<RandVar<<" "<<GoodCases.size()<<" "<<std::endl;
177 
178  Double_t beta;
179  Double_t s = 0.0;
180  Double_t mass_1[10] = {0.0};
181  Double_t *pMass;
182  Double_t M_tot = 0;
183  TLorentzVector fEnergyImpulsionLab_beam;
184  TLorentzVector fEnergyImpulsionLab_Total;
185  TLorentzVector fEnergyImpulsionLab_target;
186  TLorentzVector fEnergyImpulsionFinal;
187 
188  TVector3 fImpulsionLab_beam;
189  std::vector<TLorentzVector *> p_vector;
190  TGenPhaseSpace event1;
191 
192  fPx.clear();
193  fPy.clear();
194  fPz.clear();
195 
196  fPx.resize(fMult.at(Case));
197  fPy.resize(fMult.at(Case));
198  fPz.resize(fMult.at(Case));
199 
200  // LOG(INFO)<<" Case : "<<Case<<" with multiplicity : "<<fMult.at(Case)<<"\n";
201 
202  fIsDecay = kFALSE;
203 
204  if (fIsSequentialDecay) // NB: Decay modelled as two-step (coming from reaction generator)
205  {
206  fBeamEnergy = AtVertexPropagator::Instance()->GetTrackEnergy(0) / 1000.0;
207  TVector3 ScatP = AtVertexPropagator::Instance()->GetScatterP();
208  fPxBeam = ScatP.X();
209  fPyBeam = ScatP.Y();
210  fPzBeam = ScatP.Z();
211  excitationEnergy = ExEject; // From ejectile
212  } else { // simultaneous
213  fBeamEnergy = AtVertexPropagator::Instance()->GetEnergy() / 1000.0;
214  fPxBeam = AtVertexPropagator::Instance()->GetPx();
215  fPyBeam = AtVertexPropagator::Instance()->GetPy();
216  fPzBeam = AtVertexPropagator::Instance()->GetPz();
217  excitationEnergy = fExEnergy; // From compound nucleus
218  }
219 
220  TParticlePDG *thisPart0 = nullptr;
221  thisPart0 = TDatabasePDG::Instance()->GetParticle(
222  fIon.at(Case).at(0)->GetName()); // NB: The first particle of the list must be the decaying ion
223  int pdgType0 = thisPart0->PdgCode();
224 
225  LOG(INFO) << cBLUE << " Ejectile info : " << pdgType0 << " " << fBeamMass << " " << fPxBeam << " " << fPyBeam
226  << " " << fPzBeam << " " << fBeamEnergy << " " << excitationEnergy << cNORMAL << "\n";
227 
228  // === Phase Space Calculation
229 
230  fImpulsionLab_beam = TVector3(fPxBeam, fPyBeam, fPzBeam);
231  fEnergyImpulsionLab_beam = TLorentzVector(fImpulsionLab_beam, fBeamMass + fBeamEnergy + ExEject);
232  fEnergyImpulsionLab_target = TLorentzVector(TVector3(0, 0, 0), fTargetMass);
233 
234  if (fTargetMass > 0 && fIsSequentialDecay) {
235  LOG(INFO)
236  << cBLINKINGRED
237  << " AtTPCIonDecay - Warning, Inconsistent variables: Target Impulsion included in sequential decay. "
238  "Continue at your own risk!"
239  << cNORMAL << "\n";
240  }
241 
242  fEnergyImpulsionLab_Total = fEnergyImpulsionLab_beam + fEnergyImpulsionLab_target;
243 
244  s = fEnergyImpulsionLab_Total.M2();
245  // beta = fEnergyImpulsionLab_Total.Beta();
246 
247  for (Int_t i = 0; i < fMult.at(Case); i++) {
248  M_tot += fMasses.at(Case).at(i) * amu / 1000.0;
249  mass_1[i] = fMasses.at(Case).at(i) * amu / 1000.0;
250  std::cout << fMasses.at(Case).at(i) * amu / 1000.0 << " " << M_tot << " " << fBeamMass << " "
251  << fBeamMass - M_tot << " " << ExEject << " " << sqrt(s) << " " << (sqrt(s) - M_tot) * 1000.0
252  << std::endl;
253  }
254 
255  // std::cout<<" S : "<<s<<" Pow(M) "<<pow(M_tot,2)<<" "<<AtVertexPropagator::Instance()->GetScatterEx()<<"
256  // "<<fSepEne.at(Case)<<std::endl;
257 
258  if (s > pow(M_tot, 2)) {
259  // if(ExEject*1000.0>fSepEne){
260  fIsDecay = kTRUE;
261  event1.SetDecay(fEnergyImpulsionLab_Total, fMult.at(Case), mass_1);
262  event1.Generate(); // to generate the random final state
263  std::vector<Double_t> KineticEnergy;
264  std::vector<Double_t> ThetaLab;
265 
266  LOG(INFO) << cBLUE << " AtTPCIonDecay - Phase Space Information "
267  << "\n";
268  for (Int_t i = 0; i < fMult.at(Case); i++) {
269  p_vector.push_back(event1.GetDecay(i));
270  fPx.at(i) = p_vector.at(i)->Px();
271  fPy.at(i) = p_vector.at(i)->Py();
272  fPz.at(i) = p_vector.at(i)->Pz();
273  KineticEnergy.push_back((p_vector.at(i)->E() - mass_1[i]) * 1000.0);
274  ThetaLab.push_back(p_vector.at(i)->Theta() * 180. / TMath::Pi());
275  LOG(INFO) << " Particle " << i << " - TKE (MeV) : " << KineticEnergy.at(i)
276  << " - Lab Angle (deg) : " << ThetaLab.at(i) << "\n";
277  }
278  LOG(INFO) << cNORMAL;
279 
280  } else { // if kinematics condition
281 
282  LOG(INFO) << cYELLOW << "AtTPCIonDecay - Warning, kinematical conditions for decay not fulfilled " << cNORMAL
283  << "\n";
284  LOG(INFO) << cYELLOW << " s = " << s << " - pow(M_tot,2) = " << pow(M_tot, 2) << cNORMAL << "\n";
285  }
286 
287  // === Propagate the decay products from the vertex of the reaction
288 
289  for (Int_t i = 0; i < fMult.at(Case); i++) {
290  TParticlePDG *thisPart = nullptr;
291 
292  if (fPType.at(Case).at(i) == "Ion")
293  thisPart = TDatabasePDG::Instance()->GetParticle(fIon.at(Case).at(i)->GetName());
294  else if (fPType.at(Case).at(i) == "Proton")
295  thisPart = TDatabasePDG::Instance()->GetParticle(fParticle.at(Case).at(i)->GetName());
296  else if (fPType.at(Case).at(i) == "Neutron")
297  thisPart = TDatabasePDG::Instance()->GetParticle(fParticle.at(Case).at(i)->GetName());
298 
299  if (!thisPart) {
300  if (fPType.at(Case).at(i) == "Ion")
301  std::cout << "-W- FairIonGenerator: Ion " << fIon.at(Case).at(i)->GetName() << " not found in database!"
302  << std::endl;
303  else if (fPType.at(Case).at(i) == "Proton")
304  std::cout << "-W- FairIonGenerator: Particle " << fParticle.at(Case).at(i)->GetName()
305  << " not found in database!" << std::endl;
306  else if (fPType.at(Case).at(i) == "Neutron")
307  std::cout << "-W- FairIonGenerator: Particle " << fParticle.at(Case).at(i)->GetName()
308  << " not found in database!" << std::endl;
309  return kFALSE;
310  }
311 
312  int pdgType = thisPart->PdgCode();
313 
314  if (AtVertexPropagator::Instance()->Getd2HeEvt()) {
315  TVector3 d2HeVtx = AtVertexPropagator::Instance()->Getd2HeVtx();
316  fVx = d2HeVtx.X();
317  fVy = d2HeVtx.Y();
318  fVz = d2HeVtx.Z();
319  } else {
320  fVx = AtVertexPropagator::Instance()->GetVx();
321  fVy = AtVertexPropagator::Instance()->GetVy();
322  fVz = AtVertexPropagator::Instance()->GetVz();
323  }
324 
325  // std::cout << "-I- FairIonGenerator: Generating " <<" with mass "<<thisPart->Mass()<<" ions of type "<<
326  // fIon.at(Case).at(i)->GetName() << " (PDG code " << pdgType << ")" << std::endl; std::cout << " Momentum
327  // (" << fPx.at(i) << ", " << fPy.at(i) << ", " << fPz.at(i)
328  // << ") Gev from vertex (" << fVx << ", " << fVy
329  // << ", " << fVz << ") cm" << std::endl;
330 
331  if (fIsDecay) {
332  primGen->AddTrack(pdgType, fPx.at(i), fPy.at(i), fPz.at(i), fVx, fVy, fVz);
333  }
334 
335  } // for fMult.at(Case)
336  } // if IsGoodCase
337 
338  // if (!fIsSequentialDecay)
339  AtVertexPropagator::Instance()->IncDecayEvtCnt(); // Increase count only if no other generator is meant to do it.
340 
341  return kTRUE;
342 }
343 
344 ClassImp(AtTPCIonDecay)
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
AtVertexPropagator::GetPz
Double_t GetPz()
Definition: AtVertexPropagator.cxx:228
AtVertexPropagator::GetPy
Double_t GetPy()
Definition: AtVertexPropagator.cxx:224
ClassImp
ClassImp(AtFindVertex)
amu
constexpr float amu
Definition: AtTPCIonDecay.cxx:33
AtVertexPropagator::GetScatterEx
Double_t GetScatterEx()
Definition: AtVertexPropagator.cxx:304
AtTPCIonDecay::AtTPCIonDecay
AtTPCIonDecay()
Definition: AtTPCIonDecay.cxx:45
cNORMAL
constexpr auto cNORMAL
Definition: AtTPCIonDecay.cxx:39
cGREEN
constexpr auto cGREEN
Definition: AtTPCIonDecay.cxx:40
AtVertexPropagator::GetVy
Double_t GetVy()
Definition: AtVertexPropagator.cxx:200
AtVertexPropagator::Getd2HeVtx
TVector3 Getd2HeVtx()
Definition: AtVertexPropagator.cxx:308
AtVertexPropagator::GetEnergy
Double_t GetEnergy()
Definition: AtVertexPropagator.cxx:232
AtVertexPropagator.h
cRED
constexpr auto cRED
Definition: AtTPCIonDecay.cxx:37
cBLUE
constexpr auto cBLUE
Definition: AtTPCIonDecay.cxx:41
AtVertexPropagator::GetScatterP
TVector3 GetScatterP()
Definition: AtVertexPropagator.cxx:300
AtVertexPropagator::Instance
static AtVertexPropagator * Instance()
Definition: AtVertexPropagator.cxx:18
AtTPCIonDecay
Definition: AtTPCIonDecay.h:19
AtTPCIonDecay.h
cBLINKINGRED
constexpr auto cBLINKINGRED
Definition: AtTPCIonDecay.cxx:43
AtVertexPropagator::GetTrackEnergy
Double_t GetTrackEnergy(int trackID)
Definition: AtVertexPropagator.cxx:109
cYELLOW
constexpr auto cYELLOW
Definition: AtTPCIonDecay.cxx:38
AtVertexPropagator::IncDecayEvtCnt
void IncDecayEvtCnt()
Definition: AtVertexPropagator.cxx:321
AtTPCIonDecay::ReadEvent
virtual Bool_t ReadEvent(FairPrimaryGenerator *primGen)
Definition: AtTPCIonDecay.cxx:123
cORANGEWARNING
constexpr auto cORANGEWARNING
Definition: AtTPCIonDecay.cxx:42
Generated on Sat Nov 11 2023 15:30:38 for ATTPCROOT by   doxygen 1.8.18