readNtuple.C
Go to the documentation of this file.
1 #define readNtuple_cxx
2 #include "readNtuple.h"
3 #include <TH2.h>
4 #include <TStyle.h>
5 #include <TCanvas.h>
6 #include "TLorentzVector.h"
7 #include <iostream>
8 #include "TH3F.h"
9 #include "TProfile.h"
10 #include "TProfile2D.h"
11 #include "TTree.h"
12 #include "string.h"
13 #include "TString.h"
14 
15 using namespace std;
16 using namespace TMVA;
17 
18 void readNtuple::Loop(const char* TypeNameNu, const char* TypeNameInt)
19 {
20 
21  bool ccnumu=false;
22  bool ccnumubar=false;
23  bool ccnue=false;
24  bool ccnuebar=false;
25  bool nc=false;
26 
27  TString nutype = (TString)TypeNameNu;
28  if( nutype=="ccnumu" ){
29  ccnumu=true;
30  cout<<"Analyze the CC numu, ";
31  }
32  else if( nutype=="ccnumubar" ){
33  ccnumubar=true;
34  cout<<"Analyze the CC numubar, ";
35  }
36  else if( nutype=="ccnue" ){
37  ccnue=true;
38  cout<<"Analyze the CC nue, ";
39  }
40  else if( nutype=="ccnuebar") {
41  ccnuebar=true;
42  cout<<"Analyze the CC nuebar, ";
43  }
44  else if( nutype=="nc" ){
45  nc=true;
46  cout<<"Analyze the NC, ";
47  }
48  else{
49  cout<<"Please enter the right type of nu interactions: ccnumu, ccnumubar, ccnue, ccnuebar, nc"<<endl;
50  exit(0);
51  }
52 
53  bool inclusive=false;
54  bool qeonly=false;
55  bool disonly=false;
56  bool resonly=false;
57  bool meconly=false;
58  bool cohonly=false;
59 
60  TString protype=(TString)TypeNameInt;
61  if( protype=="all" ){
62  inclusive=true;
63  cout<<"inclusive events."<<endl;
64  }
65  else if( protype=="qe"){
66  qeonly=true;
67  cout<<"QE only events."<<endl;
68  }
69  else if( protype="res"){
70  resonly=true;
71  cout<<"RES only events."<<endl;
72  }
73  else if( protype=="dis"){
74  disonly=true;
75  cout<<"DIS only events."<<endl;
76  }
77  else if( protype="mec"){
78  meconly=true;
79  cout<<"MEC only events."<<endl;
80  }
81  else if( protype="coh"){
82  cohonly=true;
83  cout<<"COH only events."<<endl;
84  }
85  else{
86  cout<<"please enerty the right interaction process: all, qe, res, dis, mec, coh"<<endl;
87  exit(0);
88  }
89 
90  if (fChain == 0) return;
91 
92  Int_t nentries = fChain->GetEntries();
93  cout<<"Total events are "<<nentries<<endl;
94 
95  //read-in the flux spectrum
96  TFile *flux=new TFile("doc/nu_flux_v08.root");
97  TString fluxhist="";
98  if( ccnumu ) fluxhist="numu_flux";
99  else if( ccnumubar ) fluxhist="numubar_flux";
100  else if( ccnue ) fluxhist="nue_flux";
101  else if( ccnuebar ) fluxhist="nuebar_flux";
102  else fluxhist="numu_flux";
103 
104  TH1F *hflux=(TH1F*)flux->Get(fluxhist);//nu flux in unit of per E10 POT per m2
105  double integrate_flux=hflux->Integral();
106 
107  //# of nucleons in fully active detector
108  //the Z is cut off at 12.75m
109  double Nnucleons=11.5279;//unit E31
110 
111  //vertex cuts
112  double minVz=0.;//cm
113  double maxVz=1275.0;//cm
114  double minVy=-150.0;//cm
115  double maxVy=150.0; //cm
116  double minVx=-150.0;//cm
117  double maxVx=150.0;//cm
118 
119  double corr_Nnucleons= (maxVz-minVz)*(maxVy-minVy)*(maxVx-minVx)/(1275.0*400.0*400.0);
120  Nnucleons = Nnucleons*corr_Nnucleons;
121 
122  cout<<"# of nucleons are "<<Nnucleons<<endl;
123 
124  //pot of MC
125  double mcpot=0.00001;//E21
126 
127  int Nevts=0;
128  int Nfinal=0;
129 
130  for( int ii=0; ii<nentries; ii++)
131  //for( int ii=0; ii<10000; ii++)
132  {
133  fChain->GetEntry(ii);
134  ++Nevts;
135  if( Nevts%100000==0 ) cout<<Nevts<<" processed ..."<<endl;
136 
137  if( ccnumu ){
138  if( nu_ccnc != 0 ) continue;
139  if( nu_PDG != 14 ) continue;
140  if( nu_origPDG != 14 ) continue;
141  }
142  else if( ccnumubar ){
143  if( nu_ccnc != 0 ) continue;
144  if( nu_PDG != -14 ) continue;
145  if( nu_origPDG != -14 ) continue;
146  }
147  else if( ccnue ){
148  if( nu_ccnc != 0 ) continue;
149  if( nu_PDG != 12 ) continue;
150  if( nu_origPDG != 12 ) continue;
151  }
152  else if( ccnuebar ){
153  if( nu_ccnc != 0 ) continue;
154  if( nu_PDG != -12 ) continue;
155  if( nu_origPDG != -12 ) continue;
156  }
157  else if( nc ){
158  if( nu_ccnc != 1 ) continue;
159  if( nu_PDG != nu_origPDG ) continue;
160  }
161 
162  if( qeonly ){
163  if( nu_mode != 0 ) continue;
164  }
165  else if( resonly ){
166  if( nu_mode != 1 ) continue;
167  }
168  else if( disonly ){
169  if( nu_mode != 2 ) continue;
170  }
171  else if( meconly ){
172  if( nu_mode != 10 ) continue;
173  }
174  else if( cohonly ){
175  if( nu_mode != 3 ) continue;
176  }
177 
178  if( nu_Vx>maxVx || nu_Vx<minVx ) continue;
179  if( nu_Vy>maxVy || nu_Vy<minVy ) continue;
180  if( nu_Vz>maxVz || nu_Vz<minVz ) continue;
181 
182  int Xbin=hflux->FindBin(nu_trueE);
183  double flux_Ebin=0.0;
184  if( hflux->GetBinContent(Xbin)>0. ) flux_Ebin=1.0/hflux->GetBinContent(Xbin);
185  double total_weight=flux_Ebin/(mcpot*Nnucleons);
186  double diff_weight=1.0/(integrate_flux*mcpot*Nnucleons);
187 
188  myhists->Xsec_Enu->Fill(nu_trueE,total_weight);
189  myhists->diffXsec_lepE->Fill(lep_E,diff_weight);
190  TLorentzVector lepton;
191  lepton.SetPxPyPzE(lep_Px,lep_Py,lep_Pz,lep_E);
192  myhists->diffXsec_lepcostheta->Fill(cos(lepton.Eta()),diff_weight);
193 
194  ++Nfinal;
195  }//loop all events
196 
197  //dividing the bin width for diff. Xsec
198  for( int ibin=1; ibin<myhists->diffXsec_lepE->GetNbinsX()+1; ++ibin ){
199  double binwidth=myhists->diffXsec_lepE->GetBinWidth(ibin);
200  double bincon=myhists->diffXsec_lepE->GetBinContent(ibin)/binwidth;
201  double binerr=0.;
202  if( bincon>0.0 ) binerr=(myhists->diffXsec_lepE->GetBinError(ibin)/myhists->diffXsec_lepE->GetBinContent(ibin))*bincon;
203  myhists->diffXsec_lepE->SetBinContent(ibin,bincon);
204  myhists->diffXsec_lepE->SetBinError(ibin,binerr);
205  }
206 
207  cout<<"# selected events are "<<Nfinal<<endl;
208 
209 }
210 
Loaders::FluxType flux
unsigned int nutype(unsigned int u)
Definition: runWimpSim.h:122
Long64_t nentries
TH1F * hflux
Definition: Xdiff_gwt.C:121
TChain * fChain
OStream cout
Definition: OStream.cxx:6
exit(0)
enum BeamMode nc
T cos(T number)
Definition: d0nt_math.hpp:78
Definition: tmvaglob.h:28
virtual void Loop(const char *TypeNameNu, const char *TypeNameInt)
Definition: readNtuple.C:18