JetRTrack/macros/track_kinematics.C

0001 void track_kinematics()
0002 {
0003   int triggersample = 30;
0004   gStyle->SetOptStat(0);
0005   TFile* infile = TFile::Open(Form("jet_pythia_%igev/trees/0.root",triggersample));
0006
0007   //-----------------------------------------------------------------------------------------------------------------------
0008   //Weight factors to enable splicing together mb, 10 GeV, and 30 GeV MC samples
0009   //-----------------------------------------------------------------------------------------------------------------------
0010   float csMB = 4.197*pow(10,9)*pow(10,-12); //I refuse to do math
0011   float cs10GeV = 3.210*pow(10,6)*pow(10,-12);
0012   float cs30GeV = 2.178*pow(10,3)*pow(10,-12);
0013   float weight10GeV = cs10GeV/csMB;
0014   float weight30GeV = cs30GeV/csMB;
0015   float weight = weight30GeV;
0016
0017   //---------------------------------------------------------------
0018   //putting pi in a float for ease of use later on
0019   //---------------------------------------------------------------
0020   float pi = 3.1415926;
0021
0022
0023   TTree* tree = (TTree*)infile->Get("tree");
0024   //--------------------------------------------
0025   // Vectors and floats for various
0026   // quantities in the stored ttree
0027   //--------------------------------------------
0028   float zvtx;
0029   std::vector<float> *trk_pt = {0};
0030   std::vector<float> *trk_phi = {0};
0031   std::vector<float> *trk_eta = {0};
0032   std::vector<float> *trk_quality = {0};
0033   std::vector<float> *trk_nhits = {0};
0034
0035   std::vector<float> *tp_pt = {0};
0036   std::vector<float> *tp_phi = {0};
0037   std::vector<float> *tp_eta = {0};
0038   std::vector<int> *tp_pid = {0};
0039
0040   std::vector<float> *jet_pt = {0};
0041   std::vector<float> *jet_phi = {0};
0042   std::vector<float> *jet_eta = {0};
0043
0044   std::vector<float> *tjet_pt = {0};
0045   std::vector<float> *tjet_phi = {0};
0046   std::vector<float> *tjet_eta = {0};
0047   //------------------------------------------------------------------------
0048   //Set the branch address for quantities in the ttree
0049   //------------------------------------------------------------------------
0050   tree->SetBranchAddress("trk_pt",&trk_pt);
0051   tree->SetBranchAddress("trk_phi",&trk_phi);
0052   tree->SetBranchAddress("trk_eta",&trk_eta);
0053   tree->SetBranchAddress("trk_quality",&trk_quality);
0054   tree->SetBranchAddress("nmvtx_hits",&trk_nhits);
0055
0056   tree->SetBranchAddress("tjet_pt",&tjet_pt);
0057   tree->SetBranchAddress("tjet_phi",&tjet_phi);
0058   tree->SetBranchAddress("tjet_eta",&tjet_eta);
0059
0060   tree->SetBranchAddress("tjet_pt",&jet_pt);
0061   tree->SetBranchAddress("tjet_phi",&jet_phi);
0062   tree->SetBranchAddress("tjet_eta",&jet_eta);
0063
0064   tree->SetBranchAddress("tp_pt",&tp_pt);
0065   tree->SetBranchAddress("tp_phi",&tp_phi);
0066   tree->SetBranchAddress("tp_eta",&tp_eta);
0067   tree->SetBranchAddress("tp_pid",&tp_pid);
0068
0069   tree->SetBranchAddress("zvtx",&zvtx);
0070
0071   //-----------------------------------------
0072   // Define histograms for filling
0073   //-----------------------------------------
0074
0075   TH1D* h_deltar = new TH1D ("h_deltar",";delta R;",1000,0,0.4);
0076   h_deltar->Sumw2();
0077   TH2D* h_ptpart_pttrack = new TH2D("h_ptpart_pttrack",";pt_particle;pt_track/pt_particle",100,0,100,1000,0,2);
0078   h_ptpart_pttrack->Sumw2();
0079
0080   //------------------------------------------------
0081   //Loop over the events in the ttree
0082   //------------------------------------------------
0083
0084   int nentries = tree->GetEntries();
0085   for (int q = 0; q < nentries ;q++)
0086     {
0087       tree->GetEntry(q);
0088       if (abs(zvtx) > 10){continue;} // reject events at large z vertex
0089       int ntracks = trk_pt->size();
0090       int n_particles = tp_pt->size();
0091
0092       for (int j = 0;j < ntracks;j++)
0093         {
0094           float trkphi = trk_phi->at(j);
0095           float trketa = trk_eta->at(j);
0096       //-----------------------------------------
0097       //Apply some track selections
0098       // to select decent tracks
0099       //-----------------------------------------
0100       if (trk_nhits->at(j) < 4){continue;}
0101       if (trk_quality->at(j) > 6){continue;}
0102       float drmax = 0.4;
0103       float partpt = 0;
0104       for (int k = 0; k < n_particles;k++)
0105         {
0106           float deltaphi = abs(trkphi - tp_phi->at(k));
0107           if (deltaphi > pi) {deltaphi = 2*pi-deltaphi;}
0108           float deltaeta = abs(trketa-tp_eta->at(k));
0109           float deltar = TMath::Sqrt(deltaphi*deltaphi + deltaeta*deltaeta);
0110           if (deltar < drmax)
0111         {
0112           drmax = deltar;
0113           partpt = tp_pt->at(k);
0114         }
0115         }
0116       h_deltar->Fill(drmax,weight);
0117       if (drmax < 0.02)
0118         {
0119           h_ptpart_pttrack->Fill(partpt,trk_pt->at(j)/partpt,weight);
0120         }
0121     }
0122     }
0123
0124
0125   TCanvas*c1 = new TCanvas("c1","",700,500);
0126   h_ptpart_pttrack->Draw("COLZ");
0127
0128
0129   TCanvas*c2 = new TCanvas("c2","",700,500);
0130   h_deltar->Draw();
0131
0132
0133
0134
0135 }