readFlux.C
Go to the documentation of this file.
1 #define readFlux_cxx
2 #include "readFlux.h"
3 #include <TH2.h>
4 #include <TStyle.h>
5 #include <TCanvas.h>
6 #include <iostream>
7 #include <iomanip>
8 #include "TVector3.h"
9 #include "TRandom3.h"
10 #include "TMath.h"
11 
12 using namespace std;
13 
15 {
16  if (fChain == 0) return;
17 
18  int nexperiments=100;
19  bool SuperCorrelated=true;//switch to set all weights are correlated
20  bool FullyCorrelated=false;//swtich to set all weights inside each set of pZ/xF bins to be correlated
21  bool useMIPPforSecondaryOnly=false;
22  bool useTotalErrForMIPP=false;
23 
24  Int_t nentries = fChain->GetEntries();
25  cout<<"Total events are "<<nentries<<endl;
26 
27  Float_t Ndirect=0.;
28  Float_t Nindirect=0.;
29 
30  //read the cross section or yield weights
31  //From NA49 for secondary kaon
32  //cover pT from 0.05 to 1 GeV, will treat pT<0.05 same as the first bin
33  //cover pZ up to 27 GeV
34  TFile *kaonplus_na49=new TFile("../nue_weights/weight_kaonplus_NA49.root");
35  TH1F *sf_kaonplus_na49_xf0=(TH1F*)kaonplus_na49->Get("ratio_kplus_xf_0");//-0.0125<xF<0.0125
36  TH1F *sf_kaonplus_na49_xf1=(TH1F*)kaonplus_na49->Get("ratio_kplus_xf_1");//0.0125<xF<0.0375
37  TH1F *sf_kaonplus_na49_xf2=(TH1F*)kaonplus_na49->Get("ratio_kplus_xf_2");//0.0375<xF<0.0625
38  TH1F *sf_kaonplus_na49_xf3=(TH1F*)kaonplus_na49->Get("ratio_kplus_xf_3");//0.0625<xF<0.075
39  TH1F *sf_kaonplus_na49_xf4=(TH1F*)kaonplus_na49->Get("ratio_kplus_xf_4");//0.075<xF<0.125
40  TH1F *sf_kaonplus_na49_xf5=(TH1F*)kaonplus_na49->Get("ratio_kplus_xf_5");//0.125<xF<0.175
41  TH1F *sf_kaonplus_na49_xf6=(TH1F*)kaonplus_na49->Get("ratio_kplus_xf_6");//0.175<xF<0.225
42 
43  TFile *kaonminus_na49=new TFile("../nue_weights/weight_kaonminus_NA49.root");
44  TH1F *sf_kaonminus_na49_xf0=(TH1F*)kaonminus_na49->Get("ratio_kminus_xf_0");//-0.0125<xF<0.0125
45  TH1F *sf_kaonminus_na49_xf1=(TH1F*)kaonminus_na49->Get("ratio_kminus_xf_1");//0.0125<xF<0.0375
46  TH1F *sf_kaonminus_na49_xf2=(TH1F*)kaonminus_na49->Get("ratio_kminus_xf_2");//0.0375<xF<0.0625
47  TH1F *sf_kaonminus_na49_xf3=(TH1F*)kaonminus_na49->Get("ratio_kminus_xf_3");//0.0625<xF<0.075
48  TH1F *sf_kaonminus_na49_xf4=(TH1F*)kaonminus_na49->Get("ratio_kminus_xf_4");//0.075<xF<0.125
49  TH1F *sf_kaonminus_na49_xf5=(TH1F*)kaonminus_na49->Get("ratio_kminus_xf_5");//0.125<xF<0.175
50  TH1F *sf_kaonminus_na49_xf6=(TH1F*)kaonminus_na49->Get("ratio_kminus_xf_6");//0.175<xF<0.225
51 
52  //extend secondary kaon from pZ>20GeV up to 60 GeV
53  //using MIPP k/pi ratio and NA49 pion
54  //cover pT from 0 to 2 GeV
55  TFile *kaonext_na49mipp=new TFile("../nue_weights/weight_kaon_ext_NA49MIPP.root");
56  TH1F *sf_kaonplus_ext1=(TH1F*)kaonext_na49mipp->Get("kplus_ratio_inter_1");//24<pZ<31
57  TH1F *sf_kaonplus_ext2=(TH1F*)kaonext_na49mipp->Get("kplus_ratio_inter_2");//31<pZ<42
58  TH1F *sf_kaonplus_ext3=(TH1F*)kaonext_na49mipp->Get("kplus_ratio_inter_3");//42<pZ<60
59  TH1F *sf_kaonminus_ext1=(TH1F*)kaonext_na49mipp->Get("kminus_ratio_inter_1");//24<pZ<31
60  TH1F *sf_kaonminus_ext2=(TH1F*)kaonext_na49mipp->Get("kminus_ratio_inter_2");//31<pZ<42
61  TH1F *sf_kaonminus_ext3=(TH1F*)kaonext_na49mipp->Get("kminus_ratio_inter_3");//42<pZ<60
62 
63  //////pion from NA49
64  TFile *piplus_na49=new TFile("../nue_weights/weight_pionplus_NA49.root");
65  TH1F *sf_piplus_na49_xf0=(TH1F*)piplus_na49->Get("ratio_piplus_NA49_xf000");//-0.005<xF<0.005
66  TH1F *sf_piplus_na49_xf1=(TH1F*)piplus_na49->Get("ratio_piplus_NA49_xf001");//0.005<xF<0.015
67  TH1F *sf_piplus_na49_xf2=(TH1F*)piplus_na49->Get("ratio_piplus_NA49_xf002");//0.015<xF<0.025
68  TH1F *sf_piplus_na49_xf3=(TH1F*)piplus_na49->Get("ratio_piplus_NA49_xf003");//0.025<xF<0.035
69  TH1F *sf_piplus_na49_xf4=(TH1F*)piplus_na49->Get("ratio_piplus_NA49_xf005");//0.045<xF<0.055
70 
71  TFile *piminus_na49=new TFile("../nue_weights/weight_pionminus_NA49.root");
72  TH1F *sf_piminus_na49_xf0=(TH1F*)piminus_na49->Get("ratio_piminus_NA49_xf000");
73  TH1F *sf_piminus_na49_xf1=(TH1F*)piminus_na49->Get("ratio_piminus_NA49_xf001");
74  TH1F *sf_piminus_na49_xf2=(TH1F*)piminus_na49->Get("ratio_piminus_NA49_xf002");
75  TH1F *sf_piminus_na49_xf3=(TH1F*)piminus_na49->Get("ratio_piminus_NA49_xf003");
76  TH1F *sf_piminus_na49_xf4=(TH1F*)piminus_na49->Get("ratio_piminus_NA49_xf005");
77 
78  //////using difference between FLUKA and G4NuMI as systematic uncertainty
79  //////mainly for the >1 interactions nue events
80  TFile *fmc=new TFile("../nue_weights/weight_fluka_g4numi_diff.root");
81  TH1F *hpion_ftfp=(TH1F*)fmc->Get("fluka_ftfp_muon_tot");
82  TH1F *hpion_qgsp=(TH1F*)fmc->Get("fluka_qgsp_muon_tot");
83  TH1F *hkaon_ftfp=(TH1F*)fmc->Get("fluka_ftfp_kaon_tot");
84  TH1F *hkaon_qgsp=(TH1F*)fmc->Get("fluka_qgsp_kaon_tot");
85 
86  /////using MIPP pion results
87  TFile *fmipp=new TFile("../nue_weights/mipp_pion_weight.root");
88  TH1F *sf_piplus_mipp_stat[23];
89  TH1F *sf_piplus_mipp_syst[23];
90  TH1F *sf_piplus_mipp_tot[23];
91  TH1F *sf_piminus_mipp_stat[22];
92  TH1F *sf_piminus_mipp_syst[22];
93  TH1F *sf_piminus_mipp_tot[22];
94  for( int i=0; i<23; ++i ){
95  TString plus_stat=Form("mipp_piplus_weight_stat_pz_%i",i);
96  TString plus_syst=Form("mipp_piplus_weight_syst_pz_%i",i);
97  sf_piplus_mipp_stat[i]=(TH1F*)fmipp->Get(plus_stat);
98  sf_piplus_mipp_syst[i]=(TH1F*)fmipp->Get(plus_syst);
99  TString plus_tot=Form("mipp_piplus_weight_tot_pz_%i",i);
100  sf_piplus_mipp_tot[i]=(TH1F*)fmipp->Get(plus_tot);
101  }//pi+
102  for( int i=0; i<22; ++i ){
103  TString minus_stat=Form("mipp_piminus_weight_stat_pz_%i",i);
104  TString minus_syst=Form("mipp_piminus_weight_syst_pz_%i",i);
105  sf_piminus_mipp_stat[i]=(TH1F*)fmipp->Get(minus_stat);
106  sf_piminus_mipp_syst[i]=(TH1F*)fmipp->Get(minus_syst);
107  TString minus_tot=Form("mipp_piminus_weight_tot_pz_%i",i);
108  sf_piminus_mipp_tot[i]=(TH1F*)fmipp->Get(minus_tot);
109  }//pi-
110 
111  //set random number with initial seeds
112  //TRandom3 grand(111);
113  TRandom3 grand_kplus_na49_xf0[20];
114  TRandom3 grand_kplus_na49_xf1[20];
115  TRandom3 grand_kplus_na49_xf2[20];
116  TRandom3 grand_kplus_na49_xf3[20];
117  TRandom3 grand_kplus_na49_xf4[20];
118  TRandom3 grand_kplus_na49_xf5[20];
119  TRandom3 grand_kplus_na49_xf6[20];
120 
121  TRandom3 grand_kminus_na49_xf0[20];
122  TRandom3 grand_kminus_na49_xf1[20];
123  TRandom3 grand_kminus_na49_xf2[20];
124  TRandom3 grand_kminus_na49_xf3[20];
125  TRandom3 grand_kminus_na49_xf4[20];
126  TRandom3 grand_kminus_na49_xf5[20];
127  TRandom3 grand_kminus_na49_xf6[20];
128 
129  TRandom3 grand_kplus_ext1[20];
130  TRandom3 grand_kplus_ext2[20];
131  TRandom3 grand_kplus_ext3[20];
132  TRandom3 grand_kminus_ext1[20];
133  TRandom3 grand_kminus_ext2[20];
134  TRandom3 grand_kminus_ext3[20];
135 
136  TRandom3 grand_piplus_na49_xf0[20];
137  TRandom3 grand_piplus_na49_xf1[20];
138  TRandom3 grand_piplus_na49_xf2[20];
139  TRandom3 grand_piplus_na49_xf3[20];
140  TRandom3 grand_piplus_na49_xf4[20];
141  TRandom3 grand_piminus_na49_xf0[20];
142  TRandom3 grand_piminus_na49_xf1[20];
143  TRandom3 grand_piminus_na49_xf2[20];
144  TRandom3 grand_piminus_na49_xf3[20];
145  TRandom3 grand_piminus_na49_xf4[20];
146 
147  TRandom3 grand_pion_ftfp[20];
148  TRandom3 grand_pion_qgsp[20];
149  TRandom3 grand_kaon_ftfp[20];
150  TRandom3 grand_kaon_qgsp[20];
151 
152  for( int i=0; i<20; ++i ){
153  grand_kplus_na49_xf0[i].SetSeed(100+i);
154  grand_kplus_na49_xf1[i].SetSeed(200+i);
155  grand_kplus_na49_xf2[i].SetSeed(300+i);
156  grand_kplus_na49_xf3[i].SetSeed(400+i);
157  grand_kplus_na49_xf4[i].SetSeed(500+i);
158  grand_kplus_na49_xf5[i].SetSeed(600+i);
159  grand_kplus_na49_xf6[i].SetSeed(700+i);
160 
161  grand_kminus_na49_xf0[i].SetSeed(800+i);
162  grand_kminus_na49_xf1[i].SetSeed(900+i);
163  grand_kminus_na49_xf2[i].SetSeed(1000+i);
164  grand_kminus_na49_xf3[i].SetSeed(1100+i);
165  grand_kminus_na49_xf4[i].SetSeed(1200+i);
166  grand_kminus_na49_xf5[i].SetSeed(1300+i);
167  grand_kminus_na49_xf6[i].SetSeed(1400+i);
168 
169  grand_kplus_ext1[i].SetSeed(1500+i);
170  grand_kplus_ext2[i].SetSeed(1600+i);
171  grand_kplus_ext3[i].SetSeed(1700+i);
172  grand_kminus_ext1[i].SetSeed(1800+i);
173  grand_kminus_ext2[i].SetSeed(1900+i);
174  grand_kminus_ext3[i].SetSeed(2000+i);
175 
176  grand_piplus_na49_xf0[i].SetSeed(2100+i);
177  grand_piplus_na49_xf1[i].SetSeed(2200+i);
178  grand_piplus_na49_xf2[i].SetSeed(2300+i);
179  grand_piplus_na49_xf3[i].SetSeed(2400+i);
180  grand_piplus_na49_xf4[i].SetSeed(2500+i);
181  grand_piminus_na49_xf0[i].SetSeed(2600+i);
182  grand_piminus_na49_xf1[i].SetSeed(2700+i);
183  grand_piminus_na49_xf2[i].SetSeed(2800+i);
184  grand_piminus_na49_xf3[i].SetSeed(2900+i);
185  grand_piminus_na49_xf4[i].SetSeed(3000+i);
186 
187  grand_pion_ftfp[i].SetSeed(3100+i);
188  grand_pion_qgsp[i].SetSeed(3200+i);
189  grand_kaon_ftfp[i].SetSeed(3300+i);
190  grand_kaon_qgsp[i].SetSeed(3400+i);
191  }
192 
193  TRandom3 grand_piplus_mipp_stat[23][20];
194  TRandom3 grand_piplus_mipp_syst[23][20];
195  TRandom3 grand_piplus_mipp_tot[23][20];
196  TRandom3 grand_piminus_mipp_stat[22][20];
197  TRandom3 grand_piminus_mipp_syst[22][20];
198  TRandom3 grand_piminus_mipp_tot[22][20];
199 
200  for( int i=0; i<23; ++i ){
201  for( int j=0; j<20; ++j ){
202  int seed1=(i+1)*5000+j;
203  grand_piplus_mipp_stat[i][j].SetSeed(seed1);
204  int seed2=(i+1)*5000+j+100;
205  grand_piplus_mipp_syst[i][j].SetSeed(seed2);
206  int seed3=(i+1)*5000+j+200;
207  grand_piplus_mipp_tot[i][j].SetSeed(seed3);
208  }//loop 20 pt bins
209  }//loop 23 pz bins
210 
211  for( int i=0; i<22; ++i ){
212  for( int j=0; j<20; ++j ){
213  int seed1=(i+1)*5000+j+300;
214  grand_piminus_mipp_stat[i][j].SetSeed(seed1);
215  int seed2=(i+1)*5000+j+400;
216  grand_piminus_mipp_syst[i][j].SetSeed(seed2);
217  int seed3=(i+1)*5000+500;
218  grand_piminus_mipp_tot[i][j].SetSeed(seed3);
219  }//loop 20 pt bins
220  }//loop 22 pz bins
221 
222  float mipp_pzbins[25]={0.0,0.5,0.62,0.75,0.88,1.0,1.2,1.5,2.0, 4.0,6.0,8.0,10.0,12.0,14.0,17.0,20.0,24.0, 28.0,34.0,40.0,48.0,56.0,68.0,80.0};
223 
224  TRandom3 grand_other(20151010);
225  TRandom3 grand_na49_syst(20151012);
226 
227  for( int iran=0; iran<nexperiments; ++iran ){
228  int Nevts=0;
229 
230  //gaussian random weights
231  double weight_piplus_mipp_stat[23][20]={0.};
232  double weight_piplus_mipp_syst[23][20]={0.};
233  double weight_piplus_mipp_tot[23][20]={0.};
234  double weight_piminus_mipp_stat[22][20]={0.};
235  double weight_piminus_mipp_syst[22][20]={0.};
236  double weight_piminus_mipp_tot[22][20]={0.};
237  for( int i=0; i<23; ++i ){
238  for( int j=0; j<20; ++j ){
239  weight_piplus_mipp_stat[i][j]=grand_piplus_mipp_stat[i][j].Gaus(0,1);
240  weight_piplus_mipp_syst[i][j]=grand_piplus_mipp_syst[i][j].Gaus(0,1);
241  weight_piplus_mipp_tot[i][j]=grand_piplus_mipp_tot[i][j].Gaus(0,1);
242  }
243  }//pi+
244  for( int i=0; i<23; ++i ){
245  for( int j=0; j<20; ++j ){
246  weight_piminus_mipp_stat[i][j]=grand_piminus_mipp_stat[i][j].Gaus(0,1);
247  weight_piminus_mipp_syst[i][j]=grand_piminus_mipp_syst[i][j].Gaus(0,1);
248  weight_piminus_mipp_tot[i][j]=grand_piminus_mipp_tot[i][j].Gaus(0,1);
249  }
250  }//pi-
251 
252  double weight_kplus_na49_xf0[20];
253  double weight_kplus_na49_xf1[20];
254  double weight_kplus_na49_xf2[20];
255  double weight_kplus_na49_xf3[20];
256  double weight_kplus_na49_xf4[20];
257  double weight_kplus_na49_xf5[20];
258  double weight_kplus_na49_xf6[20];
259 
260  double weight_kminus_na49_xf0[20];
261  double weight_kminus_na49_xf1[20];
262  double weight_kminus_na49_xf2[20];
263  double weight_kminus_na49_xf3[20];
264  double weight_kminus_na49_xf4[20];
265  double weight_kminus_na49_xf5[20];
266  double weight_kminus_na49_xf6[20];
267 
268  double weight_kplus_ext1[20];
269  double weight_kplus_ext2[20];
270  double weight_kplus_ext3[20];
271  double weight_kminus_ext1[20];
272  double weight_kminus_ext2[20];
273  double weight_kminus_ext3[20];
274 
275  double weight_piplus_na49_xf0[20];
276  double weight_piplus_na49_xf1[20];
277  double weight_piplus_na49_xf2[20];
278  double weight_piplus_na49_xf3[20];
279  double weight_piplus_na49_xf4[20];
280  double weight_piminus_na49_xf0[20];
281  double weight_piminus_na49_xf1[20];
282  double weight_piminus_na49_xf2[20];
283  double weight_piminus_na49_xf3[20];
284  double weight_piminus_na49_xf4[20];
285 
286  double weight_pion_ftfp[20];
287  double weight_pion_qgsp[20];
288  double weight_kaon_ftfp[20];
289  double weight_kaon_qgsp[20];
290 
291  for( int i=0; i<20; ++i ){
292  weight_kplus_na49_xf0[i]=grand_kplus_na49_xf0[i].Gaus(0,1);
293  weight_kplus_na49_xf1[i]=grand_kplus_na49_xf1[i].Gaus(0,1);
294  weight_kplus_na49_xf2[i]=grand_kplus_na49_xf2[i].Gaus(0,1);
295  weight_kplus_na49_xf3[i]=grand_kplus_na49_xf3[i].Gaus(0,1);
296  weight_kplus_na49_xf4[i]=grand_kplus_na49_xf4[i].Gaus(0,1);
297  weight_kplus_na49_xf5[i]=grand_kplus_na49_xf5[i].Gaus(0,1);
298  weight_kplus_na49_xf6[i]=grand_kplus_na49_xf6[i].Gaus(0,1);
299 
300  weight_kminus_na49_xf0[i]=grand_kminus_na49_xf0[i].Gaus(0,1);
301  weight_kminus_na49_xf1[i]=grand_kminus_na49_xf1[i].Gaus(0,1);
302  weight_kminus_na49_xf2[i]=grand_kminus_na49_xf2[i].Gaus(0,1);
303  weight_kminus_na49_xf3[i]=grand_kminus_na49_xf3[i].Gaus(0,1);
304  weight_kminus_na49_xf4[i]=grand_kminus_na49_xf4[i].Gaus(0,1);
305  weight_kminus_na49_xf5[i]=grand_kminus_na49_xf5[i].Gaus(0,1);
306  weight_kminus_na49_xf6[i]=grand_kminus_na49_xf6[i].Gaus(0,1);
307 
308  weight_kplus_ext1[i]=grand_kplus_ext1[i].Gaus(0,1);
309  weight_kplus_ext2[i]=grand_kplus_ext2[i].Gaus(0,1);
310  weight_kplus_ext3[i]=grand_kplus_ext3[i].Gaus(0,1);
311  weight_kminus_ext1[i]=grand_kminus_ext1[i].Gaus(0,1);
312  weight_kminus_ext2[i]=grand_kminus_ext2[i].Gaus(0,1);
313  weight_kminus_ext3[i]=grand_kminus_ext3[i].Gaus(0,1);
314 
315  weight_piplus_na49_xf0[i]=grand_piplus_na49_xf0[i].Gaus(0,1);
316  weight_piplus_na49_xf1[i]=grand_piplus_na49_xf1[i].Gaus(0,1);
317  weight_piplus_na49_xf2[i]=grand_piplus_na49_xf2[i].Gaus(0,1);
318  weight_piplus_na49_xf3[i]=grand_piplus_na49_xf3[i].Gaus(0,1);
319  weight_piplus_na49_xf4[i]=grand_piplus_na49_xf4[i].Gaus(0,1);
320 
321  weight_piminus_na49_xf0[i]=grand_piminus_na49_xf0[i].Gaus(0,1);
322  weight_piminus_na49_xf1[i]=grand_piminus_na49_xf1[i].Gaus(0,1);
323  weight_piminus_na49_xf2[i]=grand_piminus_na49_xf2[i].Gaus(0,1);
324  weight_piminus_na49_xf3[i]=grand_piminus_na49_xf3[i].Gaus(0,1);
325  weight_piminus_na49_xf4[i]=grand_piminus_na49_xf4[i].Gaus(0,1);
326 
327  weight_pion_ftfp[i]=grand_pion_ftfp[i].Gaus(0,1);
328  weight_pion_qgsp[i]=grand_pion_qgsp[i].Gaus(0,1);
329  weight_kaon_ftfp[i]=grand_kaon_ftfp[i].Gaus(0,1);
330  weight_kaon_qgsp[i]=grand_kaon_qgsp[i].Gaus(0,1);
331  }
332 
333  double weight_other=grand_other.Gaus(0,1);
334  double weight_na49_syst=grand_na49_syst.Gaus(0,1);
335 
336  //for( int ii=0; ii<10000; ii++ ) {
337  for( int ii=0; ii<nentries; ii++ ) {
338  fChain->GetEntry(ii);
339  ++Nevts;
340 
341  if( Nevts%1000000==0 ) cout<<Nevts<<" processed ..."<<endl;
342  if( pdg != 12 ) continue;//for nue
343 
344  if( fabs(Vx)>2. ) continue;
345  if( fabs(Vy)>2. ) continue;
346 
347  //calculate the pT and pZ, as well as xF to extract the Xsec weights
348  double pT=sqrt(tpx*tpx + tpy*tpy);
349  double pZ=tpz;
350  double xF=tpz/120.;
351  if( pT<0.03 ) continue;
352 
353  double weight=1.0;
354 
355  bool fromKaon=false;
356  bool fromMuon=false;
357  bool fromPion=false;
358  bool fromOther=false;
359  if( abs(ptype)==13 ) fromMuon=true;
360  else if( abs(ptype)==211 ) fromPion=true;
361  else if( abs(ptype)==130 || abs(ptype)==311 || abs(ptype)==321 || abs(ptype)==310 ) fromKaon=true;
362  else fromOther=true;
363  if( fromOther ) cout<<"rare process from "<<ptype<<endl;
364 
365  bool fromMuonPlus=false;
366  if( ptype==-13 ) fromMuonPlus=true;
367  bool fromKaonPlus=false;
368  if( abs(ptype)==130 || abs(ptype)==311 || ptype==321 || abs(ptype)==310 ) fromKaonPlus=true;
369 
370  if( iran==0 ){
371  if( fromMuon || fromKaon ) myhists->Enorm->Fill(E,weight);
372  myhists->Enorm_total->Fill(E,weight);
373  }//save the default energy spectrum once for the first random experiment
374 
375  if( fromMuon ) Nint = Nint-2;
376  else Nint=Nint-1;
377  //cout<<ii<<" "<<ptype<<" "<<tptype<<" "<<Nint<<endl;
378 
379  myhists->nint->Fill(Nint);
380  if( fromMuon ) myhists->nint_muon->Fill(Nint);
381  else myhists->nint_kaon->Fill(Nint);
382 
383  //cross section weight
384  float Xweight=1.;
385 
386  //for pion->muon
387  if( fromMuon ){
388  if( iran==0 ){
389  myhists->Enorm_muon->Fill(E,weight);
390  myhists->muon_tgtptype->Fill(tptype,weight);
391  if( Nint==1 ) myhists->Enorm_muon_dir->Fill(E,weight);
392  else myhists->Enorm_muon_indir->Fill(E,weight);
393  }
394 
395  //phase space covered by MIPP data
396  bool useMIPP=false;
397  if( useMIPPforSecondaryOnly ){
398  if( Nint==1 ){
399  if( pZ<0.5 && pT<0.4 ) useMIPP=true;
400  if( pZ>0.5&&pZ<2.0 && pT<0.5 ) useMIPP=true;
401  if( pZ>4.0&&pZ<6.0 && pT<0.4 ) useMIPP=true;
402  if( pZ>6.0&&pZ<8.0 && pT<0.5 ) useMIPP=true;
403  if( pZ>8.0&&pZ<68.0 && pT<2.0) useMIPP=true;
404  }//for secondary pions only
405  }//MIPP weights only used for secondary pions
406  else{
407  if( pZ<0.5 && pT<0.4 ) useMIPP=true;
408  if( pZ>0.5&&pZ<2.0 && pT<0.5 ) useMIPP=true;
409  if( pZ>4.0&&pZ<6.0 && pT<0.4 ) useMIPP=true;
410  if( pZ>6.0&&pZ<8.0 && pT<0.5 ) useMIPP=true;
411  if( pZ>8.0&&pZ<68.0 && pT<2.0) useMIPP=true;
412  }//MIPP weights used for all pions
413 
414 
415  //applying MIPP weight only for secondary pions
416  if( useMIPP ){
417  //cout<<pZ<<" "<<pT<<" "<<tptype<<endl;
418  if( pZ<0.5 && pT>0.4 ) cout<<"XBBU: "<<pZ<<" "<<pT<<endl;
419  if( fromMuonPlus ){
420  for( int ipz=0; ipz<23; ++ipz ){
421  if( ipz<8 ){
422  //cout<<ipz<<" "<<mipp_pzbins[ipz]<<" "<<mipp_pzbins[ipz+1]<<endl;
423  if( pZ>mipp_pzbins[ipz] && pZ<mipp_pzbins[ipz+1] ){
424  int Nbin=sf_piplus_mipp_stat[ipz]->FindBin(pT);
425  float bincon=sf_piplus_mipp_stat[ipz]->GetBinContent(Nbin);
426  if( bincon<=0. ) cout<<"WARNING: (MIPP pi+): "<<Nbin<<" bin has "<<bincon<<" for ipz="<<ipz<<" pZ and pT is "<<pZ<<","<<pT<<endl;
427  float staterr=1.0;
428  float systerr=1.0;
429  float toterr=1.;
430  if( bincon>0. ){
431  staterr=sf_piplus_mipp_stat[ipz]->GetBinError(Nbin)/bincon;
432  systerr=sf_piplus_mipp_syst[ipz]->GetBinError(Nbin)/bincon;
433  toterr=sf_piplus_mipp_tot[ipz]->GetBinError(Nbin)/bincon;
434  }
435  float rannum_stat=staterr*weight_piplus_mipp_stat[ipz][Nbin-1];
436  float rannum_syst=0.0;
437  float rannum_tot=0.0;
438  if( SuperCorrelated ){
439  rannum_syst=systerr*weight_piplus_mipp_syst[0][0];
440  rannum_tot=toterr*weight_piplus_mipp_tot[0][0];
441  }
442  else if( FullyCorrelated ){
443  rannum_syst=systerr*weight_piplus_mipp_syst[ipz][0];
444  rannum_tot=toterr*weight_piplus_mipp_tot[ipz][0];
445  }//fully-correlated
446  else{
447  rannum_syst=systerr*weight_piplus_mipp_syst[ipz][Nbin-1];
448  rannum_tot=toterr*weight_piplus_mipp_tot[ipz][Nbin-1];
449  }//not-correlated
450  if( useTotalErrForMIPP ){
451  Xweight=bincon*(1.0+rannum_tot);
452  }//
453  else{
454  //10% syst from mock data test
455  float other_syst=0.1*weight_other;
456  Xweight=bincon*(1.0+rannum_stat+rannum_syst+other_syst);
457  }//
458  }//within pz bins
459  }//pZ<2GeV
460  else{
461  //cout<<ipz<<" "<<mipp_pzbins[ipz+1]<<" "<<mipp_pzbins[ipz+2]<<endl;
462  if( pZ>mipp_pzbins[ipz+1] && pZ<mipp_pzbins[ipz+2] ){
463  int Nbin=sf_piplus_mipp_stat[ipz]->FindBin(pT);
464  float bincon=sf_piplus_mipp_stat[ipz]->GetBinContent(Nbin);
465  if( bincon<=0. ) cout<<"WARNING: (MIPP pi+): "<<Nbin<<" bin has "<<bincon<<" for ipz="<<ipz<<" pZ and pT is "<<pZ<<","<<pT<<endl;
466  float staterr=1.0;
467  float systerr=1.0;
468  float toterr=1.;
469  if( bincon>0. ){
470  staterr=sf_piplus_mipp_stat[ipz]->GetBinError(Nbin)/bincon;
471  systerr=sf_piplus_mipp_syst[ipz]->GetBinError(Nbin)/bincon;
472  toterr=sf_piplus_mipp_tot[ipz]->GetBinError(Nbin)/bincon;
473  }
474  float rannum_stat=staterr*weight_piplus_mipp_stat[ipz][Nbin-1];
475  float rannum_syst=0.0;
476  float rannum_tot=0.;
477  if( SuperCorrelated ){
478  rannum_syst=systerr*weight_piplus_mipp_syst[0][0];
479  rannum_tot=toterr*weight_piplus_mipp_tot[0][0];
480  }
481  else if( FullyCorrelated ){
482  rannum_syst=systerr*weight_piplus_mipp_syst[ipz][0];
483  rannum_tot=toterr*weight_piplus_mipp_tot[ipz][0];
484  }//fully-correlated
485  else{
486  rannum_syst=systerr*weight_piplus_mipp_syst[ipz][Nbin-1];
487  rannum_tot=toterr*weight_piplus_mipp_tot[ipz][Nbin-1];
488  }//not-correlated
489  if( useTotalErrForMIPP ){
490  Xweight=bincon*(1.0+rannum_tot);
491  }
492  else{
493  float other_syst=0.1*weight_other;
494  Xweight=bincon*(1.0+rannum_stat+rannum_syst+other_syst);
495  }
496  }//within pz bins
497  }//pZ>4GeV
498  }//loop pz bins
499 
500  }//for pi+
501  else{
502  for( int ipz=0; ipz<22; ++ipz ){
503  if( ipz<8 ){
504  //cout<<ipz<<" "<<mipp_pzbins[ipz]<<" "<<mipp_pzbins[ipz+1]<<endl;
505  if( pZ>mipp_pzbins[ipz] && pZ<mipp_pzbins[ipz+1] ){
506  int Nbin=sf_piminus_mipp_stat[ipz]->FindBin(pT);
507  float bincon=sf_piminus_mipp_stat[ipz]->GetBinContent(Nbin);
508  if( bincon<=0. ) cout<<"WARNING: (MIPP pi-): "<<Nbin<<" bin has "<<bincon<<" for ipz="<<ipz<<endl;
509  float staterr=1.0;
510  float systerr=1.0;
511  if( bincon>0. ){
512  staterr=sf_piminus_mipp_stat[ipz]->GetBinError(Nbin)/bincon;
513  systerr=sf_piminus_mipp_syst[ipz]->GetBinError(Nbin)/bincon;
514  }
515  float rannum_stat=staterr*weight_piminus_mipp_stat[ipz][Nbin-1];
516  float rannum_syst=0.0;
517  if( SuperCorrelated ){
518  rannum_syst=systerr*weight_piminus_mipp_syst[0][0];
519  }
520  else if( FullyCorrelated ){
521  rannum_syst=systerr*weight_piminus_mipp_syst[ipz][0];
522  }//fully-correlated
523  else{
524  rannum_syst=systerr*weight_piminus_mipp_syst[ipz][Nbin-1];
525  }//not-correlated
526  float other_syst=0.1*weight_other;
527  Xweight=bincon*(1.0+rannum_stat+rannum_syst+other_syst);
528  }//within pz bins
529  }//pZ<2GeV
530  else{
531  //cout<<ipz<<" "<<mipp_pzbins[ipz+1]<<" "<<mipp_pzbins[ipz+2]<<endl;
532  if( pZ>mipp_pzbins[ipz+1] && pZ<mipp_pzbins[ipz+2] ){
533  int Nbin=sf_piminus_mipp_stat[ipz]->FindBin(pT);
534  float bincon=sf_piminus_mipp_stat[ipz]->GetBinContent(Nbin);
535  if( bincon<=0. ) cout<<"WARNING: (MIPP pi-): "<<Nbin<<" bin has "<<bincon<<" for ipz="<<ipz<<endl;
536  float staterr=1.0;
537  float systerr=1.0;
538  if( bincon>0. ){
539  staterr=sf_piminus_mipp_stat[ipz]->GetBinError(Nbin)/bincon;
540  systerr=sf_piminus_mipp_syst[ipz]->GetBinError(Nbin)/bincon;
541  }
542  float rannum_stat=staterr*weight_piminus_mipp_stat[ipz][Nbin-1];
543  float rannum_syst=0.0;
544  if( SuperCorrelated ){
545  rannum_syst=systerr*weight_piminus_mipp_syst[0][0];
546  }
547  else if( FullyCorrelated ){
548  rannum_syst=systerr*weight_piminus_mipp_syst[ipz][0];
549  }//fully-correlated
550  else{
551  rannum_syst=systerr*weight_piminus_mipp_syst[ipz][Nbin-1];
552  }//not-correlated
553  float other_syst=0.1*weight_other;
554  Xweight=bincon*(1.0+rannum_stat+rannum_syst+other_syst);
555  }//within pz bins
556  }//pZ>4GeV
557  }//loop pz bins
558  }//for pi-
559  }//use MIPP data
560 
561  bool useNA49=false;//use pion data from NA49
562  if( !useMIPP && Nint==1 ){
563  //if( 0.015<xF<0.035 && pT<0.75 ) useNA49=true;
564  //if( 0.005<xF<0.015 && 0.5<pT<0.75 ) useNA49=true;
565  //if( xF<0.005 && 0.4<pT<2.0 ) useNA49=true;
566  //if( 0.045<xF<0.055 && 0.4<pT<2.0 ) useNA49=true;
567  //conver the xF to pZ
568  if( pZ<0.5 && pT>0.4 ) useNA49=true;
569  if( pZ>0.5&&pZ<0.6 && pT>0.5 ) useNA49=true;
570  if( pZ>0.6&&pZ<2.0 && 0.5<pT<0.75 ) useNA49=true;
571  if( pZ>2.0&&pZ<4.0 && pT<0.75 ) useNA49=true;
572  if( pZ>4.0&&pZ<4.2 && 0.4<pT<0.75 ) useNA49=true;
573  if( pZ>5.4&&pZ<6.0 && pT>0.4 ) useNA49=true;
574  if( pZ>6.0&&pZ<6.6 && pT>0.5 ) useNA49=true;
575  }//for secondary pions not covered by MIPP
576 
577  if( useNA49 ){
578  float systerr=0.038*weight_na49_syst;
579  if( (pZ<0.5 && pT>0.4) || (pZ>0.5&&pZ<0.6 && pT>0.5) ){
580  if( fromMuonPlus ){
581  int Nbin=sf_piplus_na49_xf0->FindBin(pT);
582  float bincon=sf_piplus_na49_xf0->GetBinContent(Nbin);
583  if( bincon<=0. ) cout<<"WARNING: (NA49 pi+): "<<Nbin<<" bin has "<<bincon<<" for xf0"<<endl;
584  float binerr=1.0;
585  if( bincon>0. ) binerr=sf_piplus_na49_xf0->GetBinError(Nbin)/bincon;
586  binerr=binerr*weight_piplus_na49_xf0[Nbin-1];
587  Xweight=bincon*(1.0+binerr+systerr);
588  }//pi+
589  else{
590  int Nbin=sf_piminus_na49_xf0->FindBin(pT);
591  float bincon=sf_piminus_na49_xf0->GetBinContent(Nbin);
592  if( bincon<=0. ) cout<<"WARNING: (NA49 pi-): "<<Nbin<<" bin has "<<bincon<<" for xf0"<<endl;
593  float binerr=1.0;
594  if( bincon>0. ) binerr=sf_piminus_na49_xf0->GetBinError(Nbin)/bincon;
595  binerr=binerr*weight_piminus_na49_xf0[Nbin-1];
596  Xweight=bincon*(1.0+binerr+systerr);
597  }//pi-
598  }//-0.005<xF<0.005
599  if( (pZ>0.6&&pZ<2.0) && pT>0.5&&pT<0.75 ){
600  if( fromMuonPlus ){
601  int Nbin=sf_piplus_na49_xf1->FindBin(pT);
602  float bincon=sf_piplus_na49_xf1->GetBinContent(Nbin);
603  if( bincon<=0. ) cout<<"WARNING: (NA49 pi+): "<<Nbin<<" bin has "<<bincon<<" for xf1"<<endl;
604  float binerr=1.0;
605  if( bincon>0. ) binerr=sf_piplus_na49_xf1->GetBinError(Nbin)/bincon;
606  binerr=binerr*weight_piplus_na49_xf1[Nbin-1];
607  Xweight=bincon*(1.0+binerr+systerr);
608  }//pi+
609  else{
610  int Nbin=sf_piminus_na49_xf1->FindBin(pT);
611  float bincon=sf_piminus_na49_xf1->GetBinContent(Nbin);
612  if( bincon<=0. ) cout<<"WARNING: (NA49 pi-): "<<Nbin<<" bin has "<<bincon<<" for xf1"<<endl;
613  float binerr=1.0;
614  if( bincon>0. ) binerr=sf_piminus_na49_xf1->GetBinError(Nbin)/bincon;
615  binerr=binerr*weight_piminus_na49_xf1[Nbin-1];
616  Xweight=bincon*(1.0+binerr+systerr);
617  }//pi-
618  }//0.005<xF<0.015
619  if( pZ>2.0&&pZ<3.0 && pT<0.75 ){
620  if( fromMuonPlus ){
621  int Nbin=sf_piplus_na49_xf2->FindBin(pT);
622  float bincon=sf_piplus_na49_xf2->GetBinContent(Nbin);
623  if( bincon<=0. ) cout<<"WARNING: (NA49 pi+): "<<Nbin<<" bin has "<<bincon<<" for xf2"<<", pZ and pT is "<<pZ<<" "<<pT<<endl;
624  float binerr=1.0;
625  if( bincon>0. ) binerr=sf_piplus_na49_xf2->GetBinError(Nbin)/bincon;
626  binerr=binerr*weight_piplus_na49_xf2[Nbin-1];
627  Xweight=bincon*(1.0+binerr+systerr);
628  }//pi+
629  else{
630  int Nbin=sf_piminus_na49_xf2->FindBin(pT);
631  float bincon=sf_piminus_na49_xf2->GetBinContent(Nbin);
632  if( bincon<=0. ) cout<<"WARNING: (NA49 pi-): "<<Nbin<<" bin has "<<bincon<<" for xf2"<<", pZ and pT is "<<pZ<<" "<<pT<<endl;
633  float binerr=1.0;
634  if( bincon>0. ) binerr=sf_piminus_na49_xf2->GetBinError(Nbin)/bincon;
635  binerr=binerr*weight_piminus_na49_xf2[Nbin-1];
636  Xweight=bincon*(1.0+binerr+systerr);
637  }//pi-
638  }//0.015<xF<0.025
639  if( (pZ>3.0&&pZ<4.0 && pT<0.75) || (pZ>4.0&&pZ<4.2 && pZ>0.4&&pT<0.75) ){
640  if( fromMuonPlus ){
641  int Nbin=sf_piplus_na49_xf3->FindBin(pT);
642  float bincon=sf_piplus_na49_xf3->GetBinContent(Nbin);
643  if( bincon<=0. ) cout<<"WARNING: (NA49 pi+): "<<Nbin<<" bin has "<<bincon<<" for xf3"<<endl;
644  float binerr=1.0;
645  if( bincon>0. ) binerr=sf_piplus_na49_xf3->GetBinError(Nbin)/bincon;
646  binerr=binerr*weight_piplus_na49_xf3[Nbin-1];
647  Xweight=bincon*(1.0+binerr+systerr);
648  }//pi+
649  else{
650  int Nbin=sf_piminus_na49_xf3->FindBin(pT);
651  float bincon=sf_piminus_na49_xf3->GetBinContent(Nbin);
652  if( bincon<=0. ) cout<<"WARNING: (NA49 pi-): "<<Nbin<<" bin has "<<bincon<<" for xf3"<<endl;
653  float binerr=1.0;
654  if( bincon>0. ) binerr=sf_piminus_na49_xf3->GetBinError(Nbin)/bincon;
655  binerr=binerr*weight_piminus_na49_xf3[Nbin-1];
656  Xweight=bincon*(1.0+binerr+systerr);
657  }//pi-
658  }//0.025<xF<0.035
659  if( (pZ>5.4&&pZ<6.0 && pT>0.4)||(pZ>6.0&&pZ<6.6 && pT>0.5) ){
660  if( fromMuonPlus ){
661  int Nbin=sf_piplus_na49_xf4->FindBin(pT);
662  float bincon=sf_piplus_na49_xf4->GetBinContent(Nbin);
663  if( bincon<=0. ) cout<<"WARNING: (NA49 pi+): "<<Nbin<<" bin has "<<bincon<<" for xf4"<<endl;
664  float binerr=1.0;
665  if( bincon>0. ) binerr=sf_piplus_na49_xf4->GetBinError(Nbin)/bincon;
666  binerr=binerr*weight_piplus_na49_xf4[Nbin-1];
667  Xweight=bincon*(1.0+binerr+systerr);
668  }//pi+
669  else{
670  int Nbin=sf_piminus_na49_xf4->FindBin(pT);
671  float bincon=sf_piminus_na49_xf4->GetBinContent(Nbin);
672  if( bincon<=0. ) cout<<"WARNING: (NA49 pi-): "<<Nbin<<" bin has "<<bincon<<" for xf4"<<endl;
673  float binerr=1.0;
674  if( bincon>0. ) binerr=sf_piminus_na49_xf4->GetBinError(Nbin)/bincon;
675  binerr=binerr*weight_piminus_na49_xf4[Nbin-1];
676  Xweight=bincon*(1.0+binerr+systerr);
677  }//pi-
678  }//0.045<xF<0.055
679  }//using NA49 data
680 
681  bool useMC=false;//use difference between FLUKA and G4NuMI directly
682  if( !useNA49 && !useMIPP ){
683  //for the moment do not split Nint=1 and Nint>1 interactions
684  //using ftfp physics list for G4NuMI
685  int Nbin=hpion_ftfp->FindBin(E);
686  //float binerr=hpion_ftfp->GetBinContent(Nbin);
687  Xweight=1.0+0.2*weight_pion_ftfp[0];
688 
689  }//for tertiary and other pions not covered by MIPP
690 
691  myhists->Eweight_muon[iran]->Fill(E,weight*Xweight);
692  if( Nint==1 ) myhists->Eweight_muon_dir[iran]->Fill(E,weight*Xweight);
693  else myhists->Eweight_muon_indir[iran]->Fill(E,weight*Xweight);
694  //cout<<"From muon: "<<E<<" with weight "<<Xweight<<", pT,pZ and xF: "<<pT<<" "<<pZ<<" "<<xF<<endl;
695  }//from Muon
696 
697  //for kaon
698  if( fromKaon ){
699  if( iran==0 ){
700  myhists->Enorm_kaon->Fill(E,weight);
701  myhists->kaon_tgtptype->Fill(tptype,weight);
702  if( Nint==1 ) myhists->Enorm_kaon_dir->Fill(E,weight);
703  else myhists->Enorm_kaon_indir->Fill(E,weight);
704  }
705 
706  if( Nint==1 ){
707  if( xF<0.225 ){
708  if( pT<1. ){//using NA49 Kaon data with pT upto 1 GeV
709  if(xF<0.0125){
710  if( fromKaonPlus ){
711  int Nbin=sf_kaonplus_na49_xf0->FindBin(pT);
712  if(pT<0.1) Nbin=1;//use first bin for pT<0.1 GeV
713  if(pT>0.9 && pT<1. ) Nbin=9;//use last bin for 0.9<pT<1 GeV
714  float bincon=sf_kaonplus_na49_xf0->GetBinContent(Nbin);
715  if( bincon<=0. ) cout<<"WARNING: (NA49 K+): "<<Nbin<<" bin has "<<bincon<<" for xf0"<<endl;
716  float binerr=1.0;
717  if( bincon>0. ) binerr=sf_kaonplus_na49_xf0->GetBinError(Nbin)/bincon;
718  if( SuperCorrelated ){
719  binerr=binerr*weight_kplus_na49_xf0[0];
720  }
721  else if( FullyCorrelated ){
722  binerr=binerr*weight_kplus_na49_xf0[0];
723  }//fully-correlated
724  else{
725  binerr=binerr*weight_kplus_na49_xf0[Nbin-1];
726  }//no correlation
727  Xweight=bincon*(1.0+binerr);
728  }//from K+
729  else{
730  int Nbin=sf_kaonminus_na49_xf0->FindBin(pT);
731  if(pT<0.1) Nbin=1;//use first bin for pT<0.1 GeV
732  if(pT>0.9 && pT<1. ) Nbin=9;//use last bin for 0.9<pT<1 GeV
733  float bincon=sf_kaonminus_na49_xf0->GetBinContent(Nbin);
734  if( bincon<=0. ) cout<<"WARNING: (NA49 K-): "<<Nbin<<" bin has "<<bincon<<" for xf0"<<endl;
735  float binerr=1.0;
736  if( bincon>0. ) binerr=sf_kaonminus_na49_xf0->GetBinError(Nbin)/bincon;
737  if( SuperCorrelated ){
738  binerr=binerr*weight_kminus_na49_xf0[0];
739  }
740  else if( FullyCorrelated ){
741  binerr=binerr*weight_kminus_na49_xf0[0];
742  }//fully-correlated
743  else{
744  binerr=binerr*weight_kminus_na49_xf0[Nbin-1];
745  }//no correlation
746  Xweight=bincon*(1.0+binerr);
747  }//from K-
748  }//xF<0.015
749  else if(xF<0.0375){
750  if( fromKaonPlus ){
751  int Nbin=sf_kaonplus_na49_xf1->FindBin(pT);
752  if(pT<0.1) Nbin=1;//use first bin for pT<0.1 GeV
753  if(pT>0.9 && pT<1. ) Nbin=9;//use last bin for 0.9<pT<1 GeV
754  float bincon=sf_kaonplus_na49_xf1->GetBinContent(Nbin);
755  if( bincon<=0. ) cout<<"WARNING: (NA49 K+): "<<Nbin<<" bin has "<<bincon<<" for xf1"<<endl;
756  float binerr=1.0;
757  if( bincon>0. ) binerr=sf_kaonplus_na49_xf1->GetBinError(Nbin)/bincon;
758  if( SuperCorrelated ){
759  binerr=binerr*weight_kplus_na49_xf0[0];
760  }
761  else if( FullyCorrelated ){
762  binerr=binerr*weight_kplus_na49_xf1[0];
763  }//fully-correlated
764  else{
765  binerr=binerr*weight_kplus_na49_xf1[Nbin-1];
766  }//no correlation
767  Xweight=bincon*(1.0+binerr);
768  }//from K+
769  else{
770  int Nbin=sf_kaonminus_na49_xf1->FindBin(pT);
771  if(pT<0.1) Nbin=1;//use first bin for pT<0.1 GeV
772  if(pT>0.9 && pT<1. ) Nbin=9;//use last bin for 0.9<pT<1 GeV
773  float bincon=sf_kaonminus_na49_xf1->GetBinContent(Nbin);
774  if( bincon<=0. ) cout<<"WARNING: (NA49 K-): "<<Nbin<<" bin has "<<bincon<<" for xf1"<<endl;
775  float binerr=1.0;
776  if( bincon>0. ) binerr=sf_kaonminus_na49_xf1->GetBinError(Nbin)/bincon;
777  if( SuperCorrelated ){
778  binerr=binerr*weight_kminus_na49_xf0[0];
779  }
780  else if( FullyCorrelated ){
781  binerr=binerr*weight_kminus_na49_xf1[0];
782  }//fully-correlated
783  else{
784  binerr=binerr*weight_kminus_na49_xf1[Nbin-1];
785  }//no correlation
786  Xweight=bincon*(1.0+binerr);
787  }//from K-
788  }//0.0125<xF<0.0375
789  else if(xF<0.0625){
790  if( fromKaonPlus ){
791  int Nbin=sf_kaonplus_na49_xf2->FindBin(pT);
792  if(pT<0.1) Nbin=1;//use first bin for pT<0.1 GeV
793  if(pT>0.9 && pT<1. ) Nbin=9;//use last bin for 0.9<pT<1 GeV
794  float bincon=sf_kaonplus_na49_xf2->GetBinContent(Nbin);
795  if( bincon<=0. ) cout<<"WARNING: (NA49 K+): "<<Nbin<<" bin has "<<bincon<<" for xf2"<<endl;
796  float binerr=1.0;
797  if( bincon>0. ) binerr=sf_kaonplus_na49_xf2->GetBinError(Nbin)/bincon;
798  if( SuperCorrelated ){
799  binerr=binerr*weight_kplus_na49_xf0[0];
800  }
801  else if( FullyCorrelated ){
802  binerr=binerr*weight_kplus_na49_xf2[0];
803  }//fully-correlated
804  else{
805  binerr=binerr*weight_kplus_na49_xf2[Nbin-1];
806  }//no correlation
807  Xweight=bincon*(1.0+binerr);
808  }//from K+
809  else{
810  int Nbin=sf_kaonminus_na49_xf2->FindBin(pT);
811  if(pT<0.1) Nbin=1;//use first bin for pT<0.1 GeV
812  if(pT>0.9 && pT<1. ) Nbin=9;//use last bin for 0.9<pT<1 GeV
813  float bincon=sf_kaonminus_na49_xf2->GetBinContent(Nbin);
814  if( bincon<=0. ) cout<<"WARNING: (NA49 K-): "<<Nbin<<" bin has "<<bincon<<" for xf2"<<endl;
815  float binerr=1.0;
816  if( bincon>0.0 ) binerr=sf_kaonminus_na49_xf2->GetBinError(Nbin)/bincon;
817  if( SuperCorrelated ){
818  binerr=binerr*weight_kminus_na49_xf0[0];
819  }
820  else if( FullyCorrelated ){
821  binerr=binerr*weight_kminus_na49_xf2[0];
822  }//fully-correlated
823  else{
824  binerr=binerr*weight_kminus_na49_xf2[Nbin-1];
825  }//no correlation
826  Xweight=bincon*(1.0+binerr);
827  }//from K-
828  }//0.0375<xF<0.0625
829  else if(xF<0.075){
830  if( fromKaonPlus ){
831  int Nbin=sf_kaonplus_na49_xf3->FindBin(pT);
832  if(pT<0.1) Nbin=1;//use first bin for pT<0.1 GeV
833  if(pT>0.9 && pT<1. ) Nbin=9;//use last bin for 0.9<pT<1 GeV
834  float bincon=sf_kaonplus_na49_xf3->GetBinContent(Nbin);
835  if( bincon<=0. ) cout<<"WARNING: (NA49 K+): "<<Nbin<<" bin has "<<bincon<<" for xf3"<<endl;
836  float binerr=1.0;
837  if( bincon>0. ) binerr=sf_kaonplus_na49_xf3->GetBinError(Nbin)/bincon;
838  if( SuperCorrelated ){
839  binerr=binerr*weight_kplus_na49_xf0[0];
840  }
841  else if( FullyCorrelated ){
842  binerr=binerr*weight_kplus_na49_xf3[0];
843  }//fully-correlated
844  else{
845  binerr=binerr*weight_kplus_na49_xf3[Nbin-1];
846  }//no correlation
847  Xweight=bincon*(1.0+binerr);
848  }//from K+
849  else{
850  int Nbin=sf_kaonminus_na49_xf3->FindBin(pT);
851  if(pT<0.1) Nbin=1;//use first bin for pT<0.1 GeV
852  if(pT>0.9 && pT<1. ) Nbin=9;//use last bin for 0.9<pT<1 GeV
853  float bincon=sf_kaonminus_na49_xf3->GetBinContent(Nbin);
854  if( bincon<=0. ) cout<<"WARNING: (NA49 K-): "<<Nbin<<" bin has "<<bincon<<" for xf3"<<endl;
855  float binerr=1.0;
856  if( bincon>0. ) binerr=sf_kaonminus_na49_xf3->GetBinError(Nbin)/bincon;
857  if( SuperCorrelated ){
858  binerr=binerr*weight_kminus_na49_xf0[0];
859  }
860  else if( FullyCorrelated ){
861  binerr=binerr*weight_kminus_na49_xf3[0];
862  }//fully-correlated
863  else{
864  binerr=binerr*weight_kminus_na49_xf3[Nbin-1];
865  }//no correlation
866  Xweight=bincon*(1.0+binerr);
867  }//from K-
868  }//0.0625<xF<0.075
869  else if(xF<0.125){
870  if( fromKaonPlus ){
871  int Nbin=sf_kaonplus_na49_xf4->FindBin(pT);
872  if(pT<0.1) Nbin=1;//use first bin for pT<0.1 GeV
873  if(pT>0.9 && pT<1. ) Nbin=9;//use last bin for 0.9<pT<1 GeV
874  float bincon=sf_kaonplus_na49_xf4->GetBinContent(Nbin);
875  if( bincon<=0. ) cout<<"WARNING: (NA49 K+): "<<Nbin<<" bin has "<<bincon<<" for xf4"<<endl;
876  float binerr=1.0;
877  if( bincon>0. ) binerr=sf_kaonplus_na49_xf4->GetBinError(Nbin)/bincon;
878  if( SuperCorrelated ){
879  binerr=binerr*weight_kplus_na49_xf0[0];
880  }
881  else if( FullyCorrelated ){
882  binerr=binerr*weight_kplus_na49_xf4[0];
883  }//fully-correlated
884  else{
885  binerr=binerr*weight_kplus_na49_xf4[Nbin-1];
886  }//no correlation
887  Xweight=bincon*(1.0+binerr);
888  }//from K+
889  else{
890  int Nbin=sf_kaonminus_na49_xf4->FindBin(pT);
891  if(pT<0.1) Nbin=1;//use first bin for pT<0.1 GeV
892  if(pT>0.9 && pT<1. ) Nbin=9;//use last bin for 0.9<pT<1 GeV
893  float bincon=sf_kaonminus_na49_xf4->GetBinContent(Nbin);
894  if( bincon<=0. ) cout<<"WARNING: (NA49 K-): "<<Nbin<<" bin has "<<bincon<<" for xf4"<<endl;
895  float binerr=1.0;
896  if( bincon>0. ) binerr=sf_kaonminus_na49_xf4->GetBinError(Nbin)/bincon;
897  if( SuperCorrelated ){
898  binerr=binerr*weight_kminus_na49_xf0[0];
899  }
900  else if( FullyCorrelated ){
901  binerr=binerr*weight_kminus_na49_xf4[0];
902  }//fully-correlated
903  else{
904  binerr=binerr*weight_kminus_na49_xf4[Nbin-1];
905  }//no correlation
906  Xweight=bincon*(1.0+binerr);
907  }//from K-
908  }//0.075<xF<0.125
909  else if(xF<0.175){
910  if( fromKaonPlus ){
911  int Nbin=sf_kaonplus_na49_xf5->FindBin(pT);
912  if(pT<0.1) Nbin=1;//use first bin for pT<0.1 GeV
913  if(pT>0.9 && pT<1. ) Nbin=9;//use last bin for 0.9<pT<1 GeV
914  float bincon=sf_kaonplus_na49_xf5->GetBinContent(Nbin);
915  if( bincon<=0. ) cout<<"WARNING: (NA49 K+): "<<Nbin<<" bin has "<<bincon<<" for xf5"<<endl;
916  float binerr=1.0;
917  if( bincon>0. ) binerr=sf_kaonplus_na49_xf5->GetBinError(Nbin)/bincon;
918  if( SuperCorrelated ){
919  binerr=binerr*weight_kplus_na49_xf0[0];
920  }
921  else if( FullyCorrelated ){
922  binerr=binerr*weight_kplus_na49_xf5[0];
923  }//fully-correlated
924  else{
925  binerr=binerr*weight_kplus_na49_xf5[Nbin-1];
926  }//no correlation
927  Xweight=bincon*(1.0+binerr);
928  }//from K+
929  else{
930  int Nbin=sf_kaonminus_na49_xf5->FindBin(pT);
931  if(pT<0.1) Nbin=1;//use first bin for pT<0.1 GeV
932  if(pT>0.9 && pT<1. ) Nbin=9;//use last bin for 0.9<pT<1 GeV
933  float bincon=sf_kaonminus_na49_xf5->GetBinContent(Nbin);
934  if( bincon<=0. ) cout<<"WARNING: (NA49 K-): "<<Nbin<<" bin has "<<bincon<<" for xf5"<<endl;
935  float binerr=1.0;
936  if( bincon>0. ) binerr=sf_kaonminus_na49_xf5->GetBinError(Nbin);
937  if( SuperCorrelated ){
938  binerr=binerr*weight_kminus_na49_xf0[0];
939  }
940  else if( FullyCorrelated ){
941  binerr=binerr*weight_kminus_na49_xf5[0];
942  }//fully-correlated
943  else{
944  binerr=binerr*weight_kminus_na49_xf5[Nbin-1];
945  }//no correlation
946  Xweight=bincon*(1.0+binerr);
947  }//from K-
948  }//0.125<xF<0.175
949  else{
950  if( fromKaonPlus ){
951  int Nbin=sf_kaonplus_na49_xf6->FindBin(pT);
952  if(pT<0.1) Nbin=1;//use first bin for pT<0.1 GeV
953  if(pT>0.9 && pT<1. ) Nbin=9;//use last bin for 0.9<pT<1 GeV
954  float bincon=sf_kaonplus_na49_xf6->GetBinContent(Nbin);
955  if( bincon<=0. ) cout<<"WARNING: (NA49 K+): "<<Nbin<<" bin has "<<bincon<<" for xf6"<<endl;
956  float binerr=1.0;
957  if( bincon>0. ) binerr=sf_kaonplus_na49_xf6->GetBinError(Nbin)/bincon;
958  if( SuperCorrelated ){
959  binerr=binerr*weight_kplus_na49_xf0[0];
960  }
961  else if( FullyCorrelated ){
962  binerr=binerr*weight_kplus_na49_xf6[0];
963  }//fully-correlated
964  else{
965  binerr=binerr*weight_kplus_na49_xf6[Nbin-1];
966  }//no correlation
967  Xweight=bincon*(1.0+binerr);
968  }//from K+
969  else{
970  int Nbin=sf_kaonminus_na49_xf6->FindBin(pT);
971  if(pT<0.1) Nbin=1;//use first bin for pT<0.1 GeV
972  if(pT>0.9 && pT<1. ) Nbin=9;//use last bin for 0.9<pT<1 GeV
973  float bincon=sf_kaonminus_na49_xf6->GetBinContent(Nbin);
974  if( bincon<=0. ) cout<<"WARNING: (NA49 K-): "<<Nbin<<" bin has "<<bincon<<" for xf6"<<endl;
975  float binerr=1.0;
976  if( bincon>0. ) binerr=sf_kaonminus_na49_xf6->GetBinError(Nbin)/bincon;
977  if( SuperCorrelated ){
978  binerr=binerr*weight_kminus_na49_xf0[0];
979  }
980  else if( FullyCorrelated ){
981  binerr=binerr*weight_kminus_na49_xf6[0];
982  }//fully-correlated
983  else{
984  binerr=binerr*weight_kminus_na49_xf6[Nbin-1];
985  }//no correlation
986  Xweight=bincon*(1.0+binerr);
987  }//from K-
988  }//0.175<xF<0.225
989  }//pT<1.0 GeV
990  //else{
991  //using difference between FLUKA and G4NuMI
992  //}
993  }//xF<0.225 (pZ<27GeV)
994  if( xF<0.5&&pT<2.0 ){ //using extended kaon data, pT coverage upto 2GeV
995  if( pZ<31 ){
996  if( fromKaonPlus ){
997  int Nbin=sf_kaonplus_ext1->FindBin(pT);
998  float bincon=sf_kaonplus_ext1->GetBinContent(Nbin);
999  if( bincon<=0. ) cout<<"WARNING: (NA49 K+): "<<Nbin<<" bin has "<<bincon<<" for ext1"<<endl;
1000  float binerr=1.0;
1001  if( bincon>0. ) binerr=sf_kaonplus_ext1->GetBinError(Nbin)/bincon;
1002  if( SuperCorrelated ){
1003  binerr=binerr*weight_kplus_ext1[0];
1004  }
1005  else if( FullyCorrelated ){
1006  binerr=binerr*weight_kplus_ext1[0];
1007  }//fully-correlated
1008  else{
1009  binerr=binerr*weight_kplus_ext1[Nbin-1];
1010  }//no correlation
1011  Xweight=bincon*(1.0+binerr);
1012  }//from K+
1013  else{
1014  int Nbin=sf_kaonminus_ext1->FindBin(pT);
1015  float bincon=sf_kaonminus_ext1->GetBinContent(Nbin);
1016  if( bincon<=0. ) cout<<"WARNING: (NA49 K-): "<<Nbin<<" bin has "<<bincon<<" for ext1"<<endl;
1017  float binerr=1.0;
1018  if( bincon>0. ) binerr=sf_kaonminus_ext1->GetBinError(Nbin)/bincon;
1019  if( SuperCorrelated ){
1020  binerr=binerr*weight_kminus_ext1[0];
1021  }
1022  else if( FullyCorrelated ){
1023  binerr=binerr*weight_kminus_ext1[0];
1024  }//fully-correlated
1025  else{
1026  binerr=binerr*weight_kminus_ext1[Nbin-1];
1027  }//no correlation
1028  Xweight=bincon*(1.0+binerr);
1029  }//from K-
1030  }//27<pZ<31GeV
1031  else if( pZ<42 ){
1032  if( fromKaonPlus ){
1033  int Nbin=sf_kaonplus_ext2->FindBin(pT);
1034  float bincon=sf_kaonplus_ext2->GetBinContent(Nbin);
1035  if( bincon<=0. ) cout<<"WARNING: (NA49 K+): "<<Nbin<<" bin has "<<bincon<<" for ext2"<<endl;
1036  float binerr=1.0;
1037  if( bincon>0. ) binerr=sf_kaonplus_ext2->GetBinError(Nbin)/bincon;
1038  if( SuperCorrelated ){
1039  binerr=binerr*weight_kplus_ext1[0];
1040  }
1041  else if( FullyCorrelated ){
1042  binerr=binerr*weight_kplus_ext2[0];
1043  }//fully-correlated
1044  else{
1045  binerr=binerr*weight_kplus_ext2[Nbin-1];
1046  }//no correlation
1047  Xweight=bincon*(1.0+binerr);
1048  }//from K+
1049  else{
1050  int Nbin=sf_kaonminus_ext2->FindBin(pT);
1051  float bincon=sf_kaonminus_ext2->GetBinContent(Nbin);
1052  if( bincon<=0. ) cout<<"WARNING: (NA49 K-): "<<Nbin<<" bin has "<<bincon<<" for ext2"<<endl;
1053  float binerr=1.0;
1054  if( bincon>0. ) binerr=sf_kaonminus_ext2->GetBinError(Nbin)/bincon;
1055  if( SuperCorrelated ){
1056  binerr=binerr*weight_kminus_ext1[0];
1057  }
1058  else if( FullyCorrelated ){
1059  binerr=binerr*weight_kminus_ext2[0];
1060  }//fully-correlated
1061  else{
1062  binerr=binerr*weight_kminus_ext2[Nbin-1];
1063  }//no correlation
1064  Xweight=bincon*(1.0+binerr);
1065  }//from K-
1066  }//31<pZ<42GeV
1067  else{
1068  if( fromKaonPlus ){
1069  int Nbin=sf_kaonplus_ext3->FindBin(pT);
1070  float bincon=sf_kaonplus_ext3->GetBinContent(Nbin);
1071  if( bincon<=0. ) cout<<"WARNING: (NA49 K+): "<<Nbin<<" bin has "<<bincon<<" for ext3"<<endl;
1072  float binerr=1.0;
1073  if( bincon>0. ) binerr=sf_kaonplus_ext3->GetBinError(Nbin)/bincon;
1074  if( SuperCorrelated ){
1075  binerr=binerr*weight_kplus_ext1[0];
1076  }
1077  else if( FullyCorrelated ){
1078  binerr=binerr*weight_kplus_ext3[0];
1079  }//fully-correlated
1080  else{
1081  binerr=binerr*weight_kplus_ext3[Nbin-1];
1082  }//no correlation
1083  Xweight=bincon*(1.0+binerr);
1084  }//from K+
1085  else{
1086  int Nbin=sf_kaonminus_ext3->FindBin(pT);
1087  float bincon=sf_kaonminus_ext3->GetBinContent(Nbin);
1088  if( bincon<=0. ) cout<<"WARNING: (NA49 K-): "<<Nbin<<" bin has "<<bincon<<" for ext3"<<endl;
1089  float binerr=1.0;
1090  if( bincon>0. ) binerr=sf_kaonminus_ext3->GetBinError(Nbin)/bincon;
1091  if( SuperCorrelated ){
1092  binerr=binerr*weight_kminus_ext1[0];
1093  }
1094  else if( FullyCorrelated ){
1095  binerr=binerr*weight_kminus_ext3[0];
1096  }//fully-correlated
1097  else{
1098  binerr=binerr*weight_kminus_ext3[Nbin-1];
1099  }//no correlation
1100  Xweight=bincon*(1.0+binerr);
1101  }//from K-
1102  }//42<pT<60GeV
1103 
1104  }//27<pZ<60GeV && pT>2GeV
1105 
1106  }//for secondary kaons
1107 
1108  //uncovered region for kaon
1109  if( Nint!=1 || (Nint==1&&pZ>60.) ||(Nint==1&&pZ<27.0&&pT>1.0) || (Nint==1&&pZ>27.&&pZ<60.&&pT>2.0) ){
1110  //using difference between FLUKA and G4NuMI
1111  //for the moment do not split Nint=1 and Nint>1 interactions
1112  //using ftfp physics list for G4NuMI
1113  int Nbin=hkaon_ftfp->FindBin(E);
1114  float binerr=hkaon_ftfp->GetBinContent(Nbin);
1115  //if( binerr<0.1 ) binerr=0.1;
1116  if( binerr<0.2 ) binerr=0.2;
1117  Xweight=1.0+weight_kaon_ftfp[0]*binerr;
1118  }//for tertiary and other kaons
1119 
1120  myhists->Eweight_kaon[iran]->Fill(E,weight*Xweight);
1121  if( Nint==1 ) myhists->Eweight_kaon_dir[iran]->Fill(E,weight*Xweight);
1122  else myhists->Eweight_kaon_indir[iran]->Fill(E,weight*Xweight);
1123  //cout<<"From kaon: "<<E<<" with weight "<<Xweight<<", pT,pZ and xF: "<<pT<<" "<<pZ<<" "<<xF<<endl;
1124  }//from Kaon
1125 
1126  if( fromMuon || fromKaon )
1127  myhists->Eweight[iran]->Fill(E,weight*Xweight);
1128 
1129  myhists->Eweight_total[iran]->Fill(E,weight*Xweight);
1130  }//loop all events
1131  cout<<"Random experiment "<<iran<<", total events are "<<Nevts<<"."<<endl;
1132  }//end of 1000 random experiments
1133 
1134 }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const Var weight
T sqrt(T number)
Definition: d0nt_math.hpp:156
TFile * fmc
Definition: bdt_com.C:14
float abs(float number)
Definition: d0nt_math.hpp:39
Float_t E
Definition: plot.C:20
Long64_t nentries
TChain * fChain
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
virtual void Loop()
Definition: readFlux.C:14