hbcup/macro/ana_hijbkg.C

0001 // This macro reads and analyzes the TTree created in hijbkg_upc.cc
0002 //
0003
0004 #include <cmath>
0005 #include <TFile.h>
0006 #include <TTree.h>
0007 #include <TMath.h>        // Include the TMath header for trigonometric functions
0008 #include <TH1F.h>
0009 #include <TCanvas.h>
0010 #include <TParticle.h>
0011 #include <TDatabasePDG.h>
0012
0013 void ana_hijbkg(const char* filename = "hijbkg.root")
0014 {
0015   TDatabasePDG *pdg = new TDatabasePDG;
0016
0017   TFile *tfile= new TFile(filename,"READ");
0018   TTree* tree = (TTree*)tfile->Get("T"); // the TTree in the HIJING file
0019
0020   vector<float> *m_pt = 0;
0021   Int_t m_evt = 0;
0022   Float_t m_b = 0;
0023   Float_t m_cent = 0;
0024   TParticle *m_part[2] = {0};
0025
0026   //set the branches of the tree
0027   tree->SetBranchAddress("evt", &m_evt);
0028   tree->SetBranchAddress("b", &m_b);
0029   tree->SetBranchAddress("part0", &m_part[0]);
0030   tree->SetBranchAddress("part1", &m_part[1]);
0031
0032   TFile *savefile = new TFile("hijbkg_results.root","RECREATE");
0033
0034   // Plot histograms for pt, eta, phi, mass
0035   TH1F *h_pt = new TH1F("h_pt", "Pt Distribution", 200, 0, 2.0 );
0036   TH1F *h_eta = new TH1F("h_eta", "Eta Distribution", 200, -6, 6 );
0037   TH1F *h_rap = new TH1F("h_rap", "Rapidity Distribution", 200, -4, 4 );
0038   TH1F *h_phi = new TH1F("h_phi", "Phi Distribution", 100, -M_PI, M_PI );
0039   TH1F *h_dphi = new TH1F("h_dphi", "dPhi Distribution", 100, -M_PI, M_PI );
0040   TH1F *h_mass = new TH1F("h_mass", "Inv Mass", 200, 0., 6.);
0041   h_pt->SetLineColor(4);
0042   h_eta->SetLineColor(4);
0043   h_phi->SetLineColor(4);
0044   h_mass->SetLineColor(4);
0045
0046   // Assuming all tracks are electrons
0047   TH1F *he_pt = new TH1F("he_pt", "Pt, assuming e+/e-", 200, 0, 2.0 );
0048   TH1F *he_eta = new TH1F("he_eta", "Eta, assuming e+/e- ", 200, -6, 6 );
0049   TH1F *he_rap = new TH1F("he_rap", "Rapidity, assuming e+/e-", 200, -4, 4 );
0050   TH1F *he_phi = new TH1F("he_phi", "Phi, assuming e+/e- ", 100, -M_PI, M_PI );
0051   TH1F *he_dphi = new TH1F("he_dphi", "dPhi, assuming e+/e- ", 100, -M_PI, M_PI );
0052   TH1F *he_mass = new TH1F("he_mass", "Mass, assuming e+/e-", 200, 0., 6.);
0053
0054   // With cut on J/Psi mass
0055   TH1F *hej_pt = new TH1F("hej_pt", "Pt, J/Psi mass cut", 200, 0, 2.0 );
0056   TH1F *hej_eta = new TH1F("hej_eta", "Eta, J/Psi mass cut ", 200, -6, 6 );
0057   TH1F *hej_rap = new TH1F("hej_rap", "Rapidity, J/Psi mass cut", 200, -4, 4 );
0058   TH1F *hej_dphi = new TH1F("hej_dphi", "dPhi, J/Psi mass cut", 100, -M_PI, M_PI );
0059   hej_pt->SetLineColor(2);
0060   hej_eta->SetLineColor(2);
0061   hej_rap->SetLineColor(2);
0062   hej_dphi->SetLineColor(2);
0063   hej_pt->SetXTitle("#p_T (GeV/c)");
0064
0065   TLorentzVector part_v4[2];
0066   TLorentzVector sum_v4;
0067
0068   TLorentzVector epart_v4[2]; // assuming these are e+/e- pairs
0069   TLorentzVector esum_v4;
0070
0071   // Loop over the tree
0072   Int_t nEntries = tree->GetEntries();
0073   for (int iEntry = 0; iEntry < nEntries; iEntry++)
0074   {
0075     tree->GetEntry(iEntry);
0076
0077     // Print out a message every 1000 events
0078     if ( iEntry%10000 == 0 ) cout << iEntry << "\t" << m_evt << endl;
0079
0080     //cout << m_part[0]->GetPdgCode() << "\t" << m_part[0]->Px() << endl;
0081
0082     Double_t charge[2];
0083     for (int ipart=0; ipart<2; ipart++)
0084     {
0085       charge[ipart] = m_part[ipart]->GetPDG()->Charge();
0086
0087       if ( fabs( m_part[ipart]->GetPDG()->PdgCode() ) == 11 )
0088       {
0089         cout << m_part[ipart]->GetPDG()->PdgCode() << "\t" << charge[ipart] << endl;
0090       }
0091
0092       double px = m_part[ipart]->Px();
0093       double py = m_part[ipart]->Py();
0094       double pz = m_part[ipart]->Pz();
0095
0096       part_v4[ipart].SetPxPyPzE( px, py, pz, m_part[ipart]->Energy() );
0097
0098       // now assume this is an electron
0099       const double m_e = 0.511e-3;    // e- mass in GeV
0100       double e_energy = sqrt( px*px + py*py + pz*pz + m_e*m_e );
0101       epart_v4[ipart].SetPxPyPzE( px, py, pz, e_energy );
0102
0103
0104       // check for bad eta
0105       if ( fabs(part_v4[ipart].Eta()) > 1.1 )
0106       {
0107         cout << "ERROR " << iEntry << m_part[ipart]->GetPdgCode() << endl;
0108       }
0109     }
0110
0111     // skip if they are like sign pairs
0112     if ( charge[0]*charge[1] > 0 )
0113     {
0114       continue;
0115     }
0116
0117     sum_v4 = part_v4[0] + part_v4[1];
0118     esum_v4 = epart_v4[0] + epart_v4[1];
0119
0120     // Fill the histograms with the data from the TTree
0121     h_pt->Fill( sum_v4.Pt()  );
0122     h_eta->Fill( sum_v4.Eta() );
0123     h_phi->Fill( sum_v4.Phi() );
0124     h_mass->Fill( sum_v4.M() );
0125     h_rap->Fill( sum_v4.Rapidity() );
0126     h_dphi->Fill( part_v4[0].DeltaPhi( part_v4[1] ) );
0127
0128     // electron PID assumption
0129     he_pt->Fill( esum_v4.Pt()  );
0130     he_eta->Fill( esum_v4.Eta() );
0131     he_phi->Fill( esum_v4.Phi() );
0132     he_mass->Fill( esum_v4.M() );
0133     he_rap->Fill( esum_v4.Rapidity() );
0134     he_dphi->Fill( epart_v4[0].DeltaPhi( epart_v4[1] ) );
0135
0136     // in J/Psi mass window
0137     if ( fabs( esum_v4.M() - 3.1 ) < 0.3  )
0138     {
0139       hej_pt->Fill( esum_v4.Pt()  );
0140       hej_eta->Fill( esum_v4.Eta() );
0141       hej_rap->Fill( esum_v4.Rapidity() );
0142       hej_dphi->Fill( epart_v4[0].DeltaPhi( epart_v4[1] ) );
0143     }
0144   }
0145
0146   // Draw and save the histograms
0147   TCanvas *ac = new TCanvas("ac","pt",550,425);
0148   h_pt->SetMinimum(0.5);
0149   h_pt->Draw();
0150   he_pt->Draw("same");
0151   hej_pt->Draw("same");
0152   gPad->SetLogy(1);
0153   ac->SaveAs("h_pt.png");
0154
0155   TCanvas *dc = new TCanvas("dc","eta",550,425);
0156   h_eta->SetMinimum(0.5);
0157   h_eta->Draw();
0158   he_eta->Draw("same");
0159   hej_eta->Draw("same");
0160   dc->SaveAs("h_eta.png");
0161
0162   TCanvas *fc = new TCanvas("fc","rapidity",550,425);
0163   h_rap->SetMinimum(0.5);
0164   h_rap->Draw();
0165   he_rap->Draw("same");
0166   hej_rap->Draw("same");
0167   dc->SaveAs("h_eta.png");
0168
0169   TCanvas *bc = new TCanvas("bc","phi",550,425);
0170   h_phi->SetMinimum(0.5);
0171   h_phi->Draw();
0172   he_phi->Draw("same");
0173   bc->SaveAs("h_phi.png");
0174
0175   TCanvas *gc = new TCanvas("gc","dphi",550,425);
0176   h_dphi->SetMinimum(0.5);
0177   h_dphi->Draw();
0178   he_dphi->Draw("same");
0179   hej_dphi->Draw("same");
0180   dc->SaveAs("h_eta.png");
0181
0182   TCanvas *ec = new TCanvas("ec","mass",550,425);
0183   h_mass->SetMinimum(0.5);
0184   h_mass->Draw();
0185   he_mass->Draw("same");
0186   gPad->SetLogy(1);
0187   ec->SaveAs("h_mass.png");
0188
0189   savefile->Write();
0190 }