estarlight/utils/ConvertStarlightAsciiToTree.C

0001 // This macro reads a starlight output file (default name slight.out) and creates a root file
0002 // with  TLorentzVectors for the parent and a TClonesArray of TLorentzVectors for the daughter
0003 // particles.  The output is stored in a root file (default name starlight.root) with one branch
0004 // labeled "parent" and the other labeled "daughters". Any number of daughter tracks can be
0005 // accomodated.  Daughter species currently accomodated are:  electrons, muons, charged or neutral
0006 // pions, charged or neutral kaons, and protons.
0007 //
0008 // To use this macro, open root and then
0009 // type .x convertStarlightAsciiToTree.C("inputfilename", "outputfilename")
0010 //
0011 // The root file produced can be examined in a root TBrowser.
0012 //
0013 // A macro to read this root file and make some standard plots is also provided.  This macro is
0014 // called AnalyzeTree.cxx; it can be compiled and run with the anaTree.C macro by opening root
0015 // and typing .x anaTree.C()
0016
0017 #include <iostream>
0018 #include <fstream>
0019 #include <sstream>
0020 #include <string>
0021
0022 #include "TLorentzVector.h"
0023 #include "TClonesArray.h"
0024 #include "TTree.h"
0025 #include "TFile.h"
0026
0027
0028
0029 using namespace std;
0030 double IDtoMass(int particleCode);
0031
0032 struct source_electron
0033 {
0034   TLorentzVector* p;
0035   float gamma_mass;
0036   void set(float px, float py, float pz, float E, float Q2){
0037     p = new TLorentzVector(px,py,pz,E);
0038     gamma_mass = Q2;
0039   };
0040 };
0041 TLorentzVector* make_beam_particle(double lorentz, int Z, int A,int dir){
0042   double mass;
0043   if( std::abs(A) == 0 )
0044     {cout << "Considering the electron mass." <<endl;
0045     mass = 0.000510998928; //Electron GeV/c^2
0046   }
0047   else{
0048     mass = A*0.938272046; //A*proton GeV/c^2
0049   cout << "Considering the proton mass."<<endl; }
0050   double E = lorentz*mass;
0051   double pz = std::sqrt(E*E-mass*mass);
0052   cout<<"Mass = "<<mass<<endl;
0053   TLorentzVector* to_ret = new TLorentzVector(0,0,dir*pz, E);
0054   return to_ret;
0055
0056 }
0057
0058
0059 void ConvertStarlightAsciiToTree(const char* inFileName  = "slight.out",
0060                         const char* outFileName = "starlight.root")
0061 {
0062
0063     // create output tree
0064     TFile* outFile = new TFile(outFileName, "RECREATE");
0065     if (!outFile) {
0066         cerr << "    error: could not create output file '" << outFileName << "'" << endl;
0067         return;
0068     }
0069     TTree*          outTree           = new TTree("starlightTree", "starlightTree");
0070     TLorentzVector* parentParticle    = new TLorentzVector();
0071     TClonesArray*   daughterParticles = new TClonesArray("TLorentzVector");
0072     TLorentzVector* source            = new TLorentzVector();
0073     TLorentzVector* target            = new TLorentzVector();
0074     //
0075     TLorentzVector* beam1 = new TLorentzVector();
0076     TLorentzVector* beam2 = new TLorentzVector();
0077     std::map<string, int> set_up;
0078     double photon_setup[5];
0079     //
0080     Float_t        q2, W, Egamma, vertex_t;
0081     outTree->Branch("beam1","TLorentzVector", &beam1,            32000, -1);
0082     outTree->Branch("beam2","TLorentzVector", &beam2,            32000, -1);
0083     //
0084     outTree->Branch("Egamma",         &Egamma, "Egamma/F");
0085     outTree->Branch("Q2",         &q2, "q2/F");
0086     outTree->Branch("W",          &W, "W/F");
0087     outTree->Branch("t",         &vertex_t, "vertex_t/F");
0088     outTree->Branch("Target","TLorentzVector", &target,            32000, -1);
0089     outTree->Branch("source",    "TLorentzVector", &source,            32000, -1);
0090     outTree->Branch("parent",    "TLorentzVector", &parentParticle,    32000, -1);
0091     outTree->Branch("daughters", "TClonesArray",   &daughterParticles, 32000, -1);
0092     //
0093     ifstream inFile;
0094     inFile.open(inFileName);
0095     unsigned int countLines = 0;
0096     int tot_events=0;
0097     bool loaded_head = false;
0098     while (inFile.good()) {
0099         string line;
0100         string label;
0101         stringstream lineStream;
0102         // event simulation header
0103         if( loaded_head == false){
0104           if( !getline(inFile, line))
0105             break;
0106           ++countLines;
0107           lineStream.str(line);
0108           cout<<lineStream.str()<<endl;
0109           R__ASSERT( lineStream >> label >> set_up["prod_id"] >> set_up["part_id"] >> set_up["nevents"]
0110               >> set_up["qc"] >> set_up["impulse"] >> set_up["rnd_seed"] );
0111           R__ASSERT(label == "CONFIG_OPT:");
0112           // - Cout the set up for user
0113           cout << " -------------------- Simulation set up -------------------- "<<endl;
0114           cout << "prod_id: " << set_up["prod_id"] << "\t part_id: " << set_up["part_id"] << "\t nevents: " << set_up["nevents"] << endl;
0115           cout << "q_glauber: "<< set_up["qc"] << "\t impulse: "<<  set_up["impulse"] << "\t rnd_seed: "<< set_up["rnd_seed"] << endl;
0116           cout << " ___________________________________________________________ "<<endl;
0117             //
0118           lineStream.clear();
0119           // beam 1 and 2
0120           int Z,A;
0121           double lorentz;
0122           if( !getline(inFile, line))
0123             break;
0124           ++countLines;
0125           lineStream.str(line);
0126           R__ASSERT( lineStream >> label >> lorentz );
0127           R__ASSERT(label == "ELECTRON_BEAM:" );
0128           beam1 = make_beam_particle(lorentz, Z, A, 1);
0129
0130           lineStream.clear();
0131           //
0132           if( !getline(inFile, line))
0133             break;
0134           ++countLines;
0135           lineStream.str(line);
0136           R__ASSERT( lineStream >> label >> Z >> A >> lorentz );
0137           R__ASSERT(label == "TARGET_BEAM:" );
0138           beam2 = make_beam_particle(lorentz, Z, A, -1);
0139
0140           lineStream.clear();
0141           // Photon set-up
0142           if( !getline(inFile, line))
0143             break;
0144           ++countLines;
0145           lineStream.str(line);
0146           int g_bins, q_bins, fixed_q;
0147           double q2_min, q2_max;
0148           R__ASSERT( lineStream >> label >> g_bins >> fixed_q>> q_bins
0149               >> q2_min>> q2_max);
0150           R__ASSERT(label == "PHOTON:");
0151           loaded_head = true;
0152         }
0153         lineStream.clear();
0154         // read EVENT
0155         int    eventNmb, nmbTracks;
0156         if (!getline(inFile, line))
0157             break;
0158         ++countLines;
0159         ++tot_events;
0160         lineStream.str(line);
0161         R__ASSERT(lineStream >> label >> eventNmb >> nmbTracks);
0162         if (!(label == "EVENT:"))
0163             continue;
0164
0165         lineStream.clear();
0166
0167         // read vertex
0168         if (!getline(inFile, line))
0169             break;
0170         ++countLines;
0171         lineStream.str(line);
0172         R__ASSERT(lineStream >> label);
0173         R__ASSERT(label == "VERTEX:");
0174
0175         *parentParticle = TLorentzVector(0, 0, 0, 0);
0176
0177         lineStream.clear();
0178
0179         //read gamma
0180         if( !getline(inFile,line))
0181           break;
0182         lineStream.str(line);
0183         //cout<<line.c_str()<<endl;
0184         R__ASSERT(lineStream >> label >> Egamma >> q2 );
0185         R__ASSERT(label == "GAMMA:");
0186         lineStream.clear();
0187         // read t
0188         if( !getline(inFile,line))
0189           break;
0190         lineStream.str(line);
0191         //cout<<line.c_str()<<endl;
0192         R__ASSERT(lineStream >> label >> vertex_t );
0193         R__ASSERT(label == "t:");
0194         lineStream.clear();
0195         // read target
0196         if(!getline(inFile, line))
0197           break;
0198         ++countLines;
0199         lineStream.str(line);
0200         //cout<<line.c_str()<<endl;
0201         float tpx=0., tpy=0., tpz=0., tE=0.;
0202         R__ASSERT(lineStream >> label >> tpx >> tpy >> tpz >> tE) ;
0203         R__ASSERT(label == "TARGET:");
0204         *target = TLorentzVector(tpx, tpy, tpz, tE);
0205
0206         lineStream.clear();
0207         // read source
0208         if(!getline(inFile, line))
0209           break;
0210         ++countLines;
0211         lineStream.str(line);
0212         //cout<<line.c_str()<<endl;
0213         float px=0., py=0., pz=0., E=0.;
0214         R__ASSERT(lineStream >> label >> px >> py >> pz >> E) ;
0215         R__ASSERT(label == "SOURCE:");
0216         *source = TLorentzVector(px, py, pz, E);
0217
0218         lineStream.clear();
0219         // Four-momentum of the gamma-ion system started as four-momentum of the target.
0220         double pxtot = tpx;
0221         double pytot = tpy;
0222         double pztot = tpz;
0223         double Etot = tE;
0224         for (int i = 0; i < nmbTracks; ++i) {
0225             // read tracks
0226             int    particleCode;
0227             double momentum[3];
0228             if (!getline(inFile, line))
0229               break;
0230             ++countLines;
0231             lineStream.str(line);
0232             R__ASSERT(lineStream >> label >> particleCode >> momentum[0] >> momentum[1] >> momentum[2]);
0233             R__ASSERT(label == "TRACK:");
0234             Double_t daughterMass = IDtoMass(particleCode);
0235             if (daughterMass < 0) {break;}
0236             /*Adding up the momenta and energies of every daughter particle produced to get the total four-momentum.*/
0237             const double E = sqrt(  momentum[0] * momentum[0] + momentum[1] * momentum[1]
0238                                   + momentum[2] * momentum[2] + daughterMass * daughterMass);
0239             new ( (*daughterParticles)[i] ) TLorentzVector(momentum[0], momentum[1], momentum[2], E);
0240             *parentParticle += *(static_cast<TLorentzVector*>(daughterParticles->At(i)));
0241             // Distributing the components of the four-momentum of the system to 3d momentum vector and Energy.
0242             pxtot += momentum[0];
0243             pytot += momentum[1];
0244             pztot += momentum[2];
0245             Etot += E;
0246             lineStream.clear();
0247         }
0248         W = sqrt(Etot * Etot - pxtot * pxtot - pytot * pytot - pztot * pztot); //Lorentz invariant.
0249
0250         daughterParticles->Compress();
0251         outTree->Fill();
0252     }
0253     outTree->Write();
0254     if (outFile) {
0255         outFile->Close();
0256         delete outFile;
0257     }
0258     cout<<"Processed "<<tot_events<<" events"<<endl;
0259 }
0260
0261     double IDtoMass(int particleCode){
0262         double mass;
0263         if (particleCode == 2 || particleCode==3) {mass = 0.00051099907;} // electron
0264         else if (particleCode == 5 || particleCode==6) {mass = 0.105658389;} // muon
0265         else if (particleCode == 8 || particleCode==9)  {mass = 0.13956995;} // charged pion
0266         else if (particleCode == 7) {mass = 0.1345766;} // neutral pion
0267         else if (particleCode == 11|| particleCode==12) {mass = 0.493677;} // charged kaon
0268         else if (particleCode == 10 || particleCode == 16)  {mass = 0.497614;} // neutral kaon
0269         else if (particleCode == 14)    {mass = 0.93827231;} // proton
0270         else if (particleCode == 1) {mass = 0;} //photon
0271         else {
0272             cout << "unknown daughter particle (ID = " << particleCode << "), please modify code to accomodate" << endl;
0273             mass = -1.0;
0274 //          exit(0);
0275              }
0276
0277         return mass;
0278     }