6 #include "TLorentzVector.h" 10 #include "TProfile2D.h" 27 TString
nutype = (TString)TypeNameNu;
28 if( nutype==
"ccnumu" ){
30 cout<<
"Analyze the CC numu, ";
32 else if( nutype==
"ccnumubar" ){
34 cout<<
"Analyze the CC numubar, ";
36 else if( nutype==
"ccnue" ){
38 cout<<
"Analyze the CC nue, ";
40 else if( nutype==
"ccnuebar") {
42 cout<<
"Analyze the CC nuebar, ";
44 else if( nutype==
"nc" ){
46 cout<<
"Analyze the NC, ";
49 cout<<
"Please enter the right type of nu interactions: ccnumu, ccnumubar, ccnue, ccnuebar, nc"<<
endl;
60 TString protype=(TString)TypeNameInt;
65 else if( protype==
"qe"){
69 else if( protype=
"res"){
73 else if( protype==
"dis"){
77 else if( protype=
"mec"){
81 else if( protype=
"coh"){
86 cout<<
"please enerty the right interaction process: all, qe, res, dis, mec, coh"<<
endl;
93 cout<<
"Total events are "<<nentries<<
endl;
96 TFile *
flux=
new TFile(
"doc/nu_flux_v08.root");
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";
104 TH1F *
hflux=(TH1F*)flux->Get(fluxhist);
105 double integrate_flux=hflux->Integral();
109 double Nnucleons=11.5279;
119 double corr_Nnucleons= (maxVz-minVz)*(maxVy-minVy)*(maxVx-minVx)/(1275.0*400.0*400.0);
120 Nnucleons = Nnucleons*corr_Nnucleons;
122 cout<<
"# of nucleons are "<<Nnucleons<<
endl;
125 double mcpot=0.00001;
135 if( Nevts%100000==0 )
cout<<Nevts<<
" processed ..."<<
endl;
138 if( nu_ccnc != 0 )
continue;
139 if( nu_PDG != 14 )
continue;
140 if( nu_origPDG != 14 )
continue;
142 else if( ccnumubar ){
143 if( nu_ccnc != 0 )
continue;
144 if( nu_PDG != -14 )
continue;
145 if( nu_origPDG != -14 )
continue;
148 if( nu_ccnc != 0 )
continue;
149 if( nu_PDG != 12 )
continue;
150 if( nu_origPDG != 12 )
continue;
153 if( nu_ccnc != 0 )
continue;
154 if( nu_PDG != -12 )
continue;
155 if( nu_origPDG != -12 )
continue;
158 if( nu_ccnc != 1 )
continue;
159 if( nu_PDG != nu_origPDG )
continue;
163 if( nu_mode != 0 )
continue;
166 if( nu_mode != 1 )
continue;
169 if( nu_mode != 2 )
continue;
172 if( nu_mode != 10 )
continue;
175 if( nu_mode != 3 )
continue;
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;
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);
188 myhists->Xsec_Enu->Fill(nu_trueE,total_weight);
189 myhists->diffXsec_lepE->Fill(lep_E,diff_weight);
191 lepton.SetPxPyPzE(lep_Px,lep_Py,lep_Pz,lep_E);
192 myhists->diffXsec_lepcostheta->Fill(
cos(lepton.Eta()),diff_weight);
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;
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);
207 cout<<
"# selected events are "<<Nfinal<<
endl;
unsigned int nutype(unsigned int u)
virtual void Loop(const char *TypeNameNu, const char *TypeNameInt)