estarlight/utils/eTTreemaker.C

0001 // this macro takes as input an Ascii eSTARlight output
0002 // file, slight.out, and creates a standard TTREE
0003 // S. Klein, June, 2018
0004
0005
0006 #include <TFile.h>
0007 #include <TNtuple.h>
0008 #include <TMath.h>
0009 #include <cmath>
0010 #include <cassert>
0011 #include <string>
0012 #include <iostream>
0013 #include <fstream>
0014 #include <sstream>
0015
0016 using namespace std;
0017
0018 void eTTreemaker(TString slightname="slight.out")
0019 {
0020   double me=0.000511;
0021   double mfinal;
0022   // create output TTree here
0023
0024   double fspx,fspy,fspz,fse,fsrap,fspt,fsmass,p1x,p1y,p1z,p1e,p1prap;
0025   double p2x,p2y,p2z,p2e,p2prap,p1pt,p2pt;
0026   double targx, targy, targz, targE;
0027   double sx, sy, sz, sE;
0028
0029   // new variables for eSTARlight
0030   double kphoton,kperp,qsquared,xbj,ttransfer;
0031   double eTheta; //Scattered electron polar angle
0032
0033   // ttransfer will be imlemented later
0034
0035   TNtuple *esNTuple = new TNtuple("esNT","eslightNtuple","fspt:fspz:fsrap:fsmass:p1pt:p1z:p1prap:p2pt:p2z:p2prap:kg:qsq:xbj");
0036   // This is the electron, photon, final state (particle combination), particle 1 and particle 2
0037   // These are pseudorapidities, except for the final state, where it is rapidity
0038   // For the photon there is the photon energy, Q^2 the bjorken x of the target (calculated with x = Q^2/(2km_p) and t, the momentum transfer from the target
0039
0040   cout <<" Opening slight.out"<< endl;
0041
0042   ifstream inFile;
0043   cout << "FileName=" << slightname << endl;
0044   inFile.open( slightname.Data());
0045   cout << "slight.out open"<<endl;
0046   int countLines=0;
0047
0048   // Define variables for event loop
0049   int ntrk,nvtx;
0050   int i1,i2,i3,i4;  // junk variables for fscanf lines
0051
0052   string line;
0053   stringstream lineStream;
0054
0055   // event card
0056   string label;
0057   int eventNmb,nmbTracks, pcode;
0058   int nev,ntr,stopv,pdgpid1,pdgpid2;
0059
0060   //  go through input cards, which include the gammas of the electron and ion.
0061   // We expect a CONFIG_OPT, BEAM_1 (electron), BEAM_2 (ion) and PHOTON
0062   double gamma_e,gamma_ion,junk1;
0063
0064   if (!getline(inFile,line)) {cout <<" Error reading CONFIG_OPT line"<<endl;}
0065   countLines++;
0066   lineStream=stringstream(line);
0067
0068   lineStream>>label;// CONFIG_OPT
0069
0070   if (!getline(inFile,line)) {cout <<" Error reading BEAM_1 line"<<endl;}
0071   countLines++;
0072   lineStream=stringstream(line);
0073   lineStream>>label>>/*i1>>i2>>*/gamma_e;  // BEAM_1 is the electron
0074   cout<<"Electron Lorentz boost is "<<gamma_e<<endl;
0075
0076   if (!getline(inFile,line)) {cout <<" Error reading BEAM_2 line"<<endl;}
0077   countLines++;
0078   lineStream=stringstream(line);
0079   lineStream>>label>>i1>>i2>>gamma_ion;  // BEAM_2 is the ion
0080   cout << "Check... " << label << " " << i1 << " " << i2 <<  " " << gamma_ion << endl;
0081
0082   cout<<"Ion Lorentz boost is "<<gamma_ion<<endl;
0083
0084   if (!getline(inFile,line)) {cout <<" Error reading PHOTON line"<<endl;}
0085   countLines++;
0086   lineStream=stringstream(line); // PHOTON
0087
0088
0089   // begin event loop here.  eSTARlight events have the format:
0090   // Event card, Vertex card, photon card, source card, track 1 card, track 2 card
0091
0092   while (inFile.good())
0093       {
0094
0095     if (!getline(inFile,line)) {break;}
0096     countLines++;
0097     lineStream=stringstream(line);
0098     lineStream>>label>>eventNmb>>nmbTracks;
0099     if (!(label =="EVENT:")) continue;
0100     if (eventNmb < 5) {cout <<"Reached Event "<<eventNmb<<endl;}
0101
0102     // vertex card should be followed by gamma card, t card and source card.
0103     if (!getline(inFile,line)) {break;}
0104         countLines++;
0105
0106     lineStream=stringstream(line);
0107     lineStream>>label;
0108
0109     if (!getline(inFile,line)) {break;}
0110     lineStream=stringstream(line);
0111     lineStream>>label>>kphoton>>qsquared;
0112
0113     if (!getline(inFile,line)) {break;}
0114     lineStream=stringstream(line);
0115     lineStream>>label>>ttransfer;
0116
0117     if (!getline(inFile,line)) {break;}
0118     lineStream=stringstream(line);
0119     lineStream>>label >> targx >> targy >> targz >> targE;
0120
0121   if (!getline(inFile,line)) {break;}
0122     lineStream=stringstream(line);
0123     lineStream>>label >> sx >> sy >> sz >> sE;
0124
0125     // two track cards
0126     if (!getline(inFile,line)) {break;}
0127         countLines++;
0128
0129     lineStream=stringstream(line);
0130     lineStream>>label>>pcode>>p1x>>p1y>>p1z>>i1>>i2>>i3>>pdgpid1;
0131
0132     if (!getline(inFile,line)) {break;}
0133         countLines++;
0134     lineStream=stringstream(line);
0135
0136     lineStream>>label>>pcode>>p2x>>p2y>>p2z>>i1>>i2>>i3>>pdgpid2;
0137
0138     // get the final state masses  should be particle anti-particle, so pdgcodes should be opposite
0139     if(abs(pdgpid1)==22 && (pdgpid1 != pdgpid2)){
0140             cout<<"Error pdgpid codes don't match "<<pdgpid1<<" "<<pdgpid2<<endl;
0141         exit(-1);
0142     }
0143     else if (abs(pdgpid1)!=22 && pdgpid1 != -pdgpid2)
0144       { cout<<"Error pdgpid codes don't match "<<pdgpid1<<" "<<pdgpid2<<endl;
0145         exit(-1);
0146       }
0147
0148     double pdgabs=abs(pdgpid1);
0149     mfinal=0.;
0150     if (pdgabs == 211) mfinal=0.139; //pi+
0151     if (pdgabs == 11) mfinal=0.000511; //electron
0152     if (pdgabs == 13) mfinal=0.105; //muon
0153     if (pdgabs == 321) mfinal=0.494; //K+
0154     if (pdgabs == 22) mfinal=0; //photon
0155     else if (mfinal ==0)
0156       {cout<<" Error final mass=0.  pdgabs="<<pdgabs<<endl;
0157         exit(-1);
0158       }
0159
0160     // Now do needed kinematics computations
0161     double p1=sqrt(p1x*p1x+p1y*p1y+p1z*p1z);
0162     double e1=sqrt(p1*p1+mfinal*mfinal);
0163   p1prap=0.5*log((p1+p1z)/(p1-p1z));
0164     p1pt=sqrt(p1x*p1x+p1y*p1y);
0165
0166   double p2=sqrt(p2x*p2x+p2y*p2y+p2z*p2z);
0167     double e2=sqrt(p2*p2+mfinal*mfinal);
0168   p2prap=0.5*log((p2+p2z)/(p2-p2z));
0169     p2pt=sqrt(p2x*p2x+p2y*p2y);
0170
0171     eTheta = TMath::ATan(TMath::Sqrt(sx*sx+sy*sy)/sz);
0172
0173     // final state
0174         fspx=p1x+p2x;
0175         fspy=p1y+p2y;
0176         fspz=p1z+p2z;
0177     double fse=e1+e2;
0178
0179     // need to determine energy from pdgpid
0180         fspt=sqrt(fspx*fspx+fspy*fspy);
0181         fsmass=sqrt(fse*fse-(fspx*fspx+fspy*fspy+fspz*fspz));
0182         fsrap= 0.5*log((fse+fspz)/(fse-fspz));
0183     // done - now book NTuple
0184
0185     // to calculate X_bjorken, need photon energy in lab frame
0186     double kphotontargetframe=kphoton*gamma_ion;
0187     double mass = sqrt(fsmass*fsmass+qsquared*qsquared);
0188     double xbjorken=mass/(2.*gamma_ion*0.939)*exp(fsrap);
0189
0190     //Fill ntuple
0191     esNTuple->Fill(fspt,fspz,fsrap,fsmass,p1pt,p1z,p1prap,p2pt,p2z,p2prap,kphoton,qsquared,xbjorken);
0192
0193     if (eventNmb < 5)
0194        {
0195          cout<<fspx<<fspy<<fspz<<fse<<fsrap<<fspt<<fsmass<<fsrap<<endl;
0196          cout<<p1x<<p1y<<p1z<<p1prap<<" "<<p1<<"  pt1="<<p1pt<<endl;
0197          cout<<p2x<<p2y<<p2z<<p2prap<<" "<<p2<<"  pt2="<<p2pt<<endl;
0198        }
0199
0200       } //Done with event loop.
0201
0202    // write out NTuple
0203    TString baseName = slightname.ReplaceAll(".out",".root");
0204    TFile *NTfile = new TFile(TString::Format("ntuple_%s",baseName.Data()),"recreate");
0205
0206    esNTuple->Write();
0207    NTfile->Close();
0208
0209 }