ATTPCROOT  0.3.0-alpha
A ROOT-based framework for analyzing data from active target detectors
AtTPCIonPhaseSpace.cxx
Go to the documentation of this file.
1 #include "AtTPCIonPhaseSpace.h"
2 
3 #include "AtVertexPropagator.h"
4 
5 #include <FairIon.h>
6 #include <FairPrimaryGenerator.h>
7 #include <FairRunSim.h>
8 
9 #include <TDatabasePDG.h>
10 #include <TGenPhaseSpace.h>
11 #include <TLorentzVector.h>
12 #include <TMath.h>
13 #include <TParticlePDG.h>
14 #include <TString.h>
15 #include <TVector3.h>
16 
17 #include <algorithm>
18 #include <cmath>
19 #include <iostream>
20 
21 Int_t AtTPCIonPhaseSpace::fgNIon = 0;
22 
23 AtTPCIonPhaseSpace::AtTPCIonPhaseSpace()
24  : fMult(0), fPx(0.), fPy(0.), fPz(0.), fVx(0.), fVy(0.), fVz(0.), fIon(0), fQ(0)
25 {
26  // cout << "-W- AtTPCIonGenerator: "
27  // << " Please do not use the default constructor! " << endl;
28 }
29 
30 // ----- Default constructor ------------------------------------------
31 AtTPCIonPhaseSpace::AtTPCIonPhaseSpace(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,
34  std::vector<Double_t> *mass, Double_t ResEner, Int_t ZB, Int_t AB, Double_t PxB,
35  Double_t PyB, Double_t PzB, Double_t BMass, Double_t TMass)
36  : fMult(mult), fPx(0.), fPy(0.), fPz(0.), fVx(0.), fVy(0.), fVz(0.), fIon(0), fQ(0), fBeamEnergy_buff(ResEner),
37  fBeamMass(BMass), fTargetMass(TMass), fZBeam(ZB), fABeam(AB), fPxBeam(PxB), fPyBeam(PyB), fPzBeam(PzB)
38 {
39  fgNIon++;
40  fIon.reserve(fMult);
41 
42  for (Int_t i = 0; i < fMult; i++) {
43 
44  fPx.push_back(Double_t(a->at(i)) * px->at(i));
45  fPy.push_back(Double_t(a->at(i)) * py->at(i));
46  fPz.push_back(Double_t(a->at(i)) * pz->at(i));
47  Masses.push_back(mass->at(i));
48 
49  auto *IonBuff =
50  new FairIon(TString::Format("Product_Ion%d", i).Data(), z->at(i), a->at(i), q->at(i), 0.0, mass->at(i));
51  // FairIon *IonBuff = new FairIon(buffer, z->at(i), a->at(i), q->at(i));
52  // std::cout<<" Z "<<z->at(i)<<" A "<<a->at(i)<<std::endl;
53  // std::cout<<buffer<<std::endl;
54  std::cout << " Particle " << fMult << " mass " << IonBuff->GetMass() << std::endl;
55  fIon.push_back(IonBuff);
56  }
57 
58  FairRunSim *run = FairRunSim::Instance();
59  if (!run) {
60  std::cout << "-E- FairIonGenerator: No FairRun instantised!" << std::endl;
61  Fatal("FairIonGenerator", "No FairRun instantised!");
62  }
63 
64  for (Int_t i = 0; i < fMult; i++) {
65  run->AddNewIon(fIon.at(i)); // NOLINT
66  std::cout << " Z " << z->at(i) << " A " << a->at(i) << std::endl;
67  std::cout << fIon.at(i)->GetName() << std::endl;
68  }
69 }
70 
71 // ----- Public method ReadEvent --------------------------------------
72 Bool_t AtTPCIonPhaseSpace::ReadEvent(FairPrimaryGenerator *primGen)
73 {
74 
75  /* for(Int_t i=0; i<fMult; i++){
76 
77 
78  TParticlePDG* thisPart =
79  TDatabasePDG::Instance()->GetParticle(fIon.at(i)->GetName());
80 
81 
82  if ( ! thisPart ) {
83  std::cout << "-W- FairIonGenerator: Ion " << fIon.at(i)->GetName()
84  << " not found in database!" << std::endl;
85  return kFALSE;
86  }
87 
88  int pdgType = thisPart->PdgCode();
89 
90  // std::cout << "-I- FairIonGenerator: Generating " << fMult << " ions of type "
91  // << fIon.at(i)->GetName() << " (PDG code " << pdgType << ")" << std::endl;
92  // std::cout << " Momentum (" << fPx.at(i) << ", " << fPy.at(i) << ", " << fPz.at(i)
93  // << ") Gev from vertex (" << fVx << ", " << fVy
94  // << ", " << fVz << ") cm" << std::endl;
95 
96 
97  //primGen->AddTrack(pdgType, fPx.at(i), fPy.at(i), fPz.at(i), fVx, fVy, fVz);
98 
99  } */
100 
101  // === Phase Space Calculation
102  TLorentzVector fEnergyImpulsionLab_beam;
103  TLorentzVector fEnergyImpulsionLab_target;
104  TLorentzVector fEnergyImpulsionLab_Total;
105  TLorentzVector fEnergyImpulsionFinal;
106  TVector3 fImpulsionLab_beam;
107  TVector3 fImpulsionLab_target;
108  TLorentzVector *p1 = nullptr;
109  TLorentzVector *p2 = nullptr;
110  TLorentzVector *p3 = nullptr;
111  std::vector<TLorentzVector *> p_vector;
112  TGenPhaseSpace event1;
113 
114  fPx.clear();
115  fPy.clear();
116  fPx.clear();
117 
118  fPx.resize(fMult);
119  fPy.resize(fMult);
120  fPx.resize(fMult);
121 
122  fIsDecay = kFALSE;
123 
124  // AtVertexPropagator::Instance()->Test();
125 
126  // FairMCEventHeader* MCEventHeader = primGen->GetEvent();
127  // std::cout<<" Event ID : "<<MCEventHeader->GetRunID()<<std::cout;
128 
129  // gMC->CurrentMedium();
130  // TVgirtualMC* vMC =gMC->GetMC();
131  // if(!vMC->CurrentEvent()) std::cout<<" No events!"<<std::endl;
132 
133  // std::cout<<" Current Track Number : "<<stack->GetCurrentTrackNumber()<<std::endl;
134  // stack->Print(1);
135 
136  // TParticle* beam_part = stack->GetParticle(0);
137  /* TParticle* beam_part0 = stack->GetParticle(0);
138  std::cout<<" Beam particle 0 mass "<<beam_part0->GetMass()<<std::endl;
139  std::cout<<" Beam particle 0 Energy "<<beam_part0->Energy()<<std::endl;
140  std::cout<<" Beam particle 0 Pz "<<beam_part0->Pz()<<std::endl;
141  stack->Print(1);*/
142 
143  // std::cout<<" gMC Current Event : "<<gMC->CurrentEvent()<<std::cout;
144 
145  /* FairRootManager* ioMan = FairRootManager::Instance();
146 
147  TClonesArray* fPointArray = (TClonesArray*) ioMan->GetObject("AtMCPoint"); // TODO: Why this confusing name? It
148  should be fEventArray if(fPointArray) LOG(INFO)<<"-I- AtTPCIonPhaseSpace : AtMCPoint Array Found with size :
149  "<<fPointArray->GetSize()<<FairLogger::endl;
150  if(fPointArray->IsEmpty()) std::cout<<" AtMCPoint TClonesArray Empty !!!"<<std::endl;*/
151 
152  // AtMCPoint* SimPoint = (AtMCPoint*) fPointArray->At(0);
153  // SimPoint->GetXIn();
154 
155  // fBeamEnergy = fBeamEnergy_buff/1000.0; //GeV
156 
157  fBeamEnergy = AtVertexPropagator::Instance()->GetEnergy() / 1000.0;
158  std::cout << " Residual energy in AtTPCIonPhaseSpace : " << AtVertexPropagator::Instance()->GetEnergy() << std::endl;
159  fPxBeam = AtVertexPropagator::Instance()->GetPx();
160  fPyBeam = AtVertexPropagator::Instance()->GetPy();
161  fPzBeam = AtVertexPropagator::Instance()->GetPz();
162 
163  Double_t beta;
164  Double_t s = 0.0;
165  Double_t mass_1[10] = {0.0};
166  Double_t *pMass;
167 
168  Double_t M_tot = 0;
169 
170  /* mass_1[0] = (fIon.at(0)->GetMass());
171  mass_1[1] = (fIon.at(1)->GetMass());
172  mass_1[2] = (fIon.at(2)->GetMass());*/
173 
174  // std::cout<<" Beam Z momentum : "<<fABeam*fPzBeam<<std::endl;
175 
176  // fImpulsionLab_beam = TVector3(fABeam*fPxBeam,fABeam*fPyBeam,fABeam*fPzBeam);
177  fImpulsionLab_beam = TVector3(fPxBeam, fPyBeam, fPzBeam);
178  // fEnergyImpulsionLab_beam = TLorentzVector(fImpulsionLab_beam,9327.55/1000.0+fBeamEnergy);
179  fEnergyImpulsionLab_beam = TLorentzVector(fImpulsionLab_beam, fBeamMass + fBeamEnergy);
180 
181  // fEnergyImpulsionLab_target = TLorentzVector(TVector3(0,0,0),3728.40/1000.0);
182  fEnergyImpulsionLab_target = TLorentzVector(TVector3(0, 0, 0), fTargetMass);
183 
184  fEnergyImpulsionLab_Total = fEnergyImpulsionLab_beam + fEnergyImpulsionLab_target;
185  s = fEnergyImpulsionLab_Total.M2();
186 
187  std::cout << " fABeam : " << fABeam << " fPzBeam : " << fPzBeam << " fBeamEnergy : " << fBeamEnergy << std::endl;
188 
189  for (Int_t i = 0; i < fMult; i++) {
190 
191  M_tot += Masses.at(i) / 1000.0;
192  mass_1[i] = Masses.at(i) / 1000.0;
193  }
194 
195  /* mass_1[1] = Masses.at(1)/1000.0;
196  mass_1[2] = Masses.at(2)/1000.0;*/
197 
198  // std::cout<<" Mass 1 : "<<mass_1[0]<<" Mass 2 : "<<mass_1[1]<<" Mass 3 : "<<mass_1[2]<<std::endl;
199 
200  // std::cout<<" S : "<<s<<" Pow(M) "<<pow(mass_1[0]+mass_1[1]+mass_1[2],2)<<std::endl;
201  std::cout << " S : " << s << " Pow(M) " << pow(M_tot, 2) << std::endl;
202 
203  if (s > pow(M_tot, 2)) {
204 
205  fIsDecay = kTRUE;
206 
207  event1.SetDecay(fEnergyImpulsionLab_Total, fMult, mass_1);
208  // Double_t weight1 = event1.Generate();
209 
210  /* p1 = event1.GetDecay(0);
211  p2 = event1.GetDecay(1);
212  p3 = event1.GetDecay(2);*/
213 
214  std::vector<Double_t> KineticEnergy;
215  std::vector<Double_t> ThetaLab;
216 
217  std::cout << " ==== Phase Space Information ==== " << std::endl;
218  for (Int_t i = 0; i < fMult; i++) {
219 
220  p_vector.push_back(event1.GetDecay(i));
221  fPx.at(i) = p_vector.at(i)->Px();
222  fPy.at(i) = p_vector.at(i)->Py();
223  fPz.at(i) = p_vector.at(i)->Pz();
224  KineticEnergy.push_back((p_vector.at(i)->E() - mass_1[i]) * 1000);
225  ThetaLab.push_back(p_vector.at(i)->Theta() * 180. / TMath::Pi());
226  std::cout << " Particle " << i << " - TKE (MeV) : " << KineticEnergy.at(i)
227  << " - Lab Angle (deg) : " << ThetaLab.at(i) << std::endl;
228  }
229 
230  /* fPx.at(0) = p1->Px();
231  fPy.at(0) = p1->Py();
232  fPz.at(0) = p1->Pz();
233 
234  fPx.at(1) = p2->Px();
235  fPy.at(1) = p2->Py();
236  fPz.at(1) = p2->Pz();
237 
238  fPx.at(2) = p3->Px();
239  fPy.at(2) = p3->Py();
240  fPz.at(2) = p3->Pz();
241 
242  Double_t KineticEnergy_P1 = (p1->E() - mass_1[0])*1000; //MeV
243  Double_t ThetaLab_P1 = p1->Theta()*180./TMath::Pi();
244 
245  Double_t KineticEnergy_P2 = (p2->E() - mass_1[1])*1000; //MeV
246  Double_t ThetaLab_P2 = p2->Theta()*180./TMath::Pi();
247 
248  Double_t KineticEnergy_P3 = (p3->E() - mass_1[2])*1000; //MeV
249  Double_t ThetaLab_P3 = p3->Theta()*180./TMath::Pi();
250 
251  std::cout<<" ==== Phase Space Information ==== "<<std::endl;
252  std::cout<<" Particle 1 - TKE : "<<KineticEnergy_P1<<" Angle (Lab) : "<<ThetaLab_P1<<std::endl;
253  std::cout<<" Particle 2 - TKE : "<<KineticEnergy_P2<<" Angle (Lab) : "<<ThetaLab_P2<<std::endl;
254  std::cout<<" Particle 3 - TKE : "<<KineticEnergy_P3<<" Angle (Lab) : "<<ThetaLab_P3<<std::endl;*/
255 
256  } // if kinematics condition
257 
258  for (Int_t i = 0; i < fMult; i++) {
259 
260  TParticlePDG *thisPart = TDatabasePDG::Instance()->GetParticle(fIon.at(i)->GetName());
261 
262  if (!thisPart) {
263  std::cout << "-W- FairIonGenerator: Ion " << fIon.at(i)->GetName() << " not found in database!" << std::endl;
264  return kFALSE;
265  }
266 
267  int pdgType = thisPart->PdgCode();
268 
269  // Propagate the vertex of the previous event
270 
271  fVx = AtVertexPropagator::Instance()->GetVx();
272  fVy = AtVertexPropagator::Instance()->GetVy();
273  fVz = AtVertexPropagator::Instance()->GetVz();
274 
275  std::cout << "-I- FairIonGenerator: Generating " << fMult << " with mass " << thisPart->Mass() << " ions of type "
276  << fIon.at(i)->GetName() << " (PDG code " << pdgType << ")" << std::endl;
277  std::cout << " Momentum (" << fPx.at(i) << ", " << fPy.at(i) << ", " << fPz.at(i) << ") Gev from vertex ("
278  << fVx << ", " << fVy << ", " << fVz << ") cm" << std::endl;
279 
280  if (fIsDecay) {
281  primGen->AddTrack(pdgType, fPx.at(i), fPy.at(i), fPz.at(i), fVx, fVy, fVz);
282  }
283  }
284 
285  AtVertexPropagator::Instance()->IncDecayEvtCnt();
286 
287  return kTRUE;
288 }
289 
290 ClassImp(AtTPCIonPhaseSpace)
AtTPCIonPhaseSpace
Definition: AtTPCIonPhaseSpace.h:22
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)
AtTPCIonPhaseSpace.h
AtVertexPropagator::GetVy
Double_t GetVy()
Definition: AtVertexPropagator.cxx:200
AtVertexPropagator::GetEnergy
Double_t GetEnergy()
Definition: AtVertexPropagator.cxx:232
AtVertexPropagator.h
AtVertexPropagator::Instance
static AtVertexPropagator * Instance()
Definition: AtVertexPropagator.cxx:18
AtTPCIonPhaseSpace::ReadEvent
virtual Bool_t ReadEvent(FairPrimaryGenerator *primGen)
Definition: AtTPCIonPhaseSpace.cxx:72
AtVertexPropagator::IncDecayEvtCnt
void IncDecayEvtCnt()
Definition: AtVertexPropagator.cxx:321
AtTPCIonPhaseSpace::AtTPCIonPhaseSpace
AtTPCIonPhaseSpace()
Definition: AtTPCIonPhaseSpace.cxx:23
Generated on Sat Nov 11 2023 15:30:38 for ATTPCROOT by   doxygen 1.8.18