00001 #ifndef madanalysis_cxx
00002 #define madanalysis_cxx
00003 #include <iostream>
00004 #include "TClonesArray.h"
00005 #include "TMath.h"
00006 #include "Validity/VldContext.h"
00007 #include "Conventions/Detector.h"
00008 #include "Registry/Registry.h"
00009 #include "Mad/MadAnalysis.h"
00010 #include "MCReweight/MCReweight.h"
00011 #include "MCReweight/GnumiInterface.h"
00012 #include "MCReweight/NuParent.h"
00013 #include "AnalysisNtuples/ANtpHeaderInfo.h"
00014 #include "AnalysisNtuples/ANtpTruthInfoBeam.h"
00015 #include "AnalysisNtuples/ANtpAnalysisInfo.h"
00016 #include "AnalysisNtuples/ANtpRecoInfo.h"
00017 #include "AnalysisNtuples/ANtpEventInfo.h"
00018 #include "AnalysisNtuples/ANtpTrackInfo.h"
00019 #include "AnalysisNtuples/ANtpShowerInfo.h"
00020 #include "AnalysisNtuples/Module/ANtpRecoNtpManipulator.h"
00021 #include "AnalysisNtuples/Module/ANtpInfoObjectFillerBeam.h"
00022 #include "BeamDataUtil/BeamMonSpill.h"
00023 #include "BeamDataUtil/BMSpillAna.h"
00024 #include "BeamDataUtil/BDSpillAccessor.h"
00025 #include "SpillTiming/SpillTimeFinder.h"
00026
00027 using namespace std;
00028
00029
00030
00031
00032 const double FARPOTTONUMEVT = 668.4;
00033 const double NEARPOTTONUMEVT = 100000;
00034
00035
00036 MadAnalysis::MadAnalysis(TChain *chainSR,TChain *chainMC,
00037 TChain *chainTH,TChain *chainEM)
00038 :MadQuantities(chainSR,chainMC,chainTH,chainEM)
00039 {
00040
00041 NumberOfBins = 1000;
00042 RangeLowerLimit = 0.;
00043 RangeUpperLimit = 50.;
00044
00045 MCSignalHist[0] = NULL; MCSignalHist[1] = NULL;
00046 MCBkgdHist[0] = NULL; MCBkgdHist[1] = NULL;
00047 DataHist[0] = NULL; DataHist[1] = NULL;
00048 DataSignalHist[0] = NULL; DataSignalHist[1] = NULL;
00049 DataPOT[0] = 0; DataPOT[1] = 0;
00050 numDataEvts[0] = 0; numDataEvts[1] = 0;
00051 signalFracInData[0] = 0; signalFracInData[1] = 0;
00052 DataBkgdHist[0] = NULL; DataBkgdHist[1] = NULL;
00053 MCEvents[0].clear(); MCEvents[1].clear();
00054 MCPOT[0] = 0; MCPOT[1] = 0;
00055 numMCEvts[0] = 0; numMCEvts[1] = 0;
00056 BestFit[0] = NULL; BestFit[1] = NULL;
00057
00058 OscillationFunction = NULL;
00059 OscillationParameters = NULL;
00060 NumOscPars = 0;
00061 varyX = new int[2];
00062 varyX[0] = 1;
00063 varyX[1] = 2;
00064
00065 NgenBins[0] = 5; NgenBins[1] = 0; NgenBins[2] = 0;
00066 NgenMin[0] = 0.5; NgenMin[1] = 0.; NgenMin[2] = 0.;
00067 NgenMax[0] = 1.5; NgenMax[1] = 0.; NgenMax[2] = 0.;
00068 NgenMean[0] = 1.; NgenMean[1] = 0.; NgenMean[2] = 0.;
00069 NgenErr[0] = 0.1; NgenErr[1] = 0.; NgenErr[2] = 0.;
00070 NgenName[0] = "neugen:ma_qe"; NgenName[1] = "empty";
00071 NgenName[2] = "empty";
00072
00073 tag.append("NoTag");
00074
00075 Chi2Calc = new MadChi2Calc(0);
00076 Chi2Surf = NULL;
00077 ChiMin=99999.;
00078 NDOF = 0;
00079 Par1MinVal = 0.;
00080 Par2MinVal = 0.;
00081
00082 EShiftErr = 0.;
00083
00084 writeOut = true;
00085
00086 if(!chainSR) {
00087 record = 0;
00088 strecord = 0;
00089 emrecord = 0;
00090 mcrecord = 0;
00091 threcord = 0;
00092 Clear();
00093 whichSource = -1;
00094 isMC = true;
00095 isTH = true;
00096 isEM = true;
00097 Nentries = -1;
00098 return;
00099 }
00100
00101
00102 whichSource = 0;
00103
00104 }
00105
00106
00107 MadAnalysis::MadAnalysis(JobC *j,string path,int entries)
00108 : MadQuantities(j,path,entries)
00109 {
00110
00111 NumberOfBins = 1000;
00112 RangeLowerLimit = 0.;
00113 RangeUpperLimit = 50.;
00114
00115 MCSignalHist[0] = NULL; MCSignalHist[1] = NULL;
00116 MCBkgdHist[0] = NULL; MCBkgdHist[1] = NULL;
00117 DataHist[0] = NULL; DataHist[1] = NULL;
00118 DataSignalHist[0] = NULL; DataSignalHist[1] = NULL;
00119 DataPOT[0] = 0; DataPOT[1] = 0;
00120 numDataEvts[0] = 0; numDataEvts[1] = 0;
00121 signalFracInData[0] = 0; signalFracInData[1] = 0;
00122 DataBkgdHist[0] = NULL; DataBkgdHist[1] = NULL;
00123 MCEvents[0].clear(); MCEvents[1].clear();
00124 MCPOT[0] = 0; MCPOT[1] = 0;
00125 numMCEvts[0] = 0; numMCEvts[1] = 0;
00126 BestFit[0] = NULL; BestFit[1] = NULL;
00127
00128 OscillationFunction = NULL;
00129 OscillationParameters = NULL;
00130 NumOscPars = 0;
00131 varyX = new int[2];
00132 varyX[0] = 1;
00133 varyX[1] = 2;
00134
00135 tag.append("NoTag");
00136
00137 Chi2Calc = new MadChi2Calc(0);
00138 Chi2Surf = NULL;
00139 ChiMin=99999.;
00140 NDOF = 0;
00141 Par1MinVal = 0.;
00142 Par2MinVal = 0.;
00143
00144 EShiftErr = 0.;
00145
00146 writeOut = true;
00147
00148 record = 0;
00149 strecord = 0;
00150 emrecord = 0;
00151 mcrecord = 0;
00152 threcord = 0;
00153 Clear();
00154 isMC = true;
00155 isTH = true;
00156 isEM = true;
00157 Nentries = entries;
00158 whichSource = 1;
00159 jcPath = path;
00160 fJC = j;
00161 fChain = NULL;
00162
00163 }
00164
00165
00166 MadAnalysis::~MadAnalysis()
00167 {
00168
00169 if(writeOut) {
00170 EndJob();
00171 cout << "Done EndJob()" << endl;
00172 }
00173 delete Chi2Calc;
00174 delete Chi2Surf;
00175 delete MCSignalHist[0]; delete MCSignalHist[1];
00176 delete MCBkgdHist[0]; delete MCBkgdHist[1];
00177 delete DataHist[0]; delete DataHist[1];
00178 delete DataSignalHist[0]; delete DataSignalHist[1];
00179 delete DataBkgdHist[0]; delete DataBkgdHist[1];
00180 delete BestFit[0]; delete BestFit[1];
00181 delete OscillationFunction;
00182 delete [] varyX;
00183 cout << "deleting MadAnalysis object" << endl;
00184 }
00185
00186
00187 void MadAnalysis::ReInit(TChain *chainSR,TChain *chainMC,
00188 TChain *chainTH,TChain *chainEM){
00189
00190 InitChain(chainSR,chainMC,chainTH,chainEM);
00191
00192 }
00193
00194
00195 void MadAnalysis::ReInit(JobC *j,string path,int entries){
00196
00197 record = 0;
00198 strecord = 0;
00199 emrecord = 0;
00200 mcrecord = 0;
00201 threcord = 0;
00202 Clear();
00203 isMC = true;
00204 isTH = true;
00205 isEM = true;
00206 Nentries = entries;
00207 whichSource = 1;
00208 jcPath = path;
00209 fJC = j;
00210
00211 }
00212
00213
00214 void MadAnalysis::SetFileNameTag(std::string t)
00215 {
00216 tag = t;
00217 }
00218
00219
00220 void MadAnalysis::SetHistBookInfo(int bins,float low,float up){
00221
00222 NumberOfBins = bins;
00223 RangeLowerLimit = low;
00224 RangeUpperLimit = up;
00225
00226 }
00227
00228
00229 int MadAnalysis::SetDataInfo(TH1F *hist, int numev, float pot, float sigfrac,
00230 int det){
00231
00232 if(hist) {
00233 DataHist[det] = new TH1F(*hist);
00234 numDataEvts[det] = numev;
00235 DataPOT[det] = pot;
00236 signalFracInData[det] = sigfrac;
00237 return 1;
00238 }
00239 return 0;
00240
00241 }
00242
00243
00244 void MadAnalysis::CreatePAN(std::string tag){
00245
00246 this->GetEntry(0);
00247 const bool file_is_mc
00248 =(ntpHeader->GetVldContext().GetSimFlag()==SimFlag::kMC);
00249
00250 std::string savename = "PAN_" + tag + ".root";
00251 TFile save(savename.c_str(),"RECREATE");
00252 save.cd();
00253
00254 GnumiInterface *gnumi = 0;
00255 if(file_is_mc){
00256 gnumi = new GnumiInterface();
00257 if(!gnumi->Status()) {
00258 cout << "Error MadAnalysis::CreatePAN() - Flux files not open."
00259 << " Will not be filling NuParent info"
00260 << endl;
00261 cout << "Set environmental variable $GNUMIAUX to point to the "
00262 << "directory containing the gnumi files to fix this!"
00263 << endl;
00264 }
00265 else {
00266 const char* gnumidir= getenv("GNUMIAUX");
00267 cout<<"Read gnumi files from "<<gnumidir<<endl;
00268 }
00269 }
00270
00271
00272 BMSpillAna bmana;
00273
00274
00275
00276 Float_t true_enu;
00277 Float_t true_pxnu;
00278 Float_t true_pynu;
00279 Float_t true_pznu;
00280 Float_t true_etgt;
00281 Float_t true_pxtgt;
00282 Float_t true_pytgt;
00283 Float_t true_pztgt;
00284 Float_t true_emu;
00285 Float_t true_elep;
00286 Float_t true_eshw;
00287 Float_t true_x;
00288 Float_t true_y;
00289 Float_t true_q2;
00290 Float_t true_w2;
00291 Float_t true_emfrac;
00292 Float_t true_dircosneu;
00293 Float_t true_dircosz;
00294 Float_t true_vtxx;
00295 Float_t true_vtxy;
00296 Float_t true_vtxz;
00297
00298
00299 Float_t th_evt_pur;
00300 Float_t th_evt_all;
00301 Float_t th_evt_slc;
00302 Float_t th_shw_pur;
00303 Float_t th_shw_all;
00304 Float_t th_shw_slc;
00305 Float_t th_trk_pur;
00306 Float_t th_trk_all;
00307 Float_t th_trk_slc;
00308
00309
00310 Int_t flavour;
00311 Int_t nooscflavour;
00312 Int_t process;
00313 Int_t nucleus;
00314 Int_t cc_nc;
00315 Int_t initial_state;
00316 Int_t had_fs;
00317
00318
00319 NuParent *nuparent = new NuParent();
00320
00321
00322 Int_t detector;
00323 Int_t run;
00324 Int_t snarl;
00325 Int_t subrun;
00326 UInt_t trgsrc;
00327 double trgtime;
00328 Int_t spilltype;
00329
00330
00331 Int_t nevent;
00332 Int_t event;
00333 Int_t slice;
00334 Int_t mcevent;
00335 Int_t ntrack;
00336 Int_t nshower;
00337 double litime;
00338 UChar_t dmx_stat;
00339
00340 Int_t is_fid;
00341 Int_t is_cev;
00342 Int_t is_mc;
00343 Int_t pass_basic;
00344 Float_t pid0;
00345 Float_t pid1;
00346 Float_t pid2;
00347 Int_t pass;
00348
00349 Float_t reco_enu;
00350 Float_t reco_emu;
00351 Float_t reco_eshw;
00352 Float_t reco_eshw_linCC;
00353 Float_t reco_eshw_wtCC;
00354 Float_t reco_eshw_linNC;
00355 Float_t reco_eshw_wtNC;
00356 Float_t reco_eshw_em;
00357 Float_t reco_qe_enu;
00358 Float_t reco_qe_q2;
00359 Float_t reco_y;
00360 Float_t reco_dircosneu;
00361 Float_t reco_dircosz;
00362 Float_t reco_vtxx;
00363 Float_t reco_vtxy;
00364 Float_t reco_vtxz;
00365 Float_t reco_emfrac;
00366 Float_t reco_emfrac_prob;
00367
00368 Float_t evtlength;
00369 Float_t trklength;
00370 Float_t shwlength;
00371 Float_t slclength;
00372 Double_t closestevent_s;
00373 Double_t closestevent_t;
00374 Float_t fracRehitStrip;
00375
00376 Float_t trkmomrange;
00377 Float_t trkqp;
00378 Float_t trkeqp;
00379 Float_t trkphfrac;
00380 Float_t trkphplane;
00381
00382 Double_t shwstptimemin;
00383 Double_t shwstptimemax;
00384 Double_t shwstptimerms;
00385
00386
00387 Double_t hbw,vbw,hpos1,vpos1,hpos2,vpos2,hornI,closestspill;
00388
00389
00390
00391 UInt_t beamcnf;
00392
00393 Double_t nuTarZ;
00394
00395 Int_t goodbeam;
00396
00397
00398 UInt_t horn_on, target_in, n_batches;
00399 Float_t tor101, tr101d, tortgt, trtgtd;
00400 Float_t hadmon, mumon1, mumon2, mumon3;
00401
00402
00403 TTree *tree = new TTree("pan","pan");
00404
00405 tree->Branch("true_enu",&true_enu,"true_enu/F",32000);
00406 tree->Branch("true_pxnu",&true_pxnu,"true_pxnu/F",32000);
00407 tree->Branch("true_pynu",&true_pynu,"true_pynu/F",32000);
00408 tree->Branch("true_pznu",&true_pznu,"true_pznu/F",32000);
00409 tree->Branch("true_etgt",&true_etgt,"true_etgt/F",32000);
00410 tree->Branch("true_pxtgt",&true_pxtgt,"true_pxtgt/F",32000);
00411 tree->Branch("true_pytgt",&true_pytgt,"true_pytgt/F",32000);
00412 tree->Branch("true_pztgt",&true_pztgt,"true_pztgt/F",32000);
00413 tree->Branch("true_emu",&true_emu,"true_emu/F",32000);
00414 tree->Branch("true_elep",&true_elep,"true_elep/F",32000);
00415 tree->Branch("true_eshw",&true_eshw,"true_eshw/F",32000);
00416 tree->Branch("true_x",&true_x,"true_x/F",32000);
00417 tree->Branch("true_y",&true_y,"true_y/F",32000);
00418 tree->Branch("true_q2",&true_q2,"true_q2/F",32000);
00419 tree->Branch("true_w2",&true_w2,"true_w2/F",32000);
00420 tree->Branch("true_dircosneu",&true_dircosneu,"true_dircosneu/F",32000);
00421 tree->Branch("true_dircosz",&true_dircosz,"true_dircosz/F",32000);
00422 tree->Branch("true_vtxx",&true_vtxx,"true_vtxx/F",32000);
00423 tree->Branch("true_vtxy",&true_vtxy,"true_vtxy/F",32000);
00424 tree->Branch("true_vtxz",&true_vtxz,"true_vtxz/F",32000);
00425 tree->Branch("true_emfrac",&true_emfrac,"true_emfrac/F",32000);
00426
00427 tree->Branch("th_evt_pur",&th_evt_pur,"th_evt_pur/F",32000);
00428 tree->Branch("th_evt_all",&th_evt_all,"th_evt_all/F",32000);
00429 tree->Branch("th_evt_slc",&th_evt_slc,"th_evt_slc/F",32000);
00430 tree->Branch("th_shw_pur",&th_shw_pur,"th_shw_pur/F",32000);
00431 tree->Branch("th_shw_all",&th_shw_all,"th_shw_all/F",32000);
00432 tree->Branch("th_shw_slc",&th_shw_slc,"th_shw_slc/F",32000);
00433 tree->Branch("th_trk_pur",&th_trk_pur,"th_trk_pur/F",32000);
00434 tree->Branch("th_trk_all",&th_trk_all,"th_trk_all/F",32000);
00435 tree->Branch("th_trk_slc",&th_trk_slc,"th_trk_slc/F",32000);
00436
00437 tree->Branch("stdhep_ng",&fNumFSGeantino,"stdhep_ng/I",32000);
00438 tree->Branch("stdhep_eg",&fGeantinoEnergy,"stdhep_eg/F",32000);
00439 tree->Branch("stdhep_np",&fNumFSProton,"stdhep_np/I",32000);
00440 tree->Branch("stdhep_ep",&fTotFSProtonEnergy,"stdhep_ep/F",32000);
00441 tree->Branch("stdhep_mp",&fMaxFSProtonEnergy,"stdhep_mp/F",32000);
00442 tree->Branch("stdhep_nn",&fNumFSNeutron,"stdhep_nn/I",32000);
00443 tree->Branch("stdhep_en",&fTotFSNeutronEnergy,"stdhep_en/F",32000);
00444 tree->Branch("stdhep_mn",&fMaxFSNeutronEnergy,"stdhep_mn/F",32000);
00445 tree->Branch("stdhep_npp",&fNumFSPiPlus,"stdhep_npp/I",32000);
00446 tree->Branch("stdhep_epp",&fTotFSPiPlusEnergy,"stdhep_epp/F",32000);
00447 tree->Branch("stdhep_mpp",&fMaxFSPiPlusEnergy,"stdhep_mpp/F",32000);
00448 tree->Branch("stdhep_npm",&fNumFSPiMinus,"stdhep_npm/I",32000);
00449 tree->Branch("stdhep_epm",&fTotFSPiMinusEnergy,"stdhep_epm/F",32000);
00450 tree->Branch("stdhep_mpm",&fMaxFSPiMinusEnergy,"stdhep_mpm/F",32000);
00451 tree->Branch("stdhep_npz",&fNumFSPiZero,"stdhep_npz/I",32000);
00452 tree->Branch("stdhep_epz",&fTotFSPiZeroEnergy,"stdhep_epz/F",32000);
00453 tree->Branch("stdhep_mpz",&fMaxFSPiZeroEnergy,"stdhep_mpz/F",32000);
00454 tree->Branch("stdhep_nk",&fNumFSKaon,"stdhep_nk/I",32000);
00455 tree->Branch("stdhep_ek",&fTotFSKaonEnergy,"stdhep_ek/F",32000);
00456 tree->Branch("stdhep_mk",&fMaxFSKaonEnergy,"stdhep_mk/F",32000);
00457
00458 tree->Branch("flavour",&flavour,"flavour/I",32000);
00459 tree->Branch("nooscflavour",&nooscflavour,"nooscflavour/I",32000);
00460 tree->Branch("process",&process,"process/I",32000);
00461 tree->Branch("nucleus",&nucleus,"nucleus/I",32000);
00462 tree->Branch("cc_nc",&cc_nc,"cc_nc/I",32000);
00463 tree->Branch("initial_state",&initial_state,"initial_state/I",32000);
00464 tree->Branch("had_fs",&had_fs,"had_fs/I",32000);
00465
00466 tree->Branch("nuparent","NuParent",&nuparent,8000,1);
00467
00468 tree->Branch("detector",&detector,"detector/I",32000);
00469 tree->Branch("run",&run,"run/I",32000);
00470 tree->Branch("snarl",&snarl,"snarl/I",32000);
00471 tree->Branch("subrun",&subrun,"subrun/I",32000);
00472 tree->Branch("trgsrc",&trgsrc,"trgsrc/i",32000);
00473 tree->Branch("trgtime",&trgtime,"trgtime/D",32000);
00474 tree->Branch("spilltype",&spilltype,"spiltype/I",32000);
00475 tree->Branch("event",&event,"event/I",32000);
00476 tree->Branch("slice",&slice,"slice/I",32000);
00477 tree->Branch("mcevent",&mcevent,"mcevent/I",32000);
00478 tree->Branch("ntrack",&ntrack,"ntrack/I",32000);
00479 tree->Branch("nshower",&nshower,"nshower/I",32000);
00480 tree->Branch("litime",&litime,"litime/D");
00481 tree->Branch("dmx_stat",&dmx_stat,"dmx_stat/b");
00482
00483 tree->Branch("is_fid",&is_fid,"is_fid/I",32000);
00484 tree->Branch("is_cev",&is_cev,"is_cev/I",32000);
00485 tree->Branch("is_mc",&is_mc,"is_mc/I",32000);
00486 tree->Branch("pass_basic",&pass_basic,"pass_basic/I",32000);
00487 tree->Branch("pid0",&pid0,"pid0/F",32000);
00488 tree->Branch("pid1",&pid1,"pid1/F",32000);
00489 tree->Branch("pid2",&pid2,"pid2/F",32000);
00490 tree->Branch("pass",&pass,"pass/I",32000);
00491
00492 tree->Branch("reco_enu",&reco_enu,"reco_enu/F",32000);
00493 tree->Branch("reco_emu",&reco_emu,"reco_emu/F",32000);
00494 tree->Branch("reco_eshw",&reco_eshw,"reco_eshw/F",32000);
00495 tree->Branch("reco_eshw_linCC",&reco_eshw_linCC,"reco_eshw_linCC/F",32000);
00496 tree->Branch("reco_eshw_wtCC",&reco_eshw_wtCC,"reco_eshw_wtCC/F",32000);
00497 tree->Branch("reco_eshw_linNC",&reco_eshw_linNC,"reco_eshw_linNC/F",32000);
00498 tree->Branch("reco_eshw_wtNC",&reco_eshw_wtNC,"reco_eshw_wtNC/F",32000);
00499 tree->Branch("reco_eshw_em",&reco_eshw_em,"reco_eshw_em/F",32000);
00500 tree->Branch("reco_qe_enu",&reco_qe_enu,"reco_qe_enu/F",32000);
00501 tree->Branch("reco_qe_q2",&reco_qe_q2,"reco_qe_q2/F",32000);
00502 tree->Branch("reco_y",&reco_y,"reco_y/F",32000);
00503 tree->Branch("reco_dircosneu",&reco_dircosneu,"reco_dircosneu/F",32000);
00504 tree->Branch("reco_dircosz",&reco_dircosz,"reco_dircosz/F",32000);
00505 tree->Branch("reco_vtxx",&reco_vtxx,"reco_vtxx/F",32000);
00506 tree->Branch("reco_vtxy",&reco_vtxy,"reco_vtxy/F",32000);
00507 tree->Branch("reco_vtxz",&reco_vtxz,"reco_vtxz/F",32000);
00508 tree->Branch("reco_emfrac",&reco_emfrac,"reco_emfrac/F",32000);
00509 tree->Branch("reco_emfrac_prob",&reco_emfrac_prob,
00510 "reco_emfrac_prob/F",32000);
00511
00512 tree->Branch("evtlength",&evtlength,"evtlength/F",32000);
00513 tree->Branch("trklength",&trklength,"trklength/F",32000);
00514 tree->Branch("shwlength",&shwlength,"shwlength/F",32000);
00515 tree->Branch("slclength",&slclength,"slclength/F",32000);
00516 tree->Branch("closestevent_s",&closestevent_s,"closestevent_s/D",32000);
00517 tree->Branch("closestevent_t",&closestevent_t,"closestevent_t/D",32000);
00518 tree->Branch("fracRehitStrip",&fracRehitStrip,"fracRehitStrip/F",32000);
00519
00520 tree->Branch("trkmomrange",&trkmomrange,"trkmomrange/F",32000);
00521 tree->Branch("trkqp",&trkqp,"trkqp/F",32000);
00522 tree->Branch("trkeqp",&trkeqp,"trkeqp/F",32000);
00523 tree->Branch("trkphfrac",&trkphfrac,"trkphfrac/F",32000);
00524 tree->Branch("trkphplane",&trkphplane,"trkphplane/F",32000);
00525
00526 tree->Branch("shwstptimemin",&shwstptimemin,"shwstptimemin/D",32000);
00527 tree->Branch("shwstptimemax",&shwstptimemax,"shwstptimemax/D",32000);
00528 tree->Branch("shwstptimerms",&shwstptimerms,"shwstptimerms/D",32000);
00529
00530 tree->Branch("hbw",&hbw,"hbw/D");
00531 tree->Branch("vbw",&vbw,"vbw/D");
00532 tree->Branch("hpos1",&hpos1,"hpos1/D");
00533 tree->Branch("vpos1",&vpos1,"vpos1/D");
00534 tree->Branch("hpos2",&hpos2,"hpos2/D");
00535 tree->Branch("vpos2",&vpos2,"vpos2/D");
00536 tree->Branch("hornI",&hornI,"hornI/D");
00537 tree->Branch("closestspill",&closestspill,"closestspill/D");
00538 tree->Branch("beamcnf",&beamcnf,"beamcnf/i");
00539 tree->Branch("nuTarZ",&nuTarZ,"nuTarZ/D");
00540 tree->Branch("goodbeam",&goodbeam,"goodbeam/I");
00541 tree->Branch("horn_on",&horn_on,"horn_on/i");
00542 tree->Branch("target_in",&target_in,"target_in/i");
00543 tree->Branch("n_batches",&n_batches,"n_batches/i");
00544 tree->Branch("tor101",&tor101,"tor101/F");
00545 tree->Branch("tr101d",&tr101d,"tr101d/F");
00546 tree->Branch("tortgt",&tortgt,"tortgt/F");
00547 tree->Branch("trtgtd",&trtgtd,"trtgtd/F");
00548 tree->Branch("hadmon",&hadmon,"hadmon/F");
00549 tree->Branch("mumon1",&mumon1,"mumon1/F");
00550 tree->Branch("mumon2",&mumon2,"mumon2/F");
00551 tree->Branch("mumon3",&mumon3,"mumon3/F");
00552
00553 this->AddUserBranches(tree);
00554
00555
00556 float pottor101=0.;
00557 float pottortgt=0.;
00558 float pot=0.;
00559 int NRUNS=0;
00560 int NSNARLS=0;
00561 TTree *pottree = new TTree("pottree","pottree");
00562 pottree->Branch("pot",&pot,"pot/F");
00563 pottree->Branch("pottor101",&pottor101,"pottor101/F");
00564 pottree->Branch("pottortgt",&pottortgt,"pottortgt/F");
00565 pottree->Branch("nruns",&NRUNS,"nruns/I");
00566 pottree->Branch("nsnarls",&NSNARLS,"nsnarls/I");
00567
00568 int lastrun=-1;
00569 for(int i=0;i<Nentries;i++){
00570 if(i%1000==0) std::cout << float(i)*100./float(Nentries)
00571 << "% done" << std::endl;
00572
00573 if(!this->GetEntry(i)) continue;
00574 NSNARLS++;
00575 nevent = eventSummary->nevent;
00576 run = ntpHeader->GetRun();
00577 if(run!=lastrun){
00578 NRUNS++;
00579 lastrun = run;
00580 }
00581 snarl = ntpHeader->GetSnarl();
00582 trgsrc = ntpHeader->GetTrigSrc();
00583 spilltype = ntpHeader->GetRemoteSpillType();
00584 trgtime = eventSummary->trigtime;
00585 litime = eventSummary->litime;
00586 dmx_stat=0;
00587 dmx_stat+=(dmxStatus->nonphysicalfail);
00588 dmx_stat+=((dmxStatus->validplanesfail<<1));
00589 dmx_stat+=((dmxStatus->vertexplanefail<<2));
00590
00592
00593 hbw = vbw = hpos1 = vpos1 = hpos2 = vpos2 = 0.0;
00594 nuTarZ = hornI = closestspill = 0.0;
00595 beamcnf = horn_on = target_in = n_batches = 0;
00596 goodbeam=0;
00597 tor101 = tr101d = tortgt = trtgtd = 0.0;
00598 hadmon = mumon1 = mumon2 = mumon3 = 0.0;
00599
00600 if(!file_is_mc){
00601 VldContext vc = ntpHeader->GetVldContext();
00602 VldTimeStamp vts = vc.GetTimeStamp();
00603 BDSpillAccessor &sa = BDSpillAccessor::Get();
00604 const BeamMonSpill *bs = sa.LoadSpill(vts);
00605 hbw=bs->fProfWidX*1000.;
00606 vbw=bs->fProfWidY*1000.;
00607 hpos1=bs->fTargProfX*1000.;
00608 vpos1=bs->fTargProfY*1000.;
00609 n_batches=0;
00610 float btint=0.;
00611 hpos2=0.;
00612 vpos2=0.;
00613 for(int bt=0;bt<6;bt++){
00614 hpos2+=(bs->fTargBpmX[bt]*bs->fBpmInt[bt]);
00615 vpos2+=(bs->fTargBpmY[bt]*bs->fBpmInt[bt]);
00616 if(bs->fBpmInt[bt]>0.001){
00617 n_batches++;
00618 }
00619 btint+=bs->fBpmInt[bt];
00620 }
00621 if(btint!=0){
00622 hpos2=hpos2*1000./btint;
00623 vpos2=vpos2*1000./btint;
00624 }
00625 else{
00626 hpos2=-999.;
00627 vpos2=-999.;
00628 }
00629
00630 hornI=bs->fHornCur;
00631 SpillTimeFinder &sfind= SpillTimeFinder::Instance();
00632 closestspill = sfind.GetTimeToNearestSpill(vc);
00633
00634 beamcnf =bs->GetStatusBits().beam_type;
00635 nuTarZ=0.;
00636 horn_on = bs->GetStatusBits().horn_on;
00637 target_in = bs->GetStatusBits().target_in;
00638
00639 tor101=bs->fTor101;
00640 tr101d=bs->fTr101d;
00641 tortgt=bs->fTortgt;
00642 trtgtd=bs->fTrtgtd;
00643 hadmon=bs->fHadInt;
00644 mumon1=bs->fMuInt1;
00645 mumon2=bs->fMuInt2;
00646 mumon3=bs->fMuInt3;
00647
00648 goodbeam=0;
00649
00650
00651 bmana.SetSpill(*bs);
00652 bmana.SetTimeDiff(closestspill);
00653 if(bmana.SelectSpill()){
00654 goodbeam=1;
00655 }
00656 else{
00657 goodbeam=0;
00658 }
00659
00660 if(goodbeam==1&&spilltype!=3){
00661 pot+=tortgt;
00662 pottor101+=tor101;
00663 pottortgt+=tortgt;
00664 }
00665 }
00666
00667 std::vector<Double_t> evtx(nevent,0);
00668 std::vector<Double_t> evty(nevent,0);
00669 std::vector<Double_t> evtz(nevent,0);
00670 std::vector<Double_t> evtt(nevent,0);
00671
00672 for(int j=0;j<nevent;j++){
00673 if(!LoadEvent(j)) {
00674 evtx[j] = -100.;
00675 evty[j] = -100.;
00676 evtz[j] = -100.;
00677 evtt[j] = -100.;
00678 }
00679 evtx[j] = ntpEvent->vtx.x;
00680 evty[j] = ntpEvent->vtx.y;
00681 evtz[j] = ntpEvent->vtx.z;
00682 evtt[j] = ntpEvent->vtx.t;
00683 }
00684
00685
00686 Double_t stripArrayTime[500][192][2];
00687 for(int ii=0;ii<500;ii++) {
00688 for(int jj=0;jj<192;jj++) {
00689 stripArrayTime[ii][jj][0] = stripArrayTime[ii][jj][1] = -1.0;
00690 }
00691 }
00692 for(int j=0;j<int(eventSummary->nstrip);j++){
00693 if(!LoadStrip(j)) continue;
00694 Int_t pl = ntpStrip->plane;
00695 Int_t st = ntpStrip->strip;
00696 if(ntpStrip->time1<stripArrayTime[pl][st][1] ||
00697 stripArrayTime[pl][st][1]<=0) {
00698 stripArrayTime[pl][st][1] = ntpStrip->time1;
00699 }
00700 if(ntpStrip->time0<stripArrayTime[pl][st][0] ||
00701 stripArrayTime[pl][st][0]<=0) {
00702 stripArrayTime[pl][st][0] = ntpStrip->time0;
00703 }
00704 }
00705
00706 if(nevent==0) continue;
00707
00708 for(int j=0;j<nevent;j++){
00709 if(!LoadEvent(j)) continue;
00710
00711 event = j;
00712 if(LoadSlice(ntpEvent->slc)) {
00713 slice = ntpSlice->index;
00714 slclength = ntpSlice->plane.n;
00715 }
00716 ntrack = ntpEvent->ntrack;
00717 nshower = ntpEvent->nshower;
00718
00719
00720 true_enu = 0; true_emu = 0; true_elep = 0; true_eshw = 0;
00721 true_pxnu = 0; true_pynu = 0; true_pznu = 0;
00722 true_x = 0; true_y = 0; true_q2 = 0; true_w2 = 0;
00723 true_dircosneu = 0; true_dircosz = 0; true_emfrac = 0;
00724 true_vtxx = 0; true_vtxy = 0; true_vtxz = 0;
00725 th_evt_pur = 0; th_evt_all = 0; th_evt_slc = 0;
00726 th_shw_pur = 0; th_shw_all = 0; th_shw_slc = 0;
00727 th_trk_pur = 0; th_trk_all = 0; th_trk_slc = 0;
00728 flavour = 0; nooscflavour = 0; process = 0; nucleus = 0; cc_nc = 0;
00729 initial_state = 0; had_fs = 0;
00730 true_etgt = 0; true_pxtgt = 0; true_pytgt = 0; true_pztgt = 0;
00731 nuparent->Zero();
00732 detector = -1; mcevent = -1;
00733 is_fid = is_cev = is_mc = pass_basic = pass = 0;
00734 pid0 = pid1 = pid2 = 0.;
00735 reco_enu = reco_emu = reco_eshw = 0;
00736 reco_eshw_linCC = reco_eshw_wtCC = reco_eshw_linNC = 0;
00737 reco_eshw_wtNC = reco_eshw_em = 0;
00738 reco_qe_enu = reco_qe_q2 = reco_y = reco_dircosneu = 0;
00739 reco_dircosz = reco_emfrac = reco_emfrac_prob = 0;
00740 reco_vtxx = reco_vtxy = reco_vtxz = 0;
00741 evtlength = trklength = shwlength = 0;
00742 closestevent_s = closestevent_t = -1.;
00743 fracRehitStrip = 0;
00744 trkphfrac = trkphplane = trkmomrange = trkqp = trkeqp = 0;
00745 shwstptimemin = shwstptimemax = shwstptimerms = 0.0;
00746 this->SetStdHepVars(-1);
00747
00748
00749 if(ntpHeader->GetVldContext().
00750 GetDetector()==Detector::kFar) detector = 2;
00751 else if(ntpHeader->GetVldContext().
00752 GetDetector()==Detector::kNear) detector = 1;
00753
00754
00755 is_fid = 1;
00756 if(!IsFidVtxEvt(ntpEvent->index)) is_fid = 0;
00757
00758
00759 for(int k=0; k<ntpEvent->nstrip; k++){
00760 if(!LoadStrip(ntpEvent->stp[k])) continue;
00761 Int_t pl = ntpStrip->plane;
00762 Int_t st = ntpStrip->strip;
00763 if(ntpStrip->time0>stripArrayTime[pl][st][0]) {
00764 fracRehitStrip+=1;
00765 }
00766 else if(ntpStrip->time1>stripArrayTime[pl][st][1]) {
00767 fracRehitStrip+=1;
00768 }
00769 }
00770 fracRehitStrip/=ntpEvent->nstrip;
00771
00772 Bool_t good_truth = false;
00773
00774 if(LoadTruthForRecoTH(j,mcevent)) {
00775 good_truth = true;
00776 Float_t *vtx = TrueVtx(mcevent);
00777 Float_t *nu3mom = TrueNuMom(mcevent);
00778 Float_t *targ4vec = Target4Vector(mcevent);
00779
00780 is_mc = 1;
00781 nucleus = Nucleus(mcevent);
00782 flavour = Flavour(mcevent);
00783 nooscflavour = ntpTruth->inunoosc;
00784 initial_state = Initial_state(mcevent);
00785 had_fs = HadronicFinalState(mcevent);
00786 process = IResonance(mcevent);
00787 cc_nc = IAction(mcevent);
00788 true_enu = TrueNuEnergy(mcevent);
00789 true_pxnu = nu3mom[0];
00790 true_pynu = nu3mom[1];
00791 true_pznu = nu3mom[2];
00792 true_emu = TrueMuEnergy(mcevent);
00793 true_elep = TrueLeptonEnergy(mcevent);
00794 true_eshw = TrueShwEnergy(mcevent);
00795 true_x = X(mcevent);
00796 true_y = Y(mcevent);
00797 true_q2 = Q2(mcevent);
00798 true_w2 = W2(mcevent);
00799 true_dircosz = TrueMuDCosZ(mcevent);
00800 true_dircosneu = TrueMuDCosNeu(mcevent);
00801 true_vtxx = vtx[0];
00802 true_vtxy = vtx[1];
00803 true_vtxz = vtx[2];
00804 true_etgt = targ4vec[3];
00805 true_pxtgt = targ4vec[0];
00806 true_pytgt = targ4vec[1];
00807 true_pztgt = targ4vec[2];
00808 true_emfrac = ntpTruth->emfrac;
00809
00810 if(LoadTHEvent(j)) {
00811 th_evt_pur = ntpTHEvent->purity;
00812 th_evt_all = ntpTHEvent->completeall;
00813 th_evt_slc = ntpTHEvent->completeslc;
00814 }
00815
00816 this->SetStdHepVars(mcevent);
00817
00818 if(gnumi->Status()) {
00819 if(isST) gnumi->GetParent(strecord,mcevent,*nuparent);
00820 else gnumi->GetParent(mcrecord,mcevent,*nuparent);
00821 }
00822
00823 delete [] vtx;
00824 delete [] nu3mom;
00825 delete [] targ4vec;
00826 }
00827
00828
00829 this->CallUserFunctions(event);
00830
00831
00832 reco_vtxx = ntpEvent->vtx.x;
00833 reco_vtxy = ntpEvent->vtx.y;
00834 reco_vtxz = ntpEvent->vtx.z;
00835 evtlength = ntpEvent->plane.n;
00836 if(ntpHeader->GetVldContext().GetDetector()==Detector::kNear){
00837 evtlength=ntpEvent->plane.end-ntpEvent->plane.beg+1;
00838 }
00839
00840 if(nevent>1) {
00841 for(int k=0;k<nevent;k++){
00842 if(k!=j){
00843 Double_t t_sep = TMath::Abs(evtt[k]-evtt[j]);
00844 if(closestevent_t<0 || t_sep<closestevent_t){
00845 closestevent_t = t_sep;
00846 closestevent_s = TMath::Sqrt(TMath::Power(evtx[k]-evtx[j],2)+
00847 TMath::Power(evty[k]-evty[j],2)+
00848 TMath::Power(evtz[k]-evtz[j],2));
00849
00850 }
00851 }
00852 }
00853 }
00854
00855
00856 int track_index = -1;
00857 int shower_index = -1;
00858 if(LoadLargestTrackFromEvent(j,track_index)) {
00859 if(!LoadShowerAtTrackVertex(j,track_index,shower_index)) {
00860 LoadLargestShowerFromEvent(j,shower_index);
00861 }
00862 }
00863 else LoadLargestShowerFromEvent(j,shower_index);
00864
00865 if(good_truth) {
00866 if(LoadTHShower(shower_index)) {
00867 th_shw_pur = ntpTHShower->purity;
00868 th_shw_all = ntpTHShower->completeall;
00869 th_shw_slc = ntpTHShower->completeslc;
00870 }
00871 if(LoadTHTrack(track_index)) {
00872 th_trk_pur = ntpTHTrack->purity;
00873 th_trk_all = ntpTHTrack->completeall;
00874 th_trk_slc = ntpTHTrack->completeslc;
00875 }
00876 }
00877
00878
00879 reco_emu = RecoMuEnergy(0,track_index);
00880 reco_eshw_linCC = RecoShwEnergy(shower_index,0);
00881 reco_eshw_wtCC = RecoShwEnergy(shower_index,1);
00882 reco_eshw_linNC = RecoShwEnergy(shower_index,2);
00883 reco_eshw_wtNC = RecoShwEnergy(shower_index,3);
00884 reco_eshw_em = RecoShwEnergy(shower_index,4);
00885 reco_eshw = reco_eshw_linCC;
00886 reco_enu = reco_emu + reco_eshw;
00887 reco_qe_enu = RecoQENuEnergy(track_index);
00888 reco_qe_q2 = RecoQEQ2(track_index);
00889 if (reco_enu>0) reco_y = reco_eshw/reco_enu;
00890 reco_dircosz = RecoMuDCosZ(track_index);
00891 reco_dircosneu = RecoMuDCosNeu(track_index);
00892 is_cev = IsFidAll(track_index);
00893
00894 Double_t emfrac0 = 0;
00895 Double_t emfrac1 = 0;
00896 reco_emfrac = RecoEMFrac(shower_index,0,emfrac0,emfrac1);
00897 reco_emfrac_prob = RecoEMFrac(shower_index,1,emfrac0,emfrac1);
00898
00899
00900 if(PassTrackCuts(track_index)){
00901 trklength = ntpTrack->plane.n;
00902 trkmomrange = ntpTrack->momentum.range;
00903 trkqp = ntpTrack->momentum.qp;
00904 trkeqp = ntpTrack->momentum.eqp;
00905 Float_t phtrack=ntpTrack->ph.sigcor;
00906 Float_t phevent=ntpEvent->ph.sigcor;
00907 if (phevent>0) {trkphfrac=phtrack/phevent;}
00908 if(ntpTrack->plane.n>0) {
00909 trkphplane = ntpTrack->ph.sigcor/ntpTrack->plane.n;
00910 }
00911 }
00912 else {
00913 trklength = 0;
00914 trkmomrange = 0;
00915 trkqp = 0;
00916 trkeqp = 0;
00917 trkphfrac = 0;
00918 trkphplane = 0;
00919 }
00920
00921 if(shower_index!=-1){
00922 shwlength = ntpShower->plane.n;
00923 Double_t sum_x = 0;
00924 Double_t sum_x2 = 0;
00925 for(int k=0;k<ntpShower->nstrip;k++){
00926 if(!LoadStrip(ntpShower->stp[k])) continue;
00927 Double_t t0 = ntpStrip->time0;
00928 Double_t t1 = ntpStrip->time1;
00929 Double_t avet = 0;
00930 if(t0>0&&t1>0) avet = (t0+t1)/2.;
00931 else if(t0>0) avet = t0;
00932 else if(t1>0) avet = t1;
00933 if(k==0) shwstptimemin = shwstptimemax = avet;
00934 else {
00935 if(avet>shwstptimemax) shwstptimemax = avet;
00936 if(avet<shwstptimemin) shwstptimemin = avet;
00937 }
00938 sum_x += (ntpStrip->ph0.pe + ntpStrip->ph1.pe)*avet;
00939 sum_x2 += (ntpStrip->ph0.pe + ntpStrip->ph1.pe)*avet*avet;
00940 }
00941 shwstptimemin -= eventSummary->trigtime;
00942 shwstptimemax -= eventSummary->trigtime;
00943 if(ntpShower->ph.pe>0){
00944 sum_x2 /= ntpShower->ph.pe;
00945 sum_x /= ntpShower->ph.pe;
00946 shwstptimerms = sum_x2 - sum_x*sum_x;
00947 if(shwstptimerms>0) shwstptimerms = TMath::Sqrt(shwstptimerms);
00948 else shwstptimerms = 0;
00949 }
00950 }
00951
00952 if(PassBasicCuts(event)) {
00953 pass_basic = 1;
00954 pid0 = PID(event,0);
00955 pid1 = PID(event,1);
00956 pid2 = PID(event,2);
00957 if(PassAnalysisCuts(event)) pass = 1;
00958 }
00959
00960
00961 tree->Fill();
00962 }
00963 }
00964 delete nuparent;
00965
00966 save.cd();
00967 pottree->Fill();
00968 pottree->Write();
00969 tree->Write();
00970
00971 save.Close();
00972
00973 }
00974
00975
00976 void MadAnalysis::CreateANtpPAN(std::string tag){
00977
00978 std::string savename = "PAN_" + tag + ".root";
00979 TFile save(savename.c_str(),"RECREATE");
00980 save.cd();
00981
00982 GnumiInterface gnumi;
00983 if(!gnumi.Status()) {
00984 cout << "Error MadAnalysis::CreatePAN() - Flux files not open."
00985 << " Will not be filling NuParent info"
00986 << endl;
00987 cout << "Set environmental variable $GNUMIAUX to point to the "
00988 << "directory containing the gnumi files to fix this!"
00989 << endl;
00990
00991 }
00992
00993
00994 ANtpHeaderInfo *antpHeader = new ANtpHeaderInfo();
00995 ANtpTruthInfoBeam *antpTruth = new ANtpTruthInfoBeam();
00996 ANtpAnalysisInfo *antpAnalysis = new ANtpAnalysisInfo();
00997 ANtpRecoInfo *antpReco = new ANtpRecoInfo();
00998 ANtpEventInfo *antpEvent = new ANtpEventInfo();
00999 ANtpTrackInfo *antpTrack = new ANtpTrackInfo();
01000 ANtpShowerInfo *antpShower = new ANtpShowerInfo();
01001 ANtpRecoNtpManipulator *ntpManip = new ANtpRecoNtpManipulator();
01002 ANtpInfoObjectFillerBeam filla;
01003
01004
01005 Int_t mcevent;
01006 Int_t is_mc;
01007
01008
01009 TTree *tree = new TTree("pan","pan");
01010 tree->Branch("header","ANtpHeaderInfo",&antpHeader,8000,1);
01011 if(isMC) {
01012 tree->Branch("is_mc",&is_mc,"is_mc/I",32000);
01013 tree->Branch("mcevent",&mcevent,"mcevent/I",32000);
01014 tree->Branch("truth","ANtpTruthInfoBeam",&antpTruth,8000,1);
01015 }
01016 tree->Branch("reco","ANtpRecoInfo",&antpReco,8000,1);
01017 tree->Branch("analysis","ANtpAnalysisInfo",&antpAnalysis,8000,1);
01018 tree->Branch("event","ANtpEventInfo",&antpEvent,8000,1);
01019 tree->Branch("track","ANtpTrackInfo",&antpTrack,8000,1);
01020 tree->Branch("shower","ANtpShowerInfo",&antpShower,8000,1);
01021
01022 for(int i=0;i<Nentries;i++){
01023 if(i%1000==0) std::cout << float(i)*100./float(Nentries)
01024 << "% done" << std::endl;
01025
01026 if(!this->GetEntry(i)) continue;
01027 if(eventSummary->nevent==0) continue;
01028
01029
01030 for(int j=0;j<eventSummary->nevent;j++){
01031
01032 if(!LoadEvent(j)) continue;
01033
01034
01035 antpHeader->Reset();
01036 is_mc = 0;
01037 mcevent = -1;
01038 antpTruth->Reset();
01039 antpAnalysis->Reset();
01040 antpReco->Reset();
01041 antpEvent->Reset();
01042 antpTrack->Reset();
01043 antpShower->Reset();
01044 ntpManip->SetRecord(strecord);
01045
01046 antpHeader->events = eventSummary->nevent;
01047 antpHeader->run = ntpHeader->GetRun();
01048 antpHeader->subRun = ntpHeader->GetSubRun();
01049 antpHeader->snarl = ntpHeader->GetSnarl();
01050
01051 UInt_t year,month,day,hour,minute,second;
01052 ntpHeader->GetVldContext().GetTimeStamp().GetDate(kFALSE,0,&year,
01053 &month,&day);
01054 ntpHeader->GetVldContext().GetTimeStamp().GetTime(kFALSE,0,&hour,
01055 &minute,&second);
01056 antpHeader->year = year;
01057 antpHeader->month = month;
01058 antpHeader->day = day;
01059 antpHeader->hour = hour;
01060 antpHeader->minute = minute;
01061 antpHeader->second = second;
01062 antpHeader->utc = ntpHeader->GetVldContext().GetTimeStamp().GetSec();
01063
01064 antpEvent->index = j;
01065 filla.FillEventInformation(ntpManip,ntpEvent,antpEvent);
01066
01067
01068 antpHeader->detector = ntpHeader->GetVldContext().GetDetector();
01069
01070
01071 antpReco->inFiducialVolume = 1;
01072 if(!IsFidVtxEvt(ntpEvent->index))
01073 antpReco->inFiducialVolume = 0;
01074
01075
01076 if(LoadTruthForRecoTH(antpEvent->index,mcevent)) {
01077
01078 is_mc = 1;
01079 filla.FillMCTruthInformation(ntpTruth,antpTruth);
01080 if(isST) filla.FillBeamMCTruthInformation(ntpTruth,strecord->stdhep,antpTruth);
01081 else filla.FillBeamMCTruthInformation(ntpTruth,mcrecord->stdhep,antpTruth);
01082 if(gnumi.Status()) {
01083 if(isST) gnumi.FillANtpTruth(strecord,mcevent,*antpTruth);
01084 else gnumi.FillANtpTruth(mcrecord,mcevent,*antpTruth);
01085 }
01086 }
01087
01088 if(PassBasicCuts(antpEvent->index)) {
01089
01090 int track_index = -1;
01091 LoadLargestTrackFromEvent(antpEvent->index,track_index);
01092 if(track_index!=-1) filla.FillTrackInformation(ntpTrack,antpTrack);
01093
01094 int shower_index = -1;
01095 LoadLargestShowerFromEvent(antpEvent->index,shower_index);
01096 if(shower_index!=-1) filla.FillShowerInformation(ntpShower,antpShower);
01097
01098 antpReco->passesCuts = 1;
01099 antpReco->eventLength = TMath::Abs(antpEvent->endPlane -
01100 antpEvent->begPlane);
01101 antpReco->trackLength = TMath::Abs(antpTrack->endPlane -
01102 antpTrack->begPlane);
01103 antpReco->trackMomentum = antpTrack->fitMomentum;
01104 antpReco->trackRange = antpTrack->rangeMomentum;
01105 antpReco->sigmaQoverP = antpTrack->sigmaQoverP;
01106
01107 antpReco->muEnergy = RecoMuEnergy(0,track_index);
01108 antpReco->showerEnergy = RecoShwEnergy(shower_index);
01109 antpReco->nuEnergy = ( antpReco->muEnergy +
01110 antpReco->showerEnergy );
01111
01112 antpReco->qeNuEnergy = RecoQENuEnergy(track_index);
01113 antpReco->qeQ2 = RecoQEQ2(track_index);
01114 if (antpReco->nuEnergy>0) {
01115 antpReco->hadronicY = ( antpReco->showerEnergy /
01116 antpReco->nuEnergy );
01117 }
01118
01119 antpReco->muDCosZVtx = RecoMuDCosZ(track_index);
01120 antpReco->nuDCos = RecoMuDCosNeu(track_index);
01121 antpReco->vtxX = antpEvent->vtxX;
01122 antpReco->vtxY = antpEvent->vtxY;
01123 antpReco->vtxZ = antpEvent->vtxZ;
01124 antpReco->isFullyContained = IsFidAll(track_index);
01125
01126 antpAnalysis->separationParameter = PID(antpEvent->index,0);
01127 if(PassAnalysisCuts(antpEvent->index)) antpReco->pass = 1;
01128 else antpReco->pass = 0;
01129 }
01130
01131 tree->Fill();
01132 }
01133 }
01134
01135 delete antpHeader;
01136 delete antpTruth;
01137 delete antpAnalysis;
01138 delete antpReco;
01139 delete antpEvent;
01140 delete antpTrack;
01141 delete antpShower;
01142
01143 save.cd();
01144 tree->Write();
01145 save.Write();
01146 save.Close();
01147 }
01148
01149
01150 void MadAnalysis::RecoExperiment(int startnum,
01151 double NearExpectedNormToPOT,
01152 double FarExpectedNormToPOT)
01153 {
01154
01155 cout << "Reconstructing Data Energy Spectrum" << endl;
01156 if(!DataHist[0]){
01157 DataHist[0] =
01158 new TH1F("NearDataHist",
01159 "Near Detector Reconstructed Energy Spectrum",
01160 NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01161 DataHist[0]->SetXTitle("E_{reco} (GeV)");
01162 DataHist[0]->SetYTitle("Event Rate");
01163 DataHist[0]->Sumw2();
01164 }
01165 else {
01166 DataHist[0]->Reset();
01167 numDataEvts[0] = 0;
01168 signalFracInData[0] = 0;
01169 }
01170
01171 if(!DataHist[1]){
01172 DataHist[1] =
01173 new TH1F("FarDataHist",
01174 "Far Detector Oscillated Reconstructed Energy Spectrum",
01175 NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01176 DataHist[1]->SetXTitle("E_{reco} (GeV)");
01177 DataHist[1]->SetYTitle("Event Rate");
01178 DataHist[1]->Sumw2();
01179 }
01180 else {
01181 DataHist[1]->Reset();
01182 numDataEvts[1] = 0;
01183 signalFracInData[1] = 0;
01184 }
01185
01186 Float_t midupbin = DataHist[0]->GetBinCenter(NumberOfBins);
01187
01188 int i = startnum;
01189 while(i<=Nentries){
01190 if(!GetEntry(i++)) continue;
01191
01192 if((i-startnum)%500==0)
01193 std::cout << 100*float(i-startnum)/float(Nentries-startnum)
01194 << "% completed" << std::endl;
01195
01196
01197 int theDetector = ntpHeader->GetVldContext().GetDetector()-1;
01198
01199
01200 for(int j=0;j<eventSummary->nevent;j++){
01201 numDataEvts[theDetector]+=1;
01202 if(PassAnalysisCuts(j)){
01203 signalFracInData[theDetector]+=1;
01204 float ereco = RecoAnalysisEnergy(j);
01205 if(ereco<RangeUpperLimit)
01206 DataHist[theDetector]->Fill(ereco,GetWeight());
01207 else DataHist[theDetector]->Fill(midupbin,GetWeight());
01208 }
01209 }
01210 }
01211
01212 signalFracInData[0]/=float(numDataEvts[0]);
01213 signalFracInData[1]/=float(numDataEvts[1]);
01214 DataPOT[0]=numDataEvts[0]*NearExpectedNormToPOT;
01215 DataPOT[1]=numDataEvts[1]*FarExpectedNormToPOT;
01216
01217 }
01218
01219
01220 void MadAnalysis::RecoMCExperiment(int startnum,double NearPOT,double FarPOT)
01221 {
01222
01223 cout << "Reconstructing MC Data Energy Spectrum" << endl;
01224
01225 if(!DataHist[0]){
01226 DataHist[0] =
01227 new TH1F("NearDataHist",
01228 "Near Detector Reconstructed Energy Spectrum",
01229 NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01230 DataHist[0]->SetXTitle("E_{reco} (GeV)");
01231 DataHist[0]->SetYTitle("Event Rate");
01232 DataHist[0]->Sumw2();
01233
01234 DataSignalHist[0] =
01235 new TH1F("NearDataSignalHist",
01236 "Near Detector Reconstructed Energy Spectrum (Signal only)",
01237 NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01238 DataSignalHist[0]->SetXTitle("E_{reco} (GeV)");
01239 DataSignalHist[0]->SetYTitle("Event Rate");
01240 DataSignalHist[0]->Sumw2();
01241
01242 DataBkgdHist[0] =
01243 new TH1F("NearDataBkgdHist",
01244 "Near Detector Reconstructed Energy Spectrum (Background only)",
01245 NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01246 DataBkgdHist[0]->SetXTitle("E_{reco} (GeV)");
01247 DataBkgdHist[0]->SetYTitle("Event Rate");
01248 DataBkgdHist[0]->Sumw2();
01249 }
01250 else {
01251 DataHist[0]->Reset();
01252 DataSignalHist[0]->Reset();
01253 DataBkgdHist[0]->Reset();
01254 numDataEvts[0] = 0;
01255 signalFracInData[0] = 0;
01256 }
01257
01258 if(!DataHist[1]){
01259 DataHist[1] =
01260 new TH1F("FarDataHist",
01261 "Far Detector Oscillated Reconstructed Energy Spectrum",
01262 NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01263 DataHist[1]->SetXTitle("E_{reco} (GeV)");
01264 DataHist[1]->SetYTitle("Event Rate");
01265 DataHist[1]->Sumw2();
01266
01267 DataSignalHist[1] = new TH1F("FarDataSignalHist",
01268 "Far Detector Oscillated Reconstructed Energy Spectrum (Signal only)",
01269 NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01270 DataSignalHist[1]->SetXTitle("E_{reco} (GeV)");
01271 DataSignalHist[1]->SetYTitle("Event Rate");
01272 DataSignalHist[1]->Sumw2();
01273
01274 DataBkgdHist[1] = new TH1F("FarDataBkgdHist",
01275 "Far Detector Oscillated Reconstructed Energy Spectrum (Background only)",
01276 NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01277 DataBkgdHist[1]->SetXTitle("E_{reco} (GeV)");
01278 DataBkgdHist[1]->SetYTitle("Event Rate");
01279 DataBkgdHist[1]->Sumw2();
01280 }
01281 else {
01282 DataHist[1]->Reset();
01283 DataSignalHist[1]->Reset();
01284 DataBkgdHist[1]->Reset();
01285 numDataEvts[1] = 0;
01286 signalFracInData[1] = 0;
01287 }
01288
01289
01290 Float_t midupbin = DataHist[0]->GetBinCenter(NumberOfBins);
01291
01292 DataPOT[0] = NearPOT;
01293 DataPOT[1] = FarPOT;
01294 numDataEvts[0] = int(NEARPOTTONUMEVT*NearPOT);
01295 numDataEvts[1] = int(FARPOTTONUMEVT*FarPOT);
01296
01297 int counter[2] = {0};
01298
01299 int i = startnum;
01300 while((counter[0]<numDataEvts[0]||counter[1]<numDataEvts[1])&&i<=Nentries){
01301 if(!GetEntry(i++)) continue;
01302
01303
01304 int theDetector = ntpHeader->GetVldContext().GetDetector()-1;
01305
01306
01307 for(int j=0;j<eventSummary->nevent;j++){
01308
01309 int mcevent = -1;
01310 if(LoadTruthForRecoTH(j,mcevent)){
01311
01312 if(PassTruthSignal(mcevent)) counter[theDetector]+=1;
01313
01314 if(counter[theDetector]%1000==0)
01315 std::cout << 100*(float(counter[theDetector])
01316 /float(numDataEvts[theDetector]))
01317 << "% completed for Detector " << theDetector
01318 << std::endl;
01319
01320 if(PassAnalysisCuts(j)){
01321 float ereco = RecoAnalysisEnergy(j);
01322 if(PassTruthSignal(mcevent)) {
01323 signalFracInData[theDetector]+=1;
01324
01325 double oscprob = 0;
01326 if(theDetector==1)
01327 OscillationFunction->Eval(TrueNuEnergy(mcevent));
01328
01329 if(ereco<RangeUpperLimit) {
01330 DataHist[theDetector]->Fill(ereco,(1.-oscprob)*GetWeight());
01331 DataSignalHist[theDetector]->Fill(ereco,
01332 (1.-oscprob)*GetWeight());
01333 }
01334 else {
01335 DataHist[theDetector]->Fill(midupbin,(1.-oscprob)*GetWeight());
01336 DataSignalHist[theDetector]->Fill(midupbin,
01337 (1.-oscprob)*GetWeight());
01338 }
01339 }
01340 else {
01341 if(ereco<RangeUpperLimit) {
01342 DataHist[theDetector]->Fill(ereco,GetWeight());
01343 DataBkgdHist[theDetector]->Fill(ereco,GetWeight());
01344 }
01345 else {
01346 DataHist[theDetector]->Fill(midupbin,GetWeight());
01347 DataBkgdHist[theDetector]->Fill(midupbin,GetWeight());
01348 }
01349 }
01350 }
01351 }
01352 }
01353 }
01354
01355 if(counter[0]<numDataEvts[0]) numDataEvts[0]=counter[0];
01356 DataPOT[0]=numDataEvts[0]/NEARPOTTONUMEVT;
01357 signalFracInData[0]/=float(numDataEvts[0]);
01358 if(counter[1]<numDataEvts[1]) numDataEvts[1]=counter[1];
01359 DataPOT[1]=numDataEvts[1]/FARPOTTONUMEVT;
01360 signalFracInData[1]/=float(numDataEvts[1]);
01361
01362 }
01363
01364
01365 void MadAnalysis::RecoMC(int startnum,double NearPOT,double FarPOT)
01366 {
01367 cout << "Reconstructing MC Event Sample" << endl;
01368
01369 if(!MCSignalHist[0]){
01370 MCSignalHist[0] = new TH1F("NearMCSignalHist",
01371 "Near Detector Reconstructed Unoscillated Energy Spectrum (Signal only)",
01372 NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01373 MCSignalHist[0]->SetXTitle("E_{reco} (GeV)");
01374 MCSignalHist[0]->SetYTitle("Event Rate");
01375 MCSignalHist[0]->Sumw2();
01376
01377 MCBkgdHist[0] = new TH1F("NearMCBkgdHist",
01378 "Near Detector Reconstructed Unoscillated Energy Spectrum (Bkd only)",
01379 NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01380 MCBkgdHist[0]->SetXTitle("E_{reco} (GeV)");
01381 MCBkgdHist[0]->SetYTitle("Event Rate");
01382 MCBkgdHist[0]->Sumw2();
01383 }
01384 else {
01385 MCSignalHist[0]->Reset();
01386 MCBkgdHist[0]->Reset();
01387 MCPOT[0] = 0;
01388 numMCEvts[0] = 0;
01389 }
01390
01391 if(!MCSignalHist[1]){
01392 MCSignalHist[1] = new TH1F("FarMCSignalHist",
01393 "Far Detector Reconstructed Unoscillated Energy Spectrum (Signal only)",
01394 NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01395 MCSignalHist[1]->SetXTitle("E_{reco} (GeV)");
01396 MCSignalHist[1]->SetYTitle("Event Rate");
01397 MCSignalHist[1]->Sumw2();
01398
01399 MCBkgdHist[1] = new TH1F("FarMCBkgdHist",
01400 "Far Detector Reconstructed Unoscillated Energy Spectrum (Bkd only)",
01401 NumberOfBins,RangeLowerLimit,RangeUpperLimit);
01402 MCBkgdHist[1]->SetXTitle("E_{reco} (GeV)");
01403 MCBkgdHist[1]->SetYTitle("Event Rate");
01404 MCBkgdHist[1]->Sumw2();
01405 }
01406 else {
01407 MCSignalHist[1]->Reset();
01408 MCBkgdHist[1]->Reset();
01409 MCPOT[1] = 0;
01410 numMCEvts[1] = 0;
01411 }
01412
01413
01414 Float_t midupbin = MCSignalHist[0]->GetBinCenter(NumberOfBins);
01415
01416 MCPOT[0] = NearPOT; MCPOT[1] = FarPOT;
01417 numMCEvts[0] = int(NEARPOTTONUMEVT*NearPOT);
01418 numMCEvts[1] = int(FARPOTTONUMEVT*FarPOT);
01419 int counter[2] = {0};
01420 int i = startnum;
01421
01422 while((counter[0]<numMCEvts[0]||counter[1]<numMCEvts[1])&&i<=Nentries){
01423 if(!GetEntry(i++)) continue;
01424
01425
01426 int theDetector = ntpHeader->GetVldContext().GetDetector()-1;
01427
01428
01429 for(int j=0;j<eventSummary->nevent;j++){
01430
01431
01432 int mcevent = -1;
01433 if(LoadTruthForRecoTH(j,mcevent)){
01434
01435 if(PassTruthSignal(mcevent)) counter[theDetector]+=1;
01436
01437 if(counter[theDetector]%1000==0)
01438 std::cout << 100*(float(counter[theDetector])
01439 /float(numDataEvts[theDetector]))
01440 << "% completed for Detector " << theDetector
01441 << std::endl;
01442
01443 if(PassAnalysisCuts(j)){
01444 float emu = 0;
01445 float eshw = 0;
01446 float mudcosz = 0;
01447 float ereco = RecoAnalysisEnergy(j);
01448 if(PassTruthSignal(mcevent)) {
01449 if(ereco<RangeUpperLimit)
01450 MCSignalHist[theDetector]->Fill(ereco,GetWeight());
01451 else MCSignalHist[theDetector]->Fill(midupbin,GetWeight());
01452 }
01453 else {
01454 if(ereco<RangeUpperLimit)
01455 MCBkgdHist[theDetector]->Fill(ereco,GetWeight());
01456 else MCBkgdHist[theDetector]->Fill(midupbin,GetWeight());
01457 }
01458
01459 mcev_reweight mmei = {TrueNuEnergy(mcevent),TrueMuEnergy(mcevent),
01460 TrueShwEnergy(mcevent),TrueMuDCosZ(mcevent),
01461 ereco,emu,eshw,mudcosz,GetWeight(),
01462 PassTruthSignal(mcevent),INu(mcevent),
01463 IAction(mcevent),IResonance(mcevent),
01464 Initial_state(mcevent),Nucleus(mcevent),
01465 X(mcevent),Y(mcevent),Q2(mcevent),W2(mcevent)};
01466 MCEvents[theDetector].push_back(mmei);
01467 }
01468 }
01469 }
01470 }
01471
01472 if(counter[0]<numMCEvts[0]) numMCEvts[0]=counter[0];
01473 MCPOT[0]=numMCEvts[0]/NEARPOTTONUMEVT;
01474 if(counter[1]<numMCEvts[1]) numMCEvts[1]=counter[1];
01475 MCPOT[1]=numMCEvts[1]/FARPOTTONUMEVT;
01476
01477 }
01478
01479
01480 void MadAnalysis::SetOscillationFunction(Double_t (*fcn)(Double_t*,Double_t*),
01481 Float_t lowend,Float_t highend,
01482 int numpars,double *params)
01483 {
01484
01485 OscillationParameters = params;
01486 OscillationFunction = new TF1("OscFunc",fcn,lowend,
01487 highend,numpars);
01488 NumOscPars = numpars;
01489
01490 for(int i=0;i<NumOscPars;i++){
01491 OscillationFunction->SetParameter(i,OscillationParameters[i]);
01492 }
01493
01494 }
01495
01496
01497 void MadAnalysis::SetOscillationFunction(TF1 *form,double *params)
01498 {
01499
01500 OscillationParameters = params;
01501 OscillationFunction = new TF1(*form);
01502 NumOscPars = OscillationFunction->GetNpar();
01503
01504 for(int i=0;i<NumOscPars;i++){
01505 OscillationFunction->SetParameter(i,OscillationParameters[i]);
01506 }
01507
01508 }
01509
01510
01511 void MadAnalysis::SetNeugenReweightInfo(Int_t n,Int_t *bins, Float_t *min,
01512 Float_t *max, Float_t *mean,
01513 Float_t *err, std::string *name)
01514 {
01515 if(n>3) n = 3;
01516 for(int i=0;i<n;i++){
01517 NgenBins[i] = bins[i]; NgenMin[i] = min[i];
01518 NgenMax[i] = max[i]; NgenMean[i] = mean[i];
01519 NgenErr[i] = err[i]; NgenName[i] = name[i];
01520 }
01521 }
01522
01523
01524 Bool_t MadAnalysis::Do2DFit(float n_a,float min_a,float max_a,
01525 float n_b,float min_b,float max_b,
01526 float n_c,float f_c)
01527 {
01528
01529
01530
01531 if(!OscillationFunction||!DataHist[1]||!MCSignalHist[1]) {
01532 cout << "**ERROR** Did not find one of the following: " << endl;
01533 cout << "Oscillation Function,"
01534 << " Far Detector Reconstructed Data Histogram,"
01535 << " Far Detector Reconstructed MC Histogram" << endl;
01536 cout << "Ensure these are set and try again" << endl;
01537 return false;
01538 }
01539
01540
01541 Float_t midupbin = DataHist[1]->GetBinCenter(NumberOfBins);
01542
01543
01544 Float_t Par1Min=min_a;
01545 Float_t Par1Max=max_a;
01546 Int_t Par1Bins=int(n_a);
01547 Float_t Par1Width=(Par1Max-Par1Min)/Float_t(Par1Bins);
01548 Float_t Par2Min=min_b;
01549 Float_t Par2Max=max_b;
01550 Int_t Par2Bins=int(n_b);
01551 Float_t Par2Width=(Par2Max-Par2Min)/Float_t(Par2Bins);
01552
01553 Int_t EShiftBins=int(n_c);
01554 Float_t fracEShift=f_c;
01555 Float_t EShiftWidth=2.*f_c/n_c;
01556
01557 NDOF = NumberOfBins - 2;
01558
01559 Chi2Surf = new TH2F("Chi2Surf",
01560 "ChiSquare Surface",
01561 Par1Bins,Par1Min,Par1Max,
01562 Par2Bins,Par2Min,Par2Max);
01563
01564 cout << "Starting fit:" << endl;
01565
01566 for(int i=0;i<Par1Bins;i++){
01567 float par1 = Par1Min + (float(i+1)*Par1Width) - Par1Width/2.;
01568 if(i%10==0) cout << "Doing Par1 = " << par1 << endl;
01569 for(int j=0;j<Par2Bins;j++){
01570 float par2 = Par2Min + (float(j+1)*Par2Width) - Par2Width/2.;
01571 OscillationFunction->SetParameter(varyX[0],par1);
01572 OscillationFunction->SetParameter(varyX[1],par2);
01573
01574 float ChiSq = 999999;
01575 TH1F *TempBestFit = NULL;
01576
01577
01578 for(int k=0;k<EShiftBins;k++){
01579 float EShift = 1. - fracEShift + k*EShiftWidth;
01580
01581 TH1F *temp_sig = new TH1F("temp_sig","temp_sig",
01582 NumberOfBins,RangeLowerLimit,
01583 RangeUpperLimit);
01584 TH1F *temp_bkg = new TH1F("temp_bkg","temp_bkg",
01585 NumberOfBins,RangeLowerLimit,
01586 RangeUpperLimit);
01587
01588 vector<mcev_reweight>::iterator mcbeg = MCEvents[1].begin();
01589 vector<mcev_reweight>::iterator mcend = MCEvents[1].end();
01590 while(mcbeg!=mcend){
01591 Float_t Energy = mcbeg->TrueNuEnergy;
01592 Float_t Ereco = mcbeg->RecoNuEnergy*EShift;
01593 Float_t weight = mcbeg->Weight;
01594 if(mcbeg->PassTruth){
01595 Double_t oscprob = OscillationFunction->Eval(Energy);
01596 if(oscprob>1.0) oscprob=1.0;
01597 if(Ereco>RangeUpperLimit) Ereco = midupbin;
01598 temp_sig->Fill(Ereco,(1.-oscprob)*weight);
01599 }
01600 else {
01601 if(Ereco>RangeUpperLimit) Ereco = midupbin;
01602 temp_bkg->Fill(Ereco,weight);
01603 }
01604 mcbeg++;
01605 }
01606
01607
01608 float chisq = Chi2Calc->GetChi2(DataHist[1],temp_sig,temp_bkg,
01609 MCPOT[1]/DataPOT[1]);
01610
01611 if(EShiftErr>0) chisq += TMath::Power((1.-EShift)/EShiftErr,2);
01612 if(chisq<ChiSq) {
01613 delete TempBestFit;
01614 TempBestFit = temp_sig;
01615 TempBestFit->SetName("TempFarBestFit");
01616 ChiSq = chisq;
01617 }
01618 else delete temp_sig;
01619 delete temp_bkg;
01620 }
01621
01622 if(ChiSq<ChiMin){
01623 delete BestFit[1];
01624 ChiMin=ChiSq;
01625 Par1MinVal=par1;
01626 Par2MinVal=par2;
01627 TempBestFit->Scale(DataPOT[1]/MCPOT[1]);
01628 TempBestFit->SetName("FarBestFit");
01629 BestFit[1] = TempBestFit;
01630 }
01631 else delete TempBestFit;
01632 Chi2Surf->Fill(par1,par2,ChiSq);
01633 }
01634 }
01635 std::cout << "Minimum value of Chi2 is " << ChiMin << std::endl;
01636 std::cout << "at Par1 = " << Par1MinVal << " and Par2 = "
01637 << Par2MinVal << std::endl;
01638
01639 return true;
01640 }
01641
01642
01643 Bool_t MadAnalysis::Do2DFitNearFar(float n_a,float min_a,float max_a,
01644 float n_b,float min_b,float max_b,
01645 float n_c,float f_c)
01646 {
01647
01648
01649 MCReweight &fReweight = MCReweight::Instance();
01650
01651
01652
01653
01654 if(!OscillationFunction||!DataHist[0]||!MCSignalHist[0]) {
01655 cout << "**ERROR** Did not find one of the following: " << endl;
01656 cout << "Oscillation Function,"
01657 << " Reconstructed Data Histograms,"
01658 << " Reconstructed MC Histograms" << endl;
01659 cout << "Ensure these are set and try again" << endl;
01660 return false;
01661 }
01662
01663
01664 Float_t midupbin = DataHist[0]->GetBinCenter(NumberOfBins);
01665
01666
01667 Float_t Par1Min=min_a;
01668 Float_t Par1Max=max_a;
01669 Int_t Par1Bins=int(n_a);
01670 Float_t Par1Width=(Par1Max-Par1Min)/Float_t(Par1Bins);
01671 Float_t Par2Min=min_b;
01672 Float_t Par2Max=max_b;
01673 Int_t Par2Bins=int(n_b);
01674 Float_t Par2Width=(Par2Max-Par2Min)/Float_t(Par2Bins);
01675
01676 Int_t EShiftBins=int(n_c);
01677 Float_t fracEShift=f_c;
01678 Float_t EShiftWidth=2.*f_c/n_c;
01679
01680 NDOF = 2*(NumberOfBins -1) - 2;
01681
01682 Chi2Surf = new TH2F("Chi2Surf",
01683 "ChiSquare Surface",
01684 Par1Bins,Par1Min,Par1Max,
01685 Par2Bins,Par2Min,Par2Max);
01686
01687 cout << "Starting fit:" << endl;
01688
01689 for(int i=0;i<Par1Bins;i++){
01690 float par1 = Par1Min + (float(i+1)*Par1Width) - Par1Width/2.;
01691 if(i%10==0) cout << "Doing Par1 = " << par1 << endl;
01692 for(int j=0;j<Par2Bins;j++){
01693 float par2 = Par2Min + (float(j+1)*Par2Width) - Par2Width/2.;
01694 OscillationFunction->SetParameter(varyX[0],par1);
01695 OscillationFunction->SetParameter(varyX[1],par2);
01696
01697 float ChiSq = 999999;
01698 TH1F *TempBestFit[2];
01699 TempBestFit[0] = NULL; TempBestFit[1] = NULL;
01700
01701
01702 for(int k=0;k<EShiftBins;k++){
01703 float EShift = 1. - fracEShift + k*EShiftWidth;
01704
01705
01706
01707 for(int l=0;l<=NgenBins[0];l++){
01708 double NgenPar0 = NgenMin[0] + l*(NgenMax[0] -
01709 NgenMin[0])/NgenBins[0];
01710 for(int m=0;m<=NgenBins[1];m++){
01711 double NgenPar1 = NgenMin[1] + m*(NgenMax[1] -
01712 NgenMin[1])/NgenBins[1];
01713 for(int n=0;n<=NgenBins[2];n++){
01714 double NgenPar2 = NgenMin[2] + n*(NgenMax[2] -
01715 NgenMin[2])/NgenBins[2];
01716
01717 Registry rwtconfig;
01718 rwtconfig.Set(NgenName[0].data(),NgenPar0);
01719 rwtconfig.Set(NgenName[1].data(),NgenPar1);
01720 rwtconfig.Set(NgenName[2].data(),NgenPar2);
01721
01722 TH1F *temp_sig_ND = new TH1F("temp_sig_ND","temp_sig_ND",
01723 NumberOfBins,RangeLowerLimit,
01724 RangeUpperLimit);
01725 TH1F *temp_bkg_ND = new TH1F("temp_bkg_ND","temp_bkg_ND",
01726 NumberOfBins,RangeLowerLimit,
01727 RangeUpperLimit);
01728
01729 vector<mcev_reweight>::iterator mcbeg = MCEvents[0].begin();
01730 vector<mcev_reweight>::iterator mcend = MCEvents[0].end();
01731 while(mcbeg!=mcend){
01732 Float_t Ereco = mcbeg->RecoNuEnergy*EShift;
01733 Float_t weight = mcbeg->Weight;
01734 Registry event;
01735 event.Set("event:energy",mcbeg->TrueNuEnergy);
01736 event.Set("event:iaction",mcbeg->IAction);
01737 event.Set("event:inu",mcbeg->INu);
01738 event.Set("event:process",mcbeg->Process);
01739 event.Set("event:initial_state",mcbeg->InitState);
01740 event.Set("event:nucleus",mcbeg->Nucleus);
01741 event.Set("event:q2",mcbeg->KinQ2);
01742 Float_t mcreweight = fReweight.ComputeWeight(&event,
01743 &rwtconfig);
01744
01745 if(mcbeg->PassTruth){
01746 if(Ereco>RangeUpperLimit) Ereco = midupbin;
01747 temp_sig_ND->Fill(Ereco,weight*mcreweight);
01748 }
01749 else {
01750 if(Ereco>RangeUpperLimit) Ereco = midupbin;
01751 temp_bkg_ND->Fill(Ereco,weight*mcreweight);
01752 }
01753 mcbeg++;
01754 }
01755
01756 TH1F *temp_sig_FD = new TH1F("temp_sig_FD","temp_sig_FD",
01757 NumberOfBins,RangeLowerLimit,
01758 RangeUpperLimit);
01759 TH1F *temp_bkg_FD = new TH1F("temp_bkg_FD","temp_bkg_FD",
01760 NumberOfBins,RangeLowerLimit,
01761 RangeUpperLimit);
01762
01763 mcbeg = MCEvents[1].begin();
01764 mcend = MCEvents[1].end();
01765 while(mcbeg!=mcend){
01766 Float_t Energy = mcbeg->TrueNuEnergy;
01767 Float_t Ereco = mcbeg->RecoNuEnergy*EShift;
01768 Float_t weight = mcbeg->Weight;
01769 Registry event;
01770 event.Set("event:energy",mcbeg->TrueNuEnergy);
01771 event.Set("event:iaction",mcbeg->IAction);
01772 event.Set("event:inu",mcbeg->INu);
01773 event.Set("event:process",mcbeg->Process);
01774 event.Set("event:initial_state",mcbeg->InitState);
01775 event.Set("event:nucleus",mcbeg->Nucleus);
01776 event.Set("event:q2",mcbeg->KinQ2);
01777 Float_t mcreweight = fReweight.ComputeWeight(&event,
01778 &rwtconfig);
01779 if(mcbeg->PassTruth){
01780 Double_t oscprob = OscillationFunction->Eval(Energy);
01781 if(oscprob>1.0) oscprob=1.0;
01782 if(Ereco>RangeUpperLimit) Ereco = midupbin;
01783 temp_sig_FD->Fill(Ereco,(1.-oscprob)*weight*mcreweight);
01784 }
01785 else {
01786 if(Ereco>RangeUpperLimit) Ereco = midupbin;
01787 temp_bkg_FD->Fill(Ereco,weight*mcreweight);
01788 }
01789 mcbeg++;
01790 }
01791
01792
01793 float chisq = Chi2Calc->GetChi2(DataHist[0],
01794 temp_sig_ND,temp_bkg_ND,
01795 MCPOT[0]/DataPOT[0]);
01796 chisq += Chi2Calc->GetChi2(DataHist[1],
01797 temp_sig_FD,temp_bkg_FD,
01798 MCPOT[1]/DataPOT[1]);
01799
01800
01801 if(EShiftErr>0) chisq += TMath::Power((1.-EShift)/EShiftErr,2);
01802 if(NgenMax[0]!=NgenMin[0])
01803 chisq += TMath::Power((NgenMean[0]-NgenPar0)/NgenErr[0],2);
01804 if(NgenMax[1]!=NgenMin[1])
01805 chisq += TMath::Power((NgenMean[1]-NgenPar1)/NgenErr[1],2);
01806 if(NgenMax[2]!=NgenMin[2])
01807 chisq += TMath::Power((NgenMean[2]-NgenPar2)/NgenErr[2],2);
01808
01809 if(chisq<ChiSq) {
01810 delete TempBestFit[0];
01811 TempBestFit[0] = temp_sig_ND;
01812 TempBestFit[0]->SetName("TempNearBestFit");
01813 delete TempBestFit[1];
01814 TempBestFit[1] = temp_sig_FD;
01815 TempBestFit[1]->SetName("TempFarBestFit");
01816 ChiSq = chisq;
01817 }
01818 else {
01819 delete temp_sig_ND;
01820 delete temp_sig_FD;
01821 }
01822 delete temp_bkg_ND;
01823 delete temp_bkg_FD;
01824 }
01825 }
01826 }
01827 }
01828 if(ChiSq<ChiMin){
01829 delete BestFit[0];
01830 delete BestFit[1];
01831 ChiMin=ChiSq;
01832 Par1MinVal=par1;
01833 Par2MinVal=par2;
01834 TempBestFit[0]->Scale(DataPOT[0]/MCPOT[0]);
01835 TempBestFit[0]->SetName("NearBestFit");
01836 BestFit[0] = TempBestFit[0];
01837 TempBestFit[1]->Scale(DataPOT[1]/MCPOT[1]);
01838 TempBestFit[1]->SetName("FarBestFit");
01839 BestFit[1] = TempBestFit[1];
01840 }
01841 else {
01842 delete TempBestFit[0];
01843 delete TempBestFit[1];
01844 }
01845 Chi2Surf->Fill(par1,par2,ChiSq);
01846 }
01847 }
01848 std::cout << "Minimum value of Chi2 is " << ChiMin << std::endl;
01849 std::cout << "at Par1 = " << Par1MinVal << " and Par2 = "
01850 << Par2MinVal << std::endl;
01851
01852 return true;
01853 }
01854
01855
01856 void MadAnalysis::EndJob()
01857 {
01858
01859 std::string outputFileName = "MadFile_" + tag + ".root";
01860 TFile *outputFile = new TFile(outputFileName.c_str(),"RECREATE");
01861 outputFile->cd();
01862 for(int i=0;i<2;i++){
01863 if(MCSignalHist[i]) MCSignalHist[i]->Write();
01864 if(MCBkgdHist[i]) MCBkgdHist[i]->Write();
01865 if(DataHist[i]) DataHist[i]->Write();
01866 if(DataSignalHist[i]) DataSignalHist[i]->Write();
01867 if(DataBkgdHist[i]) DataBkgdHist[i]->Write();
01868 if(BestFit[i]) BestFit[i]->Write();
01869 }
01870 if(Chi2Surf) Chi2Surf->Write();
01871 TTree *varTree = new TTree("Params","Params");
01872 varTree->Branch("DataPOTNear",&DataPOT[0],"DataPOTNear/F",32000);
01873 varTree->Branch("DataPOTFar",&DataPOT[1],"DataPOTFar/F",32000);
01874 varTree->Branch("numDataEvtsNear",&numDataEvts[0],"numDataEvtsNear/I",32000);
01875 varTree->Branch("numDataEvtsFar",&numDataEvts[1],"numDataEvtsFar/I",32000);
01876 varTree->Branch("signalFracInDataNear",&signalFracInData[0],
01877 "signalFracInDataNear/F",32000);
01878 varTree->Branch("signalFracInDataFar",&signalFracInData[1],
01879 "signalFracInDataFar/F",32000);
01880
01881 varTree->Branch("MCPOTNear",&MCPOT[0],"MCPOTNear/F",32000);
01882 varTree->Branch("MCPOTFar",&MCPOT[1],"MCPOTFar/F",32000);
01883 varTree->Branch("numMCEvtsNear",&numMCEvts[0],"numMCEvtsNear/I",32000);
01884 varTree->Branch("numMCEvtsFar",&numMCEvts[1],"numMCEvtsFar/I",32000);
01885
01886 varTree->Branch("NumOscPars",&NumOscPars,"NumOscPars/I",32000);
01887 varTree->Branch("OscillationParameters",OscillationParameters,
01888 "OscillationParameters[10]/D",32000);
01889
01890 varTree->Branch("ChiMin",&ChiMin,"ChiMin/F",32000);
01891 varTree->Branch("NDOF",&NDOF,"NDOF/F",32000);
01892 varTree->Branch("Par1MinVal",&Par1MinVal,"Par1MinVal/F",32000);
01893 varTree->Branch("Par2MinVal",&Par2MinVal,"Par2MinVal/F",32000);
01894
01895 varTree->Branch("EShiftErr",&EShiftErr,"EShiftErr/F",32000);
01896
01897 double NormErr = this->GetChi2Calc()->GetNormalisationError();
01898 varTree->Branch("NormErr",&NormErr,"NormErr/D",32000);
01899 double BinErr = this->GetChi2Calc()->GetBinToBinError();
01900 varTree->Branch("BinErr",&BinErr,"BinErr/D",32000);
01901 double CalcType = this->GetChi2Calc()->GetCalcType();
01902 varTree->Branch("CalcType",&CalcType,"CalcType/D",32000);
01903 double BkdXSecErr = this->GetChi2Calc()->GetBkdXSecError();
01904 varTree->Branch("BkdXSecErr",&BkdXSecErr,"BkdXSecErr/D",32000);
01905
01906 varTree->Fill();
01907 varTree->Write();
01908 outputFile->Close();
01909
01910 }
01911 #endif // #ifdef madanalysis_cxx
