hbcup/starlight/ana_starlight.C

0001 #include <iostream>
0002 #include <fstream>
0003 #include <TLorentzVector.h>
0004 #include <TF1.h>
0005 #include <TRandom3.h>
0006
0007 using namespace std;
0008
0009 const double E_MASS = 0.000510998950;  // electron mass [Gev]
0010
0011 TRandom3 *rand3;
0012
0013 // Return the relative resolution dp/p, for a particle with momentum p
0014 // The momentum resolution was taken from sPHENIX simulation
0015 double get_dpp( double p )
0016 {
0017   // Parametrization of the relative momentum resolution (dp/p)
0018   static TF1 *f_dppfit = new TF1("f_dppfit","[0]*x + [1]/x + [2]",0.3,50.0);
0019   f_dppfit->SetParameter(0,0.0010557);
0020   f_dppfit->SetParameter(1,0.00262218);
0021   f_dppfit->SetParameter(2,0.0090893);
0022
0023   // Return dp/p for a particle of momentum p
0024   double dpp = f_dppfit->Eval( p );
0025   return dpp;
0026 }
0027
0028 // Return the smeared inv mass, for two truth particles v1 and v2
0029 Double_t get_invmass_smeared(TLorentzVector v1, TLorentzVector v2)
0030 {
0031   // Get total momentum of particles 1 and 2
0032   double p1 = v1.P();
0033   double p2 = v2.P();
0034
0035   // Get the sigma of the momentum resolution for particle 1 (dp1)
0036   double sigma_p1 = get_dpp( p1 ) * p1;
0037
0038   // Now get random value for momentum deviation (mis-measured momentum)
0039   double dp1 = rand3->Gaus( 0, sigma_p1 );
0040
0041   // now modify the true momentum to the mis-measured momentum
0042   double dpp1 = 1.0 + dp1/p1;
0043   double px1 = v1.Px()*dpp1;
0044   double py1 = v1.Py()*dpp1;
0045   double pz1 = v1.Pz()*dpp1;
0046
0047   // Use mis-measured momentum to get mis-measured E, and generate lorentz 4-vector
0048   double E1 = sqrt( E_MASS*E_MASS + px1*px1 + py1*py1 + pz1*pz1 );
0049   v1.SetPxPyPzE( px1, py1, pz1, E1 );
0050
0051   // Now get momentum smearing of 2nd particle, exactly in the same way as the 1st particle
0052   double sigma_p2 = get_dpp( p2 ) * p2;
0053   double dp2 = rand3->Gaus( 0, sigma_p2 );
0054   double dpp2 = 1.0 + dp2/p2;
0055
0056   double px2 = v2.Px()*dpp2;
0057   double py2 = v2.Py()*dpp2;
0058   double pz2 = v2.Pz()*dpp2;
0059   double E2 = sqrt( E_MASS*E_MASS + px2*px2 + py2*py2 + pz2*pz2 );
0060   v2.SetPxPyPzE( px2, py2, pz2, E2 );
0061
0062   // Once we have both tracks with momentum resolution effects included, we sum them to get
0063   // measured parent particle 4-vector and return the mass
0064   TLorentzVector sum = v1 + v2;
0065
0066   return sum.M();
0067 }
0068
0069 void ana_starlight()
0070 {
0071     gStyle->SetOptStat(0);
0072     rand3 = new TRandom3(0);
0073
0074     const double integ_lumi = 4.5e9;   // in inverse barns
0075
0076     // Get cross-section of process
0077     TString xsection_txt = gSystem->GetFromPipe("tail -1 output.txt | awk '{print $9}'");
0078     TString xsection_units_txt = gSystem->GetFromPipe("tail -1 output.txt | awk '{print $10}'");
0079     Double_t xsection = xsection_txt.Atof();    // xsection is in barns
0080     if ( xsection_units_txt == "mb." )
0081     {
0082       xsection *= 1e-3;
0083     }
0084     else if ( xsection_units_txt == "microbarn." )
0085     {
0086       xsection *= 1e-6;
0087     }
0088     cout << "xsection (b) and integ lumi (1/b) = " << xsection << "\t" << integ_lumi << endl;
0089
0090     // open ASCII file.
0091     ifstream in;
0092     //in.open(Form("output.txt"));
0093     in.open("slight.out");
0094
0095     // Check if file is open
0096
0097     if (in.is_open())
0098     {
0099         cout << "File found\n";
0100         //in.close();
0101     }
0102     else
0103     {
0104         cout << "Error opening file";
0105     }
0106
0107     // kinematics variables
0108
0109     const Int_t kMaxTrack = 1000;
0110     Int_t ntrack;
0111
0112     Float_t px1;  // momentum of first daughter in the x-direction
0113     Float_t py1;  // momentum of first daughter in the y-direction
0114     Float_t pz1;  // momentum of first daughter in the z-direction
0115     Float_t E1;   // energy of the first daughter
0116     Float_t p1;   // momentum of the first daughter
0117     Float_t pt1;  // momentum of the first daughter
0118     Float_t me1;  // mass of the first daughter
0119     Float_t rap1; // rapidity of the first daughter
0120     Float_t eta1; // pseudorapidity of the first daughter
0121
0122     Float_t px2;  // kinematic variables for the 2nd daughter...
0123     Float_t py2;
0124     Float_t pz2;
0125     Float_t E2;
0126     Float_t p2;
0127     Float_t pt2;
0128     Float_t me2;
0129     Float_t rap2;
0130     Float_t eta2;
0131
0132     Float_t pxt;      // x-direction momentum of the parent particle
0133     Float_t pyt;      // y-direction momentum  of the parent particle
0134     Float_t pzt;      // z -direction momentum of the parent atom
0135     Float_t ptt;      // transverse momentum of the parent atom
0136     Float_t ptott;    // momentum of the parent atom
0137     Float_t E;        // total energy of the parent atom
0138     Float_t rapt;     // rapidity
0139     Float_t etat;     // pseudorapidity
0140     Float_t InvMass;  // invariant mass
0141
0142     Float_t theta;    // polar angle
0143
0144     Int_t num_of_events;
0145
0146     // Create tfile to save results
0147     TFile *savefile = new TFile("upc_starlight.root","RECREATE");
0148
0149     //// Create a TTree and branches of the events
0150     TTree * tree = new TTree("t", "Analyze.tree");
0151     tree->Branch("px",&pxt,"px/F");
0152     tree->Branch("py",&pyt,"py/F");
0153     tree->Branch("pz",&pzt,"pz/F");
0154     tree->Branch("p",&ptott,"p/F");
0155     tree->Branch("pt",&ptt,"pt/F");
0156     tree->Branch("y",&rapt,"y/F");
0157     tree->Branch("eta",&etat,"eta/F");
0158     tree->Branch("m", &InvMass, "m/F");
0159     tree->Branch("p1", &p1, "p1/F");
0160     tree->Branch("pt1", &pt1, "pt1/F");
0161     tree->Branch("eta1", &eta1, "eta1/F");
0162     tree->Branch("p2", &p2, "p2/F");
0163     tree->Branch("pt2", &pt2, "pt2/F");
0164     tree->Branch("eta2", &eta2, "eta2/F");
0165
0166
0167     //tree->Print();
0168
0169
0170     const double MAX_RAP = 4;
0171     TH1F *h_rapt = new TH1F("h_rapt", "Rapidity J/#Psi#rightarrow e^{+}e^{-}", 200, -MAX_RAP, MAX_RAP);
0172     TH1F *h_rapt_sphenix = new TH1F("h_rapt_sphenix", "Rapidity J/#Psi#rightarrow e^{+}e^{-} in sPHENIX Acceptance", 200, -MAX_RAP, MAX_RAP);
0173     TH1F *h_rap1 = new TH1F("h_rap1", "Rapidity Track 1", 200, -MAX_RAP, MAX_RAP);
0174     TH1F *h_rap2 = new TH1F("h_rap2", "Rapidity Track 2", 200, -MAX_RAP, MAX_RAP);
0175
0176     const double MAX_INVMASS =6;
0177     TH1F *h_InvMass = new TH1F("h_InvMass", "Invariant Mass", 200, 0, 6.0);
0178     TH1F *h_InvMass_smeared = new TH1F("h_InvMass_smeared", "Invariant Mass", 200, 0, 6.0);
0179     TH1F *h_InvMass_sphenix = new TH1F("h_InvMass_sphenix", "Invariant Mass of J/#Psi in sPHENIX Acceptance",200, 0, 6.0);
0180     TH1F *h_InvMass_smeared_sphenix = new TH1F("h_InvMass_smeared_sphenix", "Invariant Mass of J/#Psi in sPHENIX Acceptance",200, 0, 6.0);
0181
0182     TH1F *h_pt1 = new TH1F("h_pt1", "Transverse Momentum track 1", 200, 0, 2.);
0183     TH1F *h_pt2 = new TH1F("h_pt2", "Transverse Momentum track 2", 200, 0, 2.);
0184     TH1F *h_pt = new TH1F("h_pt","Transverse Momentum of J/Psi",200, 0, 2.0);
0185     TH1F *h_pt_sphenix = new TH1F("h_pt_sphenix", "Transverse Momentum of UPC in sPHENIX Acceptance", 300, 0.0, 2.0);
0186     TH1F *h_pz = new TH1F("h_pz", "Z-Axis Momentum of J/Psi", 100,-2,6);
0187
0188     const double MAX_ETA = 6;
0189     TH1F *h_eta1 = new TH1F("h_eta1", "PseudoRapidity Track 1", 200, -MAX_ETA, MAX_ETA);
0190     TH1F *h_eta2 = new TH1F("h_eta2", "PseudoRapidity Track 2", 200, -MAX_ETA, MAX_ETA);
0191
0192     // Read in one event from starlight output
0193     string label;
0194     int nev, ntr, stopv; //nev is the event number; ntr is the number of tracks within the vertex; stopv is the vertex number where the track ends
0195     Float_t v0[4]; // 0-2 = x,y,z 3=t
0196     int nv, nproc, nparent, ndaughter;//
0197     cout << label << "\t" << ndaughter << endl;
0198     const int MAXDAUGHTERS = 10; //changing this number does not change the energy or value outputted at the end of the function
0199     int gpid[MAXDAUGHTERS], pdgpid[MAXDAUGHTERS]; //gpid is the monte carlo particle id code
0200     Float_t mom[MAXDAUGHTERS][3];
0201     int  nevtrk, stopvtrk, junk[9];//number of event tracks, number of stopping tracks
0202     /* for (int idtr=0; idtr<ndaughter; idtr++)
0203        {
0204        in >> label >> gpid[idtr] >> mom[idtr][0] >> mom[idtr][1] >> mom[idtr][2]
0205        >> nevtrk >> ntr >> stopvtrk >> pdgpid[idtr];
0206        }
0207        */
0208     TLorentzVector vp; //parent particle JPsi
0209     TLorentzVector v1; //daughter particle 1 electron
0210     TLorentzVector v2; //daughter particle 2 electron
0211
0212     /*  Double_t pxt= vp.pxt(); //truthfully idek if i need this double because I already floated pxt-E and when I run it with this commented out it works.
0213         Double_t pyt= vp.pyt(); //but I'm not getting multiple graphs.
0214         Double_t pzt= vp.pzt();
0215         Double_t E= vp.E();
0216
0217
0218         vp.SetPxPyPzE(vp->GetPx(),vp->GetPy(),vp->GetPz(),vp->GetE());
0219         v1.SetPxPyPzE(v1->GetPx(),v1->GetPy(),v1->GetPz(),v1->GetE());
0220         v2.SetPxPyPzE(v2->GetPx(),v2->GetPy(),v2->GetPz(),v2->GetE()); */
0221
0222
0223     Double_t mass = vp.M();// mass of parent
0224
0225     Double_t pt = vp.Pt();//parent pT
0226
0227     unsigned int nevents = 0;   // number of events in file
0228
0229     while ( in.good() )
0230     {
0231
0232         // read the event line
0233         in >> label >> nev >> ntr >> stopv;
0234         if ( in.eof() ) break;    // quit reading if we reached eend of file
0235
0236         // read the vertex line
0237         in >> label >> v0[0] >> v0[1] >> v0[2] >> v0[3] >> nv >> nproc >> nparent >> ndaughter;
0238
0239         // read in daughters
0240         in >> label >>gpid[0] >> mom[0][0] >> mom[0][1] >> mom[0][2]
0241             >> nevtrk >> ntr >> stopvtrk >> pdgpid[0];
0242
0243         in >> label >>gpid[1] >> mom[1][0] >> mom[1][1] >> mom[1][2]
0244             >> nevtrk >> ntr >> stopvtrk >> pdgpid[1];
0245
0246         //  cout << "\n GPID=" << gpid[idtr] << "\n Momentum 0=" << mom[idtr][0] << "\n Momentum 1=" << mom[idtr][1] << "\n Momentum 2=" << mom[idtr][2]
0247         //  << "\n # of EventTracks=" << nevtrk << "\n # of Tracks w/in Vertex=" << ntr << "\n Vertex # where Track Stops=" << stopvtrk << endl;
0248
0249         //  Calculate pT
0250         // set the values in each track
0251
0252         // get momentum from 1st daughter particle
0253         px1 = mom[0][0];
0254         py1 = mom[0][1];
0255         pz1 = mom[0][2];
0256         pt1 = TMath::Sqrt( px1*px1 + py1*py1 );
0257         p1 = TMath::Sqrt( px1*px1 + py1*py1 + pz1*pz1 );
0258
0259         // get momentum from 2nd daughter particle
0260         px2 = mom[1][0];
0261         py2 = mom[1][1];
0262         pz2 = mom[1][2];
0263         pt2 = TMath::Sqrt( px2*px2 + py2*py2 );
0264         p2 = TMath::Sqrt( px2*px2 + py2*py2 + pz2*pz2 );
0265
0266         // total momentum of the daughters
0267
0268         pxt = px1 + px2;//momentum of particle one and two in the x direction equals the total momentum in x direction
0269         pyt = py1 + py2;
0270         pzt = pz1 + pz2;
0271         //cout << pzt << endl;
0272
0273         //Calculate the tranverse momentum
0274         ptt = TMath::Sqrt(pxt*pxt + pyt*pyt); //ptt is total transverse momentum
0275         // cout << px1 << "\t" << px2 << endl;
0276         Double_t p = TMath::Sqrt(pxt*pxt + pyt*pyt + pzt*pzt); //momentum in the x,y,z plane
0277         // Fill histogram
0278         h_pt->Fill(ptt);
0279
0280         // calculate the invariant mass
0281         //  cout <<"\n"<<"InvMass  = "<< InvMass << "\t"  << endl;
0282         me1 = me2 = E_MASS;
0283
0284         E1 = TMath::Sqrt((me1 *me1) + (px1*px1) + (py1*py1) +(pz1*pz1));
0285         E2 = TMath::Sqrt((me2 *me2) + (px2*px2) + (py2*py2) +(pz2*pz2));
0286
0287         //cout <<"\n"<<"Energy = "<< E << "\t" << endl;
0288
0289         E = E1 + E2;
0290
0291         v1.SetPxPyPzE(px1, py1, pz1, E1);
0292         v2.SetPxPyPzE(px2, py2, pz2, E2);
0293
0294         //mass should always be about 3.1 GeV^2 no matter the frame
0295         InvMass =TMath::Sqrt((E * E)- (pxt * pxt) -(pyt * pyt) -(pzt*pzt));
0296         Double_t m = TMath::Sqrt(pxt*pxt + pyt*pyt + pzt*pzt);
0297
0298         Double_t InvMass_smeared = get_invmass_smeared(v1,v2);
0299
0300         //Fill Histogram
0301         h_InvMass->Fill( InvMass );
0302         h_InvMass_smeared->Fill( InvMass_smeared );
0303         // cout <<"\n"<<"InvMass  = "<< InvMass << "\t"  << endl;
0304
0305         //You have to calculate the rapitidy,eta, and transverse momentum for the daughter particles. Eta is pseudorapidity which tells you what "direction" or angle they div           ert at, rapidity is the natural unit for speed of the particle in relative spaces, and the transverse momentum tells you.... idrk tbh.
0306
0307         // calculate the rapidity
0308         rapt = 0.5 * log((E + pzt)/(E - pzt)); //rapidity of parent
0309
0310         // Rapidity of daughters
0311
0312         rap1 = 0.5 * log((E1 + pz1)/(E1 - pz1));
0313         rap2 = 0.5 * log((E2 + pz2)/(E2 - pz2));
0314         h_rap1->Fill( rap1 );
0315         h_rap2->Fill( rap2 );
0316         h_rapt->Fill( rapt );
0317
0318         // Pseudorapidity of daughters
0319         double theta1 = atan2( pt1, pz1 );
0320         eta1 = -log( tan( 0.5*theta1 ) );
0321
0322         double theta2 = atan2( pt2, pz2 );
0323         eta2 = -log( tan( 0.5*theta2 ) );
0324
0325         h_eta1->Fill( eta1 );
0326         h_eta2->Fill( eta2 );
0327
0328         double thetat = atan2( ptt, pzt );
0329         etat = -log( tan(0.5*thetat ) );
0330
0331         // Check that both daughters are in sPHENIX acceptance
0332         if ( fabs(eta1) < 1.1 && fabs(eta2) < 1.1 && pt1>0.3 && pt2>0.3 )
0333         {
0334             h_rapt_sphenix->Fill( rapt );
0335             h_InvMass_sphenix->Fill(InvMass);
0336             h_InvMass_smeared_sphenix->Fill(InvMass_smeared);
0337             h_pt_sphenix->Fill( ptt );
0338         }
0339
0340         nevents++;
0341
0342         tree->Fill();
0343     }
0344
0345     cout << "Processed " << nevents << endl;
0346
0347     // Scaling factor for integrated luminosity expected in sPHENIX Run
0348     Double_t scale = xsection*integ_lumi/nevents;
0349     cout << "Scale factor is " << endl;
0350     h_rapt->Sumw2();
0351     h_rapt_sphenix->Sumw2();
0352     h_pt->Sumw2();
0353     h_pt_sphenix->Sumw2();
0354     h_InvMass->Sumw2();
0355     h_InvMass_sphenix->Sumw2();
0356     h_InvMass_smeared->Sumw2();
0357     h_InvMass_smeared_sphenix->Sumw2();
0358     h_rapt->Scale( scale );
0359     h_rapt_sphenix->Scale( scale );
0360     h_pt->Scale( scale );
0361     h_pt_sphenix->Scale( scale );
0362     h_InvMass->Scale( scale );
0363     h_InvMass_sphenix->Scale( scale );
0364     h_InvMass_smeared->Scale( scale );
0365     h_InvMass_smeared_sphenix->Scale( scale );
0366
0367     cout << "Scale factor = " << scale << endl;
0368     cout << "Total events = " << h_rapt->Integral() << endl;
0369     cout << "Total events in sPHENIX acceptance = " << h_rapt_sphenix->Integral() << endl;
0370     cout << "Acceptance Factor = " << h_rapt_sphenix->Integral()/h_rapt->Integral() << endl;
0371
0372
0373     TCanvas *a = new TCanvas("a","Total Parent Momentum of JPsi",550,425);
0374     h_pt->SetLineColor(kBlue);
0375     h_pt->SetXTitle("p_{T} (GeV/c)");
0376     h_pt->SetYTitle("Counts");
0377     h_pt->Draw();
0378     a->SaveAs("PTofParent.png");
0379
0380     TCanvas *b = new TCanvas("b","Invariant Mass of JPsi",550,425);
0381     h_InvMass_smeared->SetLineColor(kBlack);
0382     h_InvMass_smeared->SetXTitle("mass (GeV/c^{2})");
0383     h_InvMass_smeared->SetYTitle("Counts");
0384     h_InvMass_smeared->Draw("ehist");
0385     h_InvMass_smeared_sphenix->SetLineColor(kBlue);
0386     h_InvMass_smeared_sphenix->SetMarkerColor(kBlue);
0387     h_InvMass_smeared_sphenix->SetXTitle("Mass (GeV/c^{2})");
0388     h_InvMass_smeared_sphenix->Draw("ehistsame");
0389     b->SaveAs("MassofUPCJPsi.png");
0390
0391     TCanvas *c = new TCanvas("c","Rapidity of E1 and E2",550,425);
0392     h_rap1->SetLineColor(kGreen);
0393     h_rap1->SetTitle("Rapidity e^{+},e^{-}");
0394     h_rap1->SetXTitle("Rapidity y");
0395     h_rap1->SetYTitle("Counts");
0396     h_rap1->Draw();
0397     h_rap2->SetLineColor(kRed);
0398     h_rap2->SetXTitle("Rapidity y");
0399     h_rap2->SetYTitle("Counts");
0400     h_rap2->Draw("same");
0401     c->SaveAs("RapofE1E2.png");
0402
0403     TCanvas *d = new TCanvas("d","Rapidity of UPC J/Psi",550,425);
0404     h_rapt->SetLineColor(kBlack);
0405     h_rapt->SetXTitle("Rapidity y");
0406     h_rapt->SetYTitle("Counts");
0407     h_rapt->Draw("ehist");
0408     h_rapt_sphenix->SetLineColor(kBlue);
0409     h_rapt_sphenix->SetMarkerColor(kBlue);
0410     h_rapt_sphenix->SetXTitle("Rapidity y");
0411     h_rapt_sphenix->Draw("ehistsame");
0412     d->SaveAs("RapofUPCJPsi.png");
0413
0414     TCanvas *e = new TCanvas("e","PZ",550,425);
0415     h_pz->SetLineColor(kBlue);
0416     h_pz->SetXTitle("p_{Z} GeV/c");
0417     h_pz->SetYTitle("Counts");
0418     h_pz->Draw();
0419     gPad->SetLogy(1);
0420     e->SaveAs("pz.png");
0421
0422     // Save histograms and TTree to file
0423     savefile->Write();
0424
0425 } //closing function