Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
track_kinematics.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file track_kinematics.C
1 void track_kinematics()
2 {
3  int triggersample = 30;
4  gStyle->SetOptStat(0);
5  TFile* infile = TFile::Open(Form("jet_pythia_%igev/trees/0.root",triggersample));
6 
7  //-----------------------------------------------------------------------------------------------------------------------
8  //Weight factors to enable splicing together mb, 10 GeV, and 30 GeV MC samples
9  //-----------------------------------------------------------------------------------------------------------------------
10  float csMB = 4.197*pow(10,9)*pow(10,-12); //I refuse to do math
11  float cs10GeV = 3.210*pow(10,6)*pow(10,-12);
12  float cs30GeV = 2.178*pow(10,3)*pow(10,-12);
13  float weight10GeV = cs10GeV/csMB;
14  float weight30GeV = cs30GeV/csMB;
15  float weight = weight30GeV;
16 
17  //---------------------------------------------------------------
18  //putting pi in a float for ease of use later on
19  //---------------------------------------------------------------
20  float pi = 3.1415926;
21 
22 
23  TTree* tree = (TTree*)infile->Get("tree");
24  //--------------------------------------------
25  // Vectors and floats for various
26  // quantities in the stored ttree
27  //--------------------------------------------
28  float zvtx;
29  std::vector<float> *trk_pt = {0};
30  std::vector<float> *trk_phi = {0};
31  std::vector<float> *trk_eta = {0};
32  std::vector<float> *trk_quality = {0};
33  std::vector<float> *trk_nhits = {0};
34 
35  std::vector<float> *tp_pt = {0};
36  std::vector<float> *tp_phi = {0};
37  std::vector<float> *tp_eta = {0};
38  std::vector<int> *tp_pid = {0};
39 
40  std::vector<float> *jet_pt = {0};
41  std::vector<float> *jet_phi = {0};
42  std::vector<float> *jet_eta = {0};
43 
44  std::vector<float> *tjet_pt = {0};
45  std::vector<float> *tjet_phi = {0};
46  std::vector<float> *tjet_eta = {0};
47  //------------------------------------------------------------------------
48  //Set the branch address for quantities in the ttree
49  //------------------------------------------------------------------------
50  tree->SetBranchAddress("trk_pt",&trk_pt);
51  tree->SetBranchAddress("trk_phi",&trk_phi);
52  tree->SetBranchAddress("trk_eta",&trk_eta);
53  tree->SetBranchAddress("trk_quality",&trk_quality);
54  tree->SetBranchAddress("nmvtx_hits",&trk_nhits);
55 
56  tree->SetBranchAddress("tjet_pt",&tjet_pt);
57  tree->SetBranchAddress("tjet_phi",&tjet_phi);
58  tree->SetBranchAddress("tjet_eta",&tjet_eta);
59 
60  tree->SetBranchAddress("tjet_pt",&jet_pt);
61  tree->SetBranchAddress("tjet_phi",&jet_phi);
62  tree->SetBranchAddress("tjet_eta",&jet_eta);
63 
64  tree->SetBranchAddress("tp_pt",&tp_pt);
65  tree->SetBranchAddress("tp_phi",&tp_phi);
66  tree->SetBranchAddress("tp_eta",&tp_eta);
67  tree->SetBranchAddress("tp_pid",&tp_pid);
68 
69  tree->SetBranchAddress("zvtx",&zvtx);
70 
71  //-----------------------------------------
72  // Define histograms for filling
73  //-----------------------------------------
74 
75  TH1D* h_deltar = new TH1D ("h_deltar",";delta R;",1000,0,0.4);
76  h_deltar->Sumw2();
77  TH2D* h_ptpart_pttrack = new TH2D("h_ptpart_pttrack",";pt_particle;pt_track/pt_particle",100,0,100,1000,0,2);
78  h_ptpart_pttrack->Sumw2();
79 
80  //------------------------------------------------
81  //Loop over the events in the ttree
82  //------------------------------------------------
83 
84  int nentries = tree->GetEntries();
85  for (int q = 0; q < nentries ;q++)
86  {
87  tree->GetEntry(q);
88  if (abs(zvtx) > 10){continue;} // reject events at large z vertex
89  int ntracks = trk_pt->size();
90  int n_particles = tp_pt->size();
91 
92  for (int j = 0;j < ntracks;j++)
93  {
94  float trkphi = trk_phi->at(j);
95  float trketa = trk_eta->at(j);
96  //-----------------------------------------
97  //Apply some track selections
98  // to select decent tracks
99  //-----------------------------------------
100  if (trk_nhits->at(j) < 4){continue;}
101  if (trk_quality->at(j) > 6){continue;}
102  float drmax = 0.4;
103  float partpt = 0;
104  for (int k = 0; k < n_particles;k++)
105  {
106  float deltaphi = abs(trkphi - tp_phi->at(k));
107  if (deltaphi > pi) {deltaphi = 2*pi-deltaphi;}
108  float deltaeta = abs(trketa-tp_eta->at(k));
109  float deltar = TMath::Sqrt(deltaphi*deltaphi + deltaeta*deltaeta);
110  if (deltar < drmax)
111  {
112  drmax = deltar;
113  partpt = tp_pt->at(k);
114  }
115  }
116  h_deltar->Fill(drmax,weight);
117  if (drmax < 0.02)
118  {
119  h_ptpart_pttrack->Fill(partpt,trk_pt->at(j)/partpt,weight);
120  }
121  }
122  }
123 
124 
125  TCanvas*c1 = new TCanvas("c1","",700,500);
126  h_ptpart_pttrack->Draw("COLZ");
127 
128 
129  TCanvas*c2 = new TCanvas("c2","",700,500);
130  h_deltar->Draw();
131 
132 
133 
134 
135 }