Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
event_topology_reco.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file event_topology_reco.C
1 #include "tau_commons.h"
2 
10 double Average(vector<double> v)
11 {
12  double sum=0;
13  for(int i=0;(unsigned)i<v.size();i++)
14  sum+=v[i];
15  return sum/v.size();
16 }
17 
18 
19 int event_topology_reco()
20 {
21 
22  gStyle->SetOptStat(0);
23 
24  unsigned col1 = kOrange+7;
25  unsigned col2 = kBlue+2;
26  unsigned col3 = kGreen+2;
27 
28  // labels for geant files //
29  string type[3] = {"3pion","DIScharged","DISneutral"};
30  string seed[10] = {"1","2","3","4","5","6","7","8","9","10"};
31  string R_max = "5";
32  string nevents;
33  string file;
34 
35  /* open inout files and merge trees */
36  TChain chain_CC("event");
37  TChain chain_NC("event");
38  TChain chain_tau("event");
39 
40 
41  for (int b = 0; b<3; b++){
42  for (int g = 0; g<10; g++){
43  // Skip file that is used to test identification //
44  if(b==0 & g==0) continue;
45  // Set some labels //
46  if(b==0) nevents = "1000";
47  if(b==1 || b==2) nevents = "100";
48 
49  if(b==0) file = "/gpfs/mnt/gpfs02/phenix/scratch/spjeffas/data/LeptoAna_p250_e20_"+nevents+"events_"+seed[g]+"seed_"+type[b]+"_r0"+R_max+".root";
50  if(b==1 || b==2) file = "/gpfs/mnt/gpfs02/phenix/scratch/spjeffas/data/LeptoAna_p250_e20_"+nevents+"events_"+seed[g]+"seed_"+type[b]+"_tau_r0"+R_max+".root";
51 
52  if(b==0) chain_tau.Add(file.c_str());
53  if(b==1) chain_CC.Add(file.c_str());
54  if(b==2) chain_NC.Add(file.c_str());
55  }
56  }
57 
58  // Add file that will be used to test //
59  string inputFile = "/gpfs/mnt/gpfs02/phenix/scratch/spjeffas/data/LeptoAna_p250_e20_1000events_1seed_"+type[0]+"_r0"+R_max+".root";
60 
61  TFile *f = TFile::Open(inputFile.c_str());
62  TTree *t = (TTree*)f->Get("event");
63 
64  // Create temporary canvas
65  TCanvas* ctemp = new TCanvas();
66 
67  vector< TString > observables;
68  vector< TString > observables_name;
69 
70  vector< float > plots_xmin;
71  vector< float > plots_xmax;
72 
73  // Missing energy
74  observables.push_back( "Et_miss" );
75  observables_name.push_back( "p_{T}^{miss} (GeV)" );
76  plots_xmin.push_back(0);
77  plots_xmax.push_back(50);
78 
79  // Acoplanarity
80  observables.push_back( "Et_miss_phi" );
81  observables_name.push_back( "#Delta#phi_{miss - #tau jet}" );
82  plots_xmin.push_back(0);
83  plots_xmax.push_back(180);
84 
85  TString name_CC_base("h_CC_");
86  TString name_NC_base("h_NC_");
87  TString name_tau_base("h_tau_");
88 
89  const int nvar = observables.size();
90 
91  TH1F* h_CC[nvar];
92  TH1F* h_NC[nvar];
93  TH1F* h_tau[nvar];
94 
95  TString name_CC_i = name_CC_base;
96  TString name_NC_i = name_NC_base;
97  TString name_tau_i = name_tau_base;
98 
99  // create histograms
100  for ( int l = 0; l < observables.size(); l++ ){
101  name_CC_i = name_CC_base;
102  name_NC_i = name_NC_base;
103  name_tau_i = name_tau_base;
104 
105  name_CC_i.Append(l);
106  name_NC_i.Append(l);
107  name_tau_i.Append(l);
108 
109  h_CC[l] = new TH1F( name_CC_i, "", 50, plots_xmin.at(l), plots_xmax.at(l));
110  h_NC[l] = new TH1F( name_NC_i, "", 50, plots_xmin.at(l), plots_xmax.at(l));
111  h_tau[l] = new TH1F( name_tau_i, "", 50, plots_xmin.at(l), plots_xmax.at(l));
112  }
113 
114  // Fill Histograms //
115  for(int i = 0;i < observables.size(); i++){
116 
117  ctemp->cd();
118 
119  name_CC_i = name_CC_base;
120  name_NC_i = name_NC_base;
121  name_tau_i = name_tau_base;
122 
123  name_CC_i.Append(i);
124  name_NC_i.Append(i);
125  name_tau_i.Append(i);
126 
127  // Fill Et_miss straight from tree //
128  if(i==0){
129  // Fill CC histogram //
130  int n_CC = chain_CC.Draw( observables.at(i) + " >> " + name_CC_i, "", "goff");
131  h_CC[i]->Scale(1./h_CC[i]->Integral());
132  h_CC[i]->GetXaxis()->SetTitle( observables_name.at(i) );
133  h_CC[i]->GetYaxis()->SetRangeUser(0, .13 );
134  h_CC[i]->SetLineColor(col1);
135  h_CC[i]->SetFillColor(col1);
136  h_CC[i]->SetFillStyle(0);
137  h_CC[i]->GetYaxis()->SetTitle("# entries / #Sigma entries");
138  h_CC[i]->GetYaxis()->SetTitleOffset(1.5);
139 
140  // Fill NC histogram //
141  int n_NC = chain_NC.Draw( observables.at(i) + " >> " + name_NC_i, "", "goff");
142  h_NC[i]->Scale(1./h_NC[i]->Integral());
143  h_NC[i]->SetLineColor(col3);
144  h_NC[i]->SetFillColor(col3);
145  h_NC[i]->SetFillStyle(0);
146 
147  // Fill LQ histogram //
148  int n_tau = chain_tau.Draw( observables.at(i) + " >> " + name_tau_i, "", "goff");
149  h_tau[i]->Scale(1./h_tau[i]->Integral());
150  h_tau[i]->SetLineColor(col2);
151  h_tau[i]->SetFillColor(col2);
152  h_tau[i]->SetFillStyle(0);
153 
154  TCanvas *c1 = new TCanvas();
155 
156  h_CC[i]->DrawClone("");
157  h_NC[i]->DrawClone("sames");
158  h_tau[i]->DrawClone("sames");
159  gPad->RedrawAxis();
160 
161  TLegend *legend = new TLegend(0.4,0.6,0.7,0.89);
162  legend->AddEntry(h_CC[i],"Neutral Current","l");
163  legend->AddEntry(h_NC[i],"Charged Current","l");
164  legend->AddEntry(h_tau[i],"Leptoquark","l");
165  legend->Draw();
166  }
167 
168  // Calculate acoplanarity and then fill histogram //
169  if(i==1){
170 
171  float Et_miss_phi_tau, Et_miss_phi_CC, Et_miss_phi_NC;
172  float tau_phi_tau, tau_phi_CC, tau_phi_NC;
173 
174  // Read tau phi and missing momentum phi from tree //
175  chain_CC.SetBranchAddress("Et_miss_phi",&Et_miss_phi_CC);
176  chain_CC.SetBranchAddress("reco_tau_phi",&tau_phi_CC);
177 
178  chain_NC.SetBranchAddress("Et_miss_phi",&Et_miss_phi_NC);
179  chain_NC.SetBranchAddress("reco_tau_phi",&tau_phi_NC);
180 
181  chain_tau.SetBranchAddress("Et_miss_phi",&Et_miss_phi_tau);
182  chain_tau.SetBranchAddress("reco_tau_phi",&tau_phi_tau);
183 
184  // Phi angle between reco tau and miss mom //
185  double delta_phi;
186 
187  // Calculate delta phi for each event type //
188  for(int j = 0; j<n_CC;j++){
189  chain_CC.GetEntry(j);
190  delta_phi = fabs(Et_miss_phi_CC - tau_phi_CC) * 180 / TMath::Pi();
191  if(delta_phi > 180) delta_phi = 360 - delta_phi;
192  h_CC[i]->Fill(delta_phi);
193  }
194 
195  for(int l = 0; l<n_NC;l++){
196  chain_NC.GetEntry(l);
197  delta_phi = fabs(Et_miss_phi_NC - tau_phi_NC) * 180 / TMath::Pi();
198  if(delta_phi > 180) delta_phi = 360 - delta_phi;
199  h_NC[i]->Fill(delta_phi);
200  }
201 
202  for(int k = 0; k<n_tau;k++){
203  chain_tau.GetEntry(k);
204  delta_phi = fabs(Et_miss_phi_tau - tau_phi_tau) * 180 / TMath::Pi();
205  if(delta_phi > 180) delta_phi = 360 - delta_phi;
206  h_tau[i]->Fill(delta_phi);
207  }
208 
209  //Fill histograms //
210  h_CC[i]->Scale(1./h_CC[i]->Integral());
211  h_CC[i]->GetXaxis()->SetTitle( observables_name.at(i) );
212  h_CC[i]->GetYaxis()->SetRangeUser(0, 0.27 );
213  h_CC[i]->SetLineColor(col1);
214  h_CC[i]->SetFillColor(col1);
215  h_CC[i]->SetFillStyle(0);
216  h_CC[i]->GetYaxis()->SetTitle("# entries / #Sigma entries");
217  h_CC[i]->GetYaxis()->SetTitleOffset(1.5);
218 
219  h_NC[i]->Scale(1./h_NC[i]->Integral());
220  h_NC[i]->SetLineColor(col3);
221  h_NC[i]->SetFillColor(col3);
222  h_NC[i]->SetFillStyle(0);
223 
224 
225  h_tau[i]->Scale(1./h_tau[i]->Integral());
226  h_tau[i]->SetLineColor(col2);
227  h_tau[i]->SetFillColor(col2);
228  h_tau[i]->SetFillStyle(0);
229 
230  TCanvas *c1 = new TCanvas();
231 
232  h_CC[i]->DrawClone("");
233  h_NC[i]->DrawClone("sames");
234  h_tau[i]->DrawClone("sames");
235  gPad->RedrawAxis();
236 
237  TLegend *legend = new TLegend(0.4,0.6,0.7,0.89);
238  legend->AddEntry(h_NC[i],"Neutral Current","l");
239  legend->AddEntry(h_CC[i],"Charged Current","l");
240  legend->AddEntry(h_tau[i],"Leptoquark","l");
241  legend->Draw();
242 
243  }
244 
245  }
246 
247  return 0;
248 }