EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ConvertStarlightAsciiToTree.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ConvertStarlightAsciiToTree.C
1 // This macro reads a starlight output file (default name slight.out) and creates a root file
2 // with TLorentzVectors for the parent and a TClonesArray of TLorentzVectors for the daughter
3 // particles. The output is stored in a root file (default name starlight.root) with one branch
4 // labeled "parent" and the other labeled "daughters". Any number of daughter tracks can be
5 // accomodated. Daughter species currently accomodated are: electrons, muons, charged or neutral
6 // pions, charged or neutral kaons, and protons.
7 //
8 // To use this macro, open root and then
9 // type .x convertStarlightAsciiToTree.C("inputfilename", "outputfilename")
10 //
11 // The root file produced can be examined in a root TBrowser.
12 //
13 // A macro to read this root file and make some standard plots is also provided. This macro is
14 // called AnalyzeTree.cxx; it can be compiled and run with the anaTree.C macro by opening root
15 // and typing .x anaTree.C()
16 
17 #include <iostream>
18 #include <fstream>
19 #include <sstream>
20 #include <string>
21 
22 #include "TLorentzVector.h"
23 #include "TClonesArray.h"
24 #include "TTree.h"
25 #include "TFile.h"
26 
27 
28 
29 using namespace std;
30 double IDtoMass(int particleCode);
31 
32 struct source_electron
33 {
34  TLorentzVector* p;
35  float gamma_mass;
36  void set(float px, float py, float pz, float E, float Q2){
37  p = new TLorentzVector(px,py,pz,E);
38  gamma_mass = Q2;
39  };
40 };
41 TLorentzVector* make_beam_particle(double lorentz, int Z, int A,int dir){
42  double mass;
43  if( std::abs(A) == 0 )
44  {cout << "Considering the electron mass." <<endl;
45  mass = 0.000510998928; //Electron GeV/c^2
46  }
47  else{
48  mass = A*0.938272046; //A*proton GeV/c^2
49  cout << "Considering the proton mass."<<endl; }
50  double E = lorentz*mass;
51  double pz = std::sqrt(E*E-mass*mass);
52  cout<<"Mass = "<<mass<<endl;
53  TLorentzVector* to_ret = new TLorentzVector(0,0,dir*pz, E);
54  return to_ret;
55 
56 }
57 
58 
59 void ConvertStarlightAsciiToTree(const char* inFileName = "slight.out",
60  const char* outFileName = "starlight.root")
61 {
62 
63  // create output tree
64  TFile* outFile = new TFile(outFileName, "RECREATE");
65  if (!outFile) {
66  cerr << " error: could not create output file '" << outFileName << "'" << endl;
67  return;
68  }
69  TTree* outTree = new TTree("starlightTree", "starlightTree");
70  TLorentzVector* parentParticle = new TLorentzVector();
71  TClonesArray* daughterParticles = new TClonesArray("TLorentzVector");
72  TLorentzVector* source = new TLorentzVector();
73  TLorentzVector* target = new TLorentzVector();
74  //
75  TLorentzVector* beam1 = new TLorentzVector();
76  TLorentzVector* beam2 = new TLorentzVector();
77  std::map<string, int> set_up;
78  double photon_setup[5];
79  //
80  Float_t q2, W, Egamma, vertex_t;
81  outTree->Branch("beam1","TLorentzVector", &beam1, 32000, -1);
82  outTree->Branch("beam2","TLorentzVector", &beam2, 32000, -1);
83  //
84  outTree->Branch("Egamma", &Egamma, "Egamma/F");
85  outTree->Branch("Q2", &q2, "q2/F");
86  outTree->Branch("W", &W, "W/F");
87  outTree->Branch("t", &vertex_t, "vertex_t/F");
88  outTree->Branch("Target","TLorentzVector", &target, 32000, -1);
89  outTree->Branch("source", "TLorentzVector", &source, 32000, -1);
90  outTree->Branch("parent", "TLorentzVector", &parentParticle, 32000, -1);
91  outTree->Branch("daughters", "TClonesArray", &daughterParticles, 32000, -1);
92  //
93  ifstream inFile;
94  inFile.open(inFileName);
95  unsigned int countLines = 0;
96  int tot_events=0;
97  bool loaded_head = false;
98  while (inFile.good()) {
99  string line;
100  string label;
101  stringstream lineStream;
102  // event simulation header
103  if( loaded_head == false){
104  if( !getline(inFile, line))
105  break;
106  ++countLines;
107  lineStream.str(line);
108  cout<<lineStream.str()<<endl;
109  R__ASSERT( lineStream >> label >> set_up["prod_id"] >> set_up["part_id"] >> set_up["nevents"]
110  >> set_up["qc"] >> set_up["impulse"] >> set_up["rnd_seed"] );
111  R__ASSERT(label == "CONFIG_OPT:");
112  // - Cout the set up for user
113  cout << " -------------------- Simulation set up -------------------- "<<endl;
114  cout << "prod_id: " << set_up["prod_id"] << "\t part_id: " << set_up["part_id"] << "\t nevents: " << set_up["nevents"] << endl;
115  cout << "q_glauber: "<< set_up["qc"] << "\t impulse: "<< set_up["impulse"] << "\t rnd_seed: "<< set_up["rnd_seed"] << endl;
116  cout << " ___________________________________________________________ "<<endl;
117  //
118  lineStream.clear();
119  // beam 1 and 2
120  int Z,A;
121  double lorentz;
122  if( !getline(inFile, line))
123  break;
124  ++countLines;
125  lineStream.str(line);
126  R__ASSERT( lineStream >> label >> lorentz );
127  R__ASSERT(label == "ELECTRON_BEAM:" );
128  beam1 = make_beam_particle(lorentz, Z, A, 1);
129 
130  lineStream.clear();
131  //
132  if( !getline(inFile, line))
133  break;
134  ++countLines;
135  lineStream.str(line);
136  R__ASSERT( lineStream >> label >> Z >> A >> lorentz );
137  R__ASSERT(label == "TARGET_BEAM:" );
138  beam2 = make_beam_particle(lorentz, Z, A, -1);
139 
140  lineStream.clear();
141  // Photon set-up
142  if( !getline(inFile, line))
143  break;
144  ++countLines;
145  lineStream.str(line);
146  int g_bins, q_bins, fixed_q;
147  double q2_min, q2_max;
148  R__ASSERT( lineStream >> label >> g_bins >> fixed_q>> q_bins
149  >> q2_min>> q2_max);
150  R__ASSERT(label == "PHOTON:");
151  loaded_head = true;
152  }
153  lineStream.clear();
154  // read EVENT
155  int eventNmb, nmbTracks;
156  if (!getline(inFile, line))
157  break;
158  ++countLines;
159  ++tot_events;
160  lineStream.str(line);
161  R__ASSERT(lineStream >> label >> eventNmb >> nmbTracks);
162  if (!(label == "EVENT:"))
163  continue;
164 
165  lineStream.clear();
166 
167  // read vertex
168  if (!getline(inFile, line))
169  break;
170  ++countLines;
171  lineStream.str(line);
172  R__ASSERT(lineStream >> label);
173  R__ASSERT(label == "VERTEX:");
174 
175  *parentParticle = TLorentzVector(0, 0, 0, 0);
176 
177  lineStream.clear();
178 
179  //read gamma
180  if( !getline(inFile,line))
181  break;
182  lineStream.str(line);
183  //cout<<line.c_str()<<endl;
184  R__ASSERT(lineStream >> label >> Egamma >> q2 );
185  R__ASSERT(label == "GAMMA:");
186  lineStream.clear();
187  // read t
188  if( !getline(inFile,line))
189  break;
190  lineStream.str(line);
191  //cout<<line.c_str()<<endl;
192  R__ASSERT(lineStream >> label >> vertex_t );
193  R__ASSERT(label == "t:");
194  lineStream.clear();
195  // read target
196  if(!getline(inFile, line))
197  break;
198  ++countLines;
199  lineStream.str(line);
200  //cout<<line.c_str()<<endl;
201  float tpx=0., tpy=0., tpz=0., tE=0.;
202  R__ASSERT(lineStream >> label >> tpx >> tpy >> tpz >> tE) ;
203  R__ASSERT(label == "TARGET:");
204  *target = TLorentzVector(tpx, tpy, tpz, tE);
205 
206  lineStream.clear();
207  // read source
208  if(!getline(inFile, line))
209  break;
210  ++countLines;
211  lineStream.str(line);
212  //cout<<line.c_str()<<endl;
213  float px=0., py=0., pz=0., E=0.;
214  R__ASSERT(lineStream >> label >> px >> py >> pz >> E) ;
215  R__ASSERT(label == "SOURCE:");
216  *source = TLorentzVector(px, py, pz, E);
217 
218  lineStream.clear();
219  // Four-momentum of the gamma-ion system started as four-momentum of the target.
220  double pxtot = tpx;
221  double pytot = tpy;
222  double pztot = tpz;
223  double Etot = tE;
224  for (int i = 0; i < nmbTracks; ++i) {
225  // read tracks
226  int particleCode;
227  double momentum[3];
228  if (!getline(inFile, line))
229  break;
230  ++countLines;
231  lineStream.str(line);
232  R__ASSERT(lineStream >> label >> particleCode >> momentum[0] >> momentum[1] >> momentum[2]);
233  R__ASSERT(label == "TRACK:");
234  Double_t daughterMass = IDtoMass(particleCode);
235  if (daughterMass < 0) {break;}
236  /*Adding up the momenta and energies of every daughter particle produced to get the total four-momentum.*/
237  const double E = sqrt( momentum[0] * momentum[0] + momentum[1] * momentum[1]
238  + momentum[2] * momentum[2] + daughterMass * daughterMass);
239  new ( (*daughterParticles)[i] ) TLorentzVector(momentum[0], momentum[1], momentum[2], E);
240  *parentParticle += *(static_cast<TLorentzVector*>(daughterParticles->At(i)));
241  // Distributing the components of the four-momentum of the system to 3d momentum vector and Energy.
242  pxtot += momentum[0];
243  pytot += momentum[1];
244  pztot += momentum[2];
245  Etot += E;
246  lineStream.clear();
247  }
248  W = sqrt(Etot * Etot - pxtot * pxtot - pytot * pytot - pztot * pztot); //Lorentz invariant.
249 
250  daughterParticles->Compress();
251  outTree->Fill();
252  }
253  outTree->Write();
254  if (outFile) {
255  outFile->Close();
256  delete outFile;
257  }
258  cout<<"Processed "<<tot_events<<" events"<<endl;
259 }
260 
261  double IDtoMass(int particleCode){
262  double mass;
263  if (particleCode == 2 || particleCode==3) {mass = 0.00051099907;} // electron
264  else if (particleCode == 5 || particleCode==6) {mass = 0.105658389;} // muon
265  else if (particleCode == 8 || particleCode==9) {mass = 0.13956995;} // charged pion
266  else if (particleCode == 7) {mass = 0.1345766;} // neutral pion
267  else if (particleCode == 11|| particleCode==12) {mass = 0.493677;} // charged kaon
268  else if (particleCode == 10 || particleCode == 16) {mass = 0.497614;} // neutral kaon
269  else if (particleCode == 14) {mass = 0.93827231;} // proton
270  else if (particleCode == 1) {mass = 0;} //photon
271  else {
272  cout << "unknown daughter particle (ID = " << particleCode << "), please modify code to accomodate" << endl;
273  mass = -1.0;
274 // exit(0);
275  }
276 
277  return mass;
278  }