loon: NCExtrapolationNS.cxx Source File

00001 
00002 //$Id: NCExtrapolationNS.cxx,v 1.43 2008/08/25 16:23:44 bckhouse Exp $
00003 //
00004 //NCExtrapolationNS.cxx
00005 //
00006 //class to hold pdfs for doing nccc separation
00007 //
00008 //B. Rebel 2/2007
00010 
00011 #include "NCUtils/Extrapolation/NCExtrapolationNS.h"
00012 #include "NCUtils/Extrapolation/NCEnergyBin.h"
00013 #include "NCUtils/Cuts/NCAnalysisCutsNC.h"
00014 #include "MessageService/MsgService.h"
00015 #include "AnalysisNtuples/ANtpAnalysisInfo.h"
00016 #include "AnalysisNtuples/ANtpBeamInfo.h"
00017 #include "AnalysisNtuples/ANtpRecoInfo.h"
00018 #include "AnalysisNtuples/ANtpHeaderInfo.h"
00019 #include "AnalysisNtuples/ANtpTruthInfoBeam.h"
00020 
00021 #include "TH2.h"
00022 #include "TSystem.h"
00023 #include "TDirectory.h"
00024 #include "TMath.h"
00025 #include <cassert>
00026 #include <fstream>
00027 
00028 using namespace NCType;
00029 
00030 // These three values correspond to shower_shift loops//
00031 const int     n1     = 21;
00032 const int     n0     = 0;
00033 const double  les    = 0.80;
00034 const double  hes    = 1.20;
00035 const int     n1d    = 20;
00036 const double size_e=(hes-les)/n1d;
00037 // These three values correspond to nc shift loops//
00038 const int     nc1    = 21;
00039 const int     nc0    = 0;
00040 const double  lnc    = 0.20;
00041 const double  hnc    = 1.80;
00042 const int     ncd    = 20;
00043 const double size_nc=(hnc-lnc)/ncd;
00044 // These three values correspond to normalization shift loops//
00045 const int     nr1    = 11;//21;
00046 const int     nr0    = 10;//0;
00047 const double  lnr    = 0.95;
00048 const double  hnr    = 1.05;
00049 const int     nrd    = 20;
00050 const double size_nr=(hnr-lnr)/nrd;
00051 
00052 // These are the sigmas//
00053 const double signorm = NCType::kParameterSigmas[NCType::kNormalization];
00054 const double siges   = NCType::kParameterSigmas[NCType::kShowerEnergy];
00055 const double signc   = NCType::kParameterSigmas[NCType::kNCBackground];
00056 
00057 //far and near detector masses
00058 const double kmass_f    = 52.97*486*0.0254*7.874+52.97*484*0.011*1.19879;
00059 const double kmass_n_cc = 3.141592654*1.*1.*(1.710210*7.874+0.746800*1.20317); //reference to cc analysis
00060 const double kmass_n    = 2.4930144*(1.710210*7.874+0.746800*1.20317)*(4.7368-1.728)/4.; //for nc analysis
00061 const double kscalef    = 1.00;
00062 
00063 CVSID("$Id: NCExtrapolationNS.cxx,v 1.43 2008/08/25 16:23:44 bckhouse Exp $");
00064 
00065 //......................................................................
00066 NCExtrapolationNS::NCExtrapolationNS() :
00067   NCExtrapolation::NCExtrapolation()
00068 {
00069 }
00070 
00071 //.......................................................................
00072 NCExtrapolationNS::NCExtrapolationNS(NCAnalysisUtils *utils,
00073                                      std::map<Int_t,double>& dataNear,
00074                                      std::map<Int_t,double>& mcNear,
00075                                      std::map<Int_t,double>& dataFar,
00076                                      std::map<Int_t,double>& mcFar,
00077                                      bool useCC,
00078                                      bool useNC) :
00079   NCExtrapolation(utils, dataNear, mcNear, dataFar, mcFar, useCC, useNC)
00080 {
00081   //set up the chi^2 histograms
00082   fchisqall   = new TH2F("chisqall","",
00083                          NCType::kNumUMu3SqrBins,
00084                          NCType::kUMu3SqrStart,
00085                          NCType::kUMu3SqrEnd,
00086                          NCType::kNumDeltaMSqrBins,
00087                          NCType::kDeltaMSqrStart,
00088                          NCType::kDeltaMSqrEnd);
00089   fchisqallnc = new TH2F("chisqallnc","",
00090                          NCType::kNumUS3SqrBins,
00091                          NCType::kUS3SqrStart,
00092                          NCType::kUS3SqrEnd,
00093                          NCType::kNumDeltaMSqrBins,
00094                          NCType::kDeltaMSqrStart,
00095                          NCType::kDeltaMSqrEnd);
00096 
00097 }
00098 
00099 //----------------------------------------------------------------------
00100 NCExtrapolationNS::~NCExtrapolationNS()
00101 {
00102 }
00103 
00104 //----------------------------------------------------------------------
00105 //this method allows the user to pass in appropriate settings for the
00106 //extrapolation such as locations of files, whether to generate the files
00107 //again, etc.  the registry comes from the job module GetConfig method
00108 void NCExtrapolationNS::Prepare(const Registry& r)
00109 {
00110   MSG("NCExtrapolationNS", Msg::kWarning) << "NCExtrapolationNS::Prepare(const Registry& r)" << endl;
00111   int         tmpb;  // a temp bool. See comment under DefaultConfig...
00112   int         tmpi;  // a temp int.
00113   double      tmpd;  // a temp double.
00114   const char* tmps;  // a temp string
00115 
00116   NCExtrapolation::Prepare(r);
00117 
00118   if(r.Get("Decay",                    tmpi)) fDecay                    = tmpi;
00119   if(r.Get("DoMC",                     tmpi)) fDoMC                     = tmpi;
00120   if(r.Get("BeamFileName",             tmps)) fBeamFileName             = tmps;
00121   if(r.Get("BeamFileNamePostShutDown", tmps)) fBeamFileNamePostShutDown = tmps;
00122   if(r.Get("TrueCorrectionHistNear",   tmps)) fTrueCorrectionHistNear   = tmps;
00123   if(r.Get("TrueCorrectionHistFar",    tmps)) fTrueCorrectionHistFar    = tmps;
00124   if(r.Get("FarTrueCorrectionFile",    tmps)) fFarTrueCorrectionFile    = tmps;
00125   if(r.Get("NearTrueCorrectionFile",   tmps)) fNearTrueCorrectionFile   = tmps;
00126 
00127   TString topDir="NCUtils/data";
00128   TString base="";
00129   base=getenv("SRT_PRIVATE_CONTEXT");
00130   if(base!="" && base!="."){
00131     // check if directory exists in SRT_PRIVATE_CONTEXT
00132     TString path = base + "/" + topDir;
00133     void *dir_ptr = gSystem->OpenDirectory(path);
00134     if(!dir_ptr) base=getenv("SRT_PUBLIC_CONTEXT");
00135   }
00136   // if it doesn't exist then use SRT_PUBLIC_CONTEXT
00137   else base=getenv("SRT_PUBLIC_CONTEXT");
00138 
00139   if(base==""){
00140     MSG("NCExtrapolationNS",Msg::kFatal) << "NCExtrapolationNS: No SRT_PUBLIC_CONTEXT set"
00141                                          << endl;
00142     assert(false);
00143   }
00144   topDir += "/";
00145   topDir.Prepend(base+"/");
00146   MSG("NCExtrapoationNS", Msg::kInfo)
00147     << "NCExtrapoationNS reading files from: " << topDir << endl;
00148 
00149   fBeamFileName.Prepend(topDir);
00150   fBeamFileNamePostShutDown.Prepend(topDir);
00151   fFarTrueCorrectionFile.Prepend(topDir);
00152   fNearTrueCorrectionFile.Prepend(topDir);
00153 
00154   MakeHistograms();
00155 
00156   tmpb = 0;
00157   tmpi = 0;
00158   tmpd = 0.;
00159   tmps = "";
00160 
00161   return;
00162 }
00163 
00164 //----------------------------------------------------------------------
00165 void NCExtrapolationNS::MakeHistograms()
00166 {
00167   fbinvect[0]=0.;
00168   fbinvect[1]=0.5;
00169   for(int i=2;i<200;i++){
00170     fbinvect[i]=0.25+(i)*0.25;
00171     if(i==199) fbinvect[i]=200.;
00172   }
00173 
00174   //save the current directory before opening some files
00175   TDirectory *savedir = gDirectory;
00176 
00177   TFile beamfile(fBeamFileName,"READ");
00178   fbmat_norm_r2 = dynamic_cast<TH2D *>(beamfile.Get("bmat_norm_r2"));
00179   fbmat_norm_r2->SetDirectory(savedir);
00180   TFile beamfile_post(fBeamFileNamePostShutDown,"READ");
00181   fbmat_norm_r2_post = dynamic_cast<TH2D *>(beamfile.Get("bmat_norm_r2"));
00182   fbmat_norm_r2_post->SetDirectory(savedir);
00183 
00184   savedir->cd();
00185 
00186   //these histograms will be made by the method now
00187   fSingleMCNearNuCC = new TH1D("singleMCNearNuCC", "", fnbins, fbinvect);
00188   fSingleMCFarNuCC = new TH1D("singleMCFarNuCC", "", fnbins, fbinvect);
00189   fSingleMCNearAllCC = new TH1D("singleMCNearAllCC", "", fnbins, fbinvect);
00190   fSingleMCFarAllCC = new TH1D("singleMCFarAllCC", "", fnbins, fbinvect);
00191   fMCNearNuCC = new TH1D("mcNearNuCC", "", fnbins, fbinvect);
00192   fMCFarNuCC = new TH1D("mcFarNuCC", "", fnbins, fbinvect);
00193   fMCNearAllCC = new TH1D("mcNearAllCC", "", fnbins, fbinvect);
00194   fMCFarAllCC = new TH1D("mcFarAllCC", "", fnbins, fbinvect);
00195   //LLH deprecating fcorf = new TH1D(fTrueCorrectionHistFar, "", fnbins, fbinvect); //for cc only - LLH
00196   //LLH deprecating fcor = new TH1D(fTrueCorrectionHistNear, "", fnbins, fbinvect); //for cc only - LLH
00197   fcorall = new TH1D("fcorall","",fnbins,fbinvect);  //for cc and nc? - LLH
00198   fcorfall = new TH1D("fcorfall","",fnbins,fbinvect);  //for cc and nc? - LLH
00199 
00200 //   TFile corrf(fFarTrueCorrectionFile,"READ");
00201 //   fcorf = dynamic_cast<TH1D *>(corrf.Get(fTrueCorrectionHistFar));
00202 //   fcorf->SetDirectory(savedir);
00203 //   TFile corr(fNearTrueCorrectionFile,"READ");
00204 //   fcor = dynamic_cast<TH1D *>(corr.Get(fTrueCorrectionHistNear));
00205 //   fcor->SetDirectory(savedir);
00206 
00207   fTrueND = new TH1F("trueND", "", fnbins, fbinvect);
00208   fTrueNDAfterRT = new TH1F("trueNDAfterRT", "", fnbins, fbinvect);
00209   fTrueNDAfterRT_noeffcorr = new TH1F("trueNDAfterRT_noeffcorr", "", fnbins, fbinvect); //LLH adding
00210   fRecoND = new TH1F("recoND", "", fnbins, fbinvect);
00211   fTrueFD = new TH1F("trueFD", "", fnbins, fbinvect);
00212   fTrueFD_nc = new TH1F("trueFD_nc", "", fnbins, fbinvect);
00213   fRecoFD = new TH1F("recoFD", "", fnbins, fbinvect);
00214   fRecoFD_nc = new TH1F("recoFD_nc", "", fnbins, fbinvect);
00215   fTrueFDAfterBM = new TH1F("trueFDAfterBM", "", fnbins, fbinvect);
00216   fTrueFDAfterBM_noacccorr = new TH1F("trueFDAfterBM_noacccorr", "", fnbins, fbinvect); //LLH adding
00217   fTrueFDAfterBM_posteffcorr = new TH1F("trueFDAfterBM_posteffcorr", "", fnbins, fbinvect); //LLH adding
00218   fTrueFDAfterOsc = new TH1F("trueFDAfterOsc", "", fnbins, fbinvect); //LLH adding
00219   fRecoFDAfterBM = new TH1F("recoFDAfterBM", "", fnbins, fbinvect);
00220 //  fdatsel = new TH1F("datselstore", "", fnbins, fbinvect); //LLH adding, store before POT scaling
00221 
00222   //used in efficiency calculation in ND
00223   fallccmc = new TH1F("allccmc","",fnbins,fbinvect);
00224   fallccmcf = new TH1F("allccmcf","",fnbins,fbinvect); //LLH adding for proper eff calc
00225   fallncmcf = new TH1F("allncmcf","",fnbins,fbinvect); //LLH adding for proper eff calc
00226 
00227   // FAR CC AND NC TRUE-RECO
00228 //not used-LLH ftruerecomat = new TH2F("truerecomat", "", fnbins, fbinvect, fnbins, fbinvect);
00229 //not used-LLH ftruerecomatnc = new TH2F("truerecomatnc", "", fnbins, fbinvect, fnbins, fbinvect);
00230 //not used-LLH ftruerecomatb = new TH2F("truerecomatb",  "", fnbins, fbinvect, fnbins, fbinvect);
00231 //not used-LLH ftruerecomatncb = new TH2F("truerecomatncb","", fnbins, fbinvect, fnbins, fbinvect);
00232 
00233   // FAR CC AND NC SELECTION EFFICIENCIES (SIGNAL AND BACKGROUND)
00234   feffic = new TH1F("effic",   "",fnbins, fbinvect);
00235   fefficnc = new TH1F("efficnc", "",fnbins, fbinvect);
00236   fefficbg = new TH1F("efficb",  "",fnbins, fbinvect);
00237   fefficncbg = new TH1F("efficncb","",fnbins, fbinvect);
00238   fratccnc = new TH1F("ratccnc", "",fnbins, fbinvect);
00239   feffic_numcheck = new TH1F("effic_numcheck",   "",fnbins, fbinvect); //LLH adding to check eff calc
00240 
00241   // TAUS
00242   ftruerecomatt = new TH2F("truerecomatt", "", fnbins, fbinvect, fnbins, fbinvect);
00243   ftruerecomatbt = new TH2F("truerecomatbt", "", fnbins, fbinvect, fnbins, fbinvect);
00244   feffict = new TH1F("effict", "", fnbins, fbinvect);
00245   fefficbgt = new TH1F("efficbt", "", fnbins, fbinvect);
00246   frattau = new TH1F("rattau", "", fnbins, fbinvect);
00247 
00248   // NEAR
00249   frecotruemat            = new TH2F("recotruemat","",fnbins,fbinvect,fnbins,fbinvect);
00250   fpurityn                = new TH1F("purityn","",fnbins,fbinvect);
00251   fefficiencyn            = new TH1F("efficiencyn","",fnbins,fbinvect);
00252 
00253   // FAR CC AND NC TRUE-RECO  POST
00254 //not used-LLH  ftruerecomat_post            = new TH2F("truerecomat_post",   "",fnbins,fbinvect,fnbins,fbinvect);
00255 //not used-LLH  ftruerecomatnc_post          = new TH2F("truerecomatnc_post", "",fnbins,fbinvect,fnbins,fbinvect);
00256 //not used-LLH  ftruerecomatb_post           = new TH2F("truerecomatb_post",  "",fnbins,fbinvect,fnbins,fbinvect);
00257 //not used-LLH  ftruerecomatncb_post         = new TH2F("truerecomatncb_post","",fnbins,fbinvect,fnbins,fbinvect);
00258 
00259   // FAR CC AND NC SELECTION EFFICIENCIES (SIGNAL AND BACKGROUND)
00260   feffic_post                  = new TH1F("effic_post",   "",fnbins,fbinvect);
00261   fefficnc_post                = new TH1F("efficnc_post", "",fnbins,fbinvect);
00262   fefficbg_post                 = new TH1F("efficb_post",  "",fnbins,fbinvect);
00263   fefficncbg_post               = new TH1F("efficncb_post","",fnbins,fbinvect);
00264   fratccnc_post                = new TH1F("ratccnc_post", "",fnbins,fbinvect);
00265   // TAUS
00266   ftruerecomatt_post           = new TH2F("truerecomatt_post",   "",fnbins,fbinvect,fnbins,fbinvect);
00267   ftruerecomatbt_post          = new TH2F("truerecomatbt_post",  "",fnbins,fbinvect,fnbins,fbinvect);
00268   feffict_post                 = new TH1F("effict_post",   "",fnbins,fbinvect);
00269   fefficbgt_post                = new TH1F("efficbt_post",  "",fnbins,fbinvect);
00270   frattau_post                 = new TH1F("rattau_post","",fnbins,fbinvect);
00271 
00272   // Near
00273   frecotruemat_post            = new TH2F("recotruemat_post","",fnbins,fbinvect,fnbins,fbinvect);
00274   fpurityn_post                = new TH1F("purityn_post","",fnbins,fbinvect);
00275   fefficiencyn_post            = new TH1F("efficiencyn_post","",fnbins,fbinvect);
00276 
00277   //PREDICTED SPECTRUM CC - bestfit sys unosc
00278   //fnoosctruereal           = new TH1F("noosctruereal","",fnbins,fbinvect); //not being used - LLH
00279   fnoosctrue               = new TH1F("noosctrue","",fnbins,fbinvect); //unosc prediction cc, trueE
00280   fnooscreco               = new TH1F("nooscreco","",fnbins,fbinvect); //unosc prediction cc, recoE
00281   fnooscreconc             = new TH1F("nooscreconc","",fnbins,fbinvect); //unosc prediction nc bg, recoE
00282   fnooscrecotau            = new TH1F("nooscrecotau","",fnbins,fbinvect);
00283 
00284   //PREDICTED SPECTRUM NC - bestfit sys unosc
00285   fnoosctrue_ncfit             = new TH1F("noosctrue_ncfit","",fnbins,fbinvect); //unosc prediction nc, trueE, no sel eff!!!
00286   fnooscreco_ncfit            = new TH1F("nooscreco_ncfit","",fnbins,fbinvect); //unosc preditcion nc, recoE
00287   fnooscrecocc_ncfit             = new TH1F("nooscrecocc_ncfit","",fnbins,fbinvect); //onosc precition cc bg, recoE
00288   fnooscrecotau_ncfit          = new TH1F("nooscrecotau_ncfit","",fnbins,fbinvect);
00289 
00290   //PREDICTED SPECTRUM CC POST
00291   //fnoosctruereal_post          = new TH1F("noosctruereal_post","",fnbins,fbinvect);  //not being used - LLH
00292   fnoosctrue_post              = new TH1F("noosctrue_post","",fnbins,fbinvect);
00293   fnooscreco_post              = new TH1F("nooscreco_post","",fnbins,fbinvect);
00294   fnooscreconc_post            = new TH1F("nooscreconc_post","",fnbins,fbinvect);
00295   fnooscrecotau_post           = new TH1F("nooscrecotau_post","",fnbins,fbinvect);
00296 
00297   //PREDICTED SPECTRUM NC post
00298   fnoosctrue_ncfit_post            = new TH1F("noosctrue_ncfit_post","",fnbins,fbinvect);
00299   fnooscreco_ncfit_post           = new TH1F("nooscreco_ncfit_post","",fnbins,fbinvect);
00300   fnooscrecocc_ncfit_post            = new TH1F("nooscrecocc_ncfit_post","",fnbins,fbinvect);
00301   fnooscrecotau_ncfit_post         = new TH1F("nooscrecotau_ncfit_post","",fnbins,fbinvect);
00302 
00303   // auxilary spectra
00304   fnooscar              = new TH1F("nooscar","",fnbins,fbinvect);
00305   fnooscar_ccselcc      = new TH1F("nooscar_ccselcc","",fnbins,fbinvect);
00306   fnooscar_ccselnc         = new TH1F("nooscar_ccselnc","",fnbins,fbinvect);
00307   fnooscar_ncselnc      = new TH1F("nooscar_ncselnc","",fnbins,fbinvect);
00308   fnooscar_ncselcc         = new TH1F("nooscar_ncselcc","",fnbins,fbinvect);
00309   //Taus
00310   fnooscar_ccselcct        = new TH1F("nooscar_ccselcct","",fnbins,fbinvect);
00311   fnooscar_ccselnct     = new TH1F("nooscar_ccselnct","",fnbins,fbinvect);
00312 
00313   //POST
00314   fnooscar_post                = new TH1F("nooscar_post","",fnbins,fbinvect);
00315   fnooscar_ccselcc_post         = new TH1F("nooscar_ccselcc_post","",fnbins,fbinvect);
00316   fnooscar_ccselnc_post         = new TH1F("nooscar_ccselnc_post","",fnbins,fbinvect);
00317   fnooscar_ncselnc_post         = new TH1F("nooscar_ncselnc_post","",fnbins,fbinvect);
00318   fnooscar_ncselcc_post         = new TH1F("nooscar_ncselcc_post","",fnbins,fbinvect);
00319   //Taus
00320   fnooscar_ccselcct_post       = new TH1F("nooscar_ccselcct_post","",fnbins,fbinvect);
00321   fnooscar_ccselnct_post       = new TH1F("nooscar_ccselnct_post","",fnbins,fbinvect);
00322 
00323   // BEST FIT TRUE AND RECO FROM FIT  CC
00324 
00325   fbestfitreco            = new TH1F("bestfitreco","",fnbins,fbinvect);
00326   fbestfitreconc          = new TH1F("bestfitreconc","",fnbins,fbinvect);
00327   fbestfitrecotau         = new TH1F("bestfitrecotau","",fnbins,fbinvect);
00328 
00329   // BEST FIT TRUE AND RECO FROM FIT NC
00330 
00331   fbestfitreco_ncfit         = new TH1F("bestfitreco_ncfit","",fnbins,fbinvect);
00332   fbestfitrecocc_ncfit          = new TH1F("bestfitrecocc_ncfit","",fnbins,fbinvect);
00333   fbestfitrecotau_ncfit      = new TH1F("bestfitrecotau_ncfit","",fnbins,fbinvect);
00334 
00335   // BEST FIT TRUE AND RECO FROM FIT CC
00336 
00337   fbestfitreco_post         = new TH1F("bestfitreco_post","",fnbins,fbinvect);
00338   fbestfitreconc_post       = new TH1F("bestfitreconc_post","",fnbins,fbinvect);
00339   fbestfitrecotau_post      = new TH1F("bestfitrecotau_post","",fnbins,fbinvect);
00340 
00341   // BEST FIT TRUE AND RECO FROM FIT NC
00342 
00343   fbestfitreco_ncfit_post      = new TH1F("bestfitreco_ncfit_post","",fnbins,fbinvect);
00344   fbestfitrecocc_ncfit_post       = new TH1F("bestfitrecocc_ncfit_post","",fnbins,fbinvect);
00345   fbestfitrecotau_ncfit_post   = new TH1F("bestfitrecotau_ncfit_post","",fnbins,fbinvect);
00346 
00347   // HELP ON FITTING CC
00348   ftemptrue               = new TH1F("temptrue","",fnbins,fbinvect);
00349   ftempreco               = new TH1F("tempreco","",fnbins,fbinvect);
00350   ftempreconc             = new TH1F("tempreconc","",fnbins,fbinvect);
00351   ftemprecotau            = new TH1F("temprecotau","",fnbins,fbinvect);
00352 
00353   // HELP ON FITTING NC
00354   ftemptrue_ncfit            = new TH1F("temptrue_ncfit","",fnbins,fbinvect);
00355   ftempreco_ncfit            = new TH1F("tempreco_ncfit","",fnbins,fbinvect);
00356   ftemprecocc_ncfit             = new TH1F("temprecocc_ncfit","",fnbins,fbinvect);
00357   ftemprecotau_ncfit         = new TH1F("temprecotau_ncfit","",fnbins,fbinvect);
00358 
00359   // HELP ON FITTING CC
00360   ftemptrue_post            = new TH1F("temptrue_post","",fnbins,fbinvect);
00361   ftempreco_post            = new TH1F("tempreco_post","",fnbins,fbinvect);
00362   ftempreconc_post          = new TH1F("tempreconc_post","",fnbins,fbinvect);
00363   ftemprecotau_post         = new TH1F("temprecotau_post","",fnbins,fbinvect);
00364 
00365   // HELP ON FITTING NC
00366   ftemptrue_ncfit_post         = new TH1F("temptrue_ncfit_post","",fnbins,fbinvect);
00367   ftempreco_ncfit_post         = new TH1F("tempreco_ncfit_post","",fnbins,fbinvect);
00368   ftemprecocc_ncfit_post          = new TH1F("temprecocc_ncfit_post","",fnbins,fbinvect);
00369   ftemprecotau_ncfit_post      = new TH1F("temprecotau_ncfit_post","",fnbins,fbinvect);
00370 
00371   int nbins1 = 15;
00372   if(fDoMC) nbins1 = 39;
00373 
00374   float binvect1[40];
00375   binvect1[0]=0;
00376   binvect1[1]=1.;
00377   for (int i=2; i < nbins1+1; i++){
00378     if(fDoMC > 0 && i <= 19)  binvect1[i]=0.5+(i)*0.5;
00379     if(fDoMC > 0 && i > 19)   binvect1[i]=10.+(i-19)*1.;
00380     if(fDoMC < 1 && i <= 10)  binvect1[i]=(i)*1.;
00381     if(fDoMC < 1 && i > 10)   binvect1[i]=10.+(i-10)*2.;
00382   }
00383 
00384   fbestfitrecotau_rebinned        = new TH1F("bestfitrecotau_rebinned","",
00385                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00386   fbestfitreco_rebinned           = new TH1F("bestfitreco_rebinned","",
00387                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00388   fbestfitreconc_rebinned         = new TH1F("bestfitreconc_rebinned","",
00389                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00390   foscreco_rebinned                  = new TH1F("oscreco_rebinned","",
00391                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00392   ftempreco_rebinned              = new TH1F("tempreco_rebinned","",
00393                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00394   ftempreconc_rebinned            = new TH1F("tempreconc_rebinned","",
00395                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00396   fnooscreco_rebinned             = new TH1F("nooscreco_rebinned","",
00397                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00398   fnooscreconc_rebinned           = new TH1F("nooscreconc_rebinned","",
00399                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00400   fnooscrecotau_rebinned          = new TH1F("nooscrecotau_rebinned","",
00401                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00402 
00403   fbestfitrecotau_rebinned_post       = new TH1F("bestfitrecotau_rebinned_post","",
00404                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00405   fbestfitreco_rebinned_post          = new TH1F("bestfitreco_rebinned_post","",
00406                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00407   fbestfitreconc_rebinned_post        = new TH1F("bestfitreconc_rebinned_post","",
00408                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00409   foscreco_rebinned_post                     = new TH1F("oscreco_rebinned_post","",
00410                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00411   ftempreco_rebinned_post             = new TH1F("tempreco_rebinned_post","",
00412                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00413   ftempreconc_rebinned_post           = new TH1F("tempreconc_rebinned_post","",
00414                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00415   fnooscreco_rebinned_post            = new TH1F("nooscreco_rebinned_post","",
00416                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00417   fnooscreconc_rebinned_post          = new TH1F("nooscreconc_rebinned_post","",
00418                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00419   fnooscrecotau_rebinned_post         = new TH1F("nooscrecotau_rebinned_post","",
00420                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00421 
00422   //NC
00423   fbestfitrecotau_ncfit_rebinned     = new TH1F("bestfitrecotau_ncfit_rebinned","",
00424                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00425   fbestfitreco_ncfit_rebinned        = new TH1F("bestfitreco_ncfit_rebinned","",
00426                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00427   fbestfitrecocc_ncfit_rebinned         = new TH1F("bestfitrecocc_ncfit_rebinned","",
00428                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00429   foscreco_ncfit_rebinned            = new TH1F("oscreco_ncfit_rebinned","",
00430                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00431   ftempreco_ncfit_rebinned           = new TH1F("tempreco_ncfit_rebinned","",
00432                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00433   ftemprecocc_ncfit_rebinned            = new TH1F("temprecocc_ncfit_rebinned","",
00434                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00435   fnooscreco_ncfit_rebinned          = new TH1F("nooscreco_ncfit_rebinned","",
00436                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00437   fnooscrecocc_ncfit_rebinned           = new TH1F("nooscrecocc_ncfit_rebinned","",
00438                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00439   fnooscrecotau_ncfit_rebinned        = new TH1F("nooscrecotau_ncfit_rebinned","",
00440                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00441 
00442 
00443   fbestfitrecotau_ncfit_rebinned_post    = new TH1F("bestfitrecotau_ncfit_rebinned_post","",
00444                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00445   fbestfitreco_ncfit_rebinned_post       = new TH1F("bestfitreco_ncfit_rebinned_post","",
00446                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00447   fbestfitrecocc_ncfit_rebinned_post        = new TH1F("bestfitrecovc_rebinned_post","",
00448                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00449   foscreco_ncfit_rebinned_post       = new TH1F("oscreco_ncfit_rebinned_post","",
00450                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00451   ftempreco_ncfit_rebinned_post          = new TH1F("tempreco_ncfit_rebinned_post","",
00452                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00453   ftemprecocc_ncfit_rebinned_post           = new TH1F("temprecocc_ncfit_rebinned_post","",
00454                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00455   fnooscreco_ncfit_rebinned_post         = new TH1F("nooscreco_ncfit_rebinned_post","",
00456                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00457   fnooscrecocc_ncfit_rebinned_post          = new TH1F("nooscrecocc_ncfit_rebinned_post","",
00458                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00459   fnooscrecotau_ncfit_rebinned_post       = new TH1F("nooscrecotau_ncfit_rebinned_post","",
00460                                         NCType::kNumEnergyBinsFar, NCType::kEnergyBinsFar);
00461 
00462   return;
00463 }
00464 
00465 //----------------------------------------------------------------------
00466 //add events to the histograms used for the weighting
00467 void NCExtrapolationNS::AddEvent(ANtpHeaderInfo    *headerInfo,
00468                                  ANtpBeamInfo      *beamInfo,
00469                                  ANtpRecoInfo      *recoInfo,
00470                                  ANtpAnalysisInfo  *analysisInfo,
00471                                  ANtpTruthInfoBeam *truthInfo,
00472                                  NCType::ERunType runType,
00473                                  bool useMCAsData,
00474                                  NCType::EFileType fileType)
00475 {
00476 
00477   //++++Energies are set to the correct NC/CC values in the NCExtrapolationModule
00478   //before being passed into the extrapolations
00479 
00480   //only use the L010z185i beam
00481   if(beamInfo->beamType != BeamType::ToZarko(BeamType::kL010z185i))
00482     return;
00483 
00484   if(headerInfo->detector == (int)Detector::kNear)
00485     AddNearEvent(headerInfo, beamInfo, recoInfo, analysisInfo,
00486                  truthInfo, runType, useMCAsData, fileType);
00487 
00488   if(headerInfo->detector == (int)Detector::kFar)
00489     AddFarEvent(headerInfo, beamInfo, recoInfo, analysisInfo,
00490                 truthInfo, runType, useMCAsData, fileType);
00491 
00492   return;
00493 }
00494 
00495 //----------------------------------------------------------------------
00496 void NCExtrapolationNS::AddNearEvent(ANtpHeaderInfo *headerInfo,
00497                                      ANtpBeamInfo *beamInfo,
00498                                      ANtpRecoInfo *recoInfo,
00499                                      ANtpAnalysisInfo *analysisInfo,
00500                                      ANtpTruthInfoBeam *truthInfo,
00501                                      NCType::ERunType runType,
00502                                      bool useMCAsData,
00503                                      NCType::EFileType fileType)
00504 {
00505 
00506   //fill histogram for all MC events with a true vertex in the fiducial volume
00507 
00508   //Note: if this code is run on stripped uDST files, it will miss events that were reco'd
00509   //reco'd outside the FV!  Make sure to run on unstripped files or find a way to apply
00510   //a correction... This should be a small effect - LLH
00511   if(headerInfo->dataType == (int)SimFlag::kMC && !useMCAsData){
00512     NCAnalysisCutsNC cuts;
00513     cuts.SetInfoObjects(headerInfo, beamInfo, 0, 0, 0, truthInfo);
00514     if( cuts.InFiducialVolumeTrue() && headerInfo->events != 0
00515         && truthInfo->interactionType == 1
00516         && TMath::Abs(truthInfo->nuFlavor) == 14){
00517 
00518       if(truthInfo->nuFlavor > 0) fallccmc->Fill(truthInfo->nuEnergy, recoInfo->weight);
00519 
00520       //check if this is a new snarl
00521 //       if(headerInfo->newSnarl > 0) AddSingleMCEventsToHistogram(Detector::kNear);
00522 
00524       //Following lines were for fcorf calculation, since we can't do
00525       //this from uDST files, comment out for now - LLH
00527 
00528 //       fTrueNuEnergy.push_back(truthInfo->nuEnergy);
00529 //       fTrueNuFlavor.push_back(truthInfo->nuFlavor);
00530 //       fTrueNuComplete.push_back(truthInfo->eventCompleteness);
00531     }//end if true vertex is in fiducial volume
00532 
00533   }//end if mc
00534 
00535   //Here apply FV and cleaning cuts!
00536   NCExtrapolation::AddEvent(headerInfo, beamInfo, recoInfo, analysisInfo,
00537                             truthInfo, runType, useMCAsData, fileType);
00538 
00539   //dont bother with events outside the fiducial volume
00540 //   if(recoInfo->inFiducialVolume < 1) return;
00541 
00542 //   if(headerInfo->dataType == (int)SimFlag::kMC && !useMCAsData){
00543 
00544 //     feshnd.push_back(recoInfo->showerEnergyCC);
00545 //     femund.push_back(recoInfo->muEnergy);
00546 //     fenutnd.push_back(truthInfo->nuEnergy);
00547 //     fccnd.push_back(truthInfo->interactionType);
00548 //     if(analysisInfo->isNC > 0)
00549 //       fpidpnd.push_back(NCType::kNC);
00550 //     else if(analysisInfo->isCC > 0)
00551 //       fpidpnd.push_back(NCType::kCC);
00552 //     fflavnd.push_back(truthInfo->nuFlavor);
00553 //     finisnd.push_back(truthInfo->initialState);
00554 //     fweind.push_back(recoInfo->weight);
00555 
00556 //   }//end if MC
00557 
00558   return;
00559 }
00560 
00561 //----------------------------------------------------------------------
00562 //all events passed to this method are garaunteed to be cc events, have
00563 //their true interaction vertex within the fiducial volume and are either
00564 //numu or numubar
00565 
00566 //Note (9/10/07) This is deprecated for the moment since uDST's don't save
00567 //info on generated MC events that weren't found by reconstruction - LLH
00568 void NCExtrapolationNS::AddSingleMCEventsToHistogram(Detector::Detector_t det)
00569 {
00570   //double loop over mc events in snarl and see if there are any repeats
00571   //for each entry with duplicates, store the number of duplicates found with j>i
00572   //store only instances where number of additional duplicates is 0
00573   //Because we are saving true energies, there's no need to worry about event completeness
00574 
00575   //note that if cleaning cut are applied to all events in the strip file,
00576   //then this method will not find any duplicates! - LLH
00577   vector<int> repeated;
00578 
00579 
00580   for(uint i = 0; i < fTrueNuEnergy.size(); ++i){
00581     repeated.push_back(-1);
00582     uint duplicatesCtr = 0; //for each duplicate found for entry i, increment this by 1
00583 
00584     for(uint j = i+1; j < fTrueNuEnergy.size(); ++j) //only search through j > i
00585     {
00586        if(fTrueNuEnergy[i] == fTrueNuEnergy[j])
00587           duplicatesCtr++;
00588 
00589     }//end inner loop
00590 
00591     repeated[i] = duplicatesCtr; //storing number of duplicates found after entry i
00592 
00593     MSG("NCExtrapolationNS", Msg::kInfo) << "i= " << i  //LLH changed from kDebug, 100
00594                                          << ", repeated[i] = " << repeated[i]
00595                                          << ", duplicatesCtr = " << duplicatesCtr
00596                                          <<", energy = " << fTrueNuEnergy[i]
00597                                          << endl;
00598 
00599   }//end outer loop
00600 
00601   //loop over the array one more time and fill the histogram
00602   //with single events or the event with the largest completeness
00603   //if it is a duplicate
00604   for(uint i = 0; i < fTrueNuEnergy.size(); ++i){
00605 
00606     if(det == Detector::kNear){
00607        if(repeated[i] ==0) fSingleMCNearAllCC->Fill(fTrueNuEnergy[i]);  //all cc-like
00608        fMCNearAllCC->Fill(fTrueNuEnergy[i]);
00609        if(fTrueNuFlavor[i] == 14){
00610           if(repeated[i] == 0) fSingleMCNearNuCC->Fill(fTrueNuEnergy[i]);  //true numu cc-like
00611           fMCNearNuCC->Fill(fTrueNuEnergy[i]);
00612       }//end if numu
00613     }//end if near
00614 
00615     else if(det == Detector::kFar){
00616       if(repeated[i] ==0) fSingleMCFarAllCC->Fill(fTrueNuEnergy[i]);  //all cc-like
00617       fSingleMCFarAllCC->Fill(fTrueNuEnergy[i]);
00618       if(fTrueNuFlavor[i] == 14){
00619         if(repeated[i] == 0) fSingleMCFarNuCC->Fill(fTrueNuEnergy[i]);  //true numu cc-like
00620         fMCFarNuCC->Fill(fTrueNuEnergy[i]);
00621       }//end if numu
00622     }//end if far
00623   }//end loop over events
00624 
00625   //reset the arrays (for next snarl)
00626   fTrueNuEnergy.clear();
00627   fTrueNuFlavor.clear();
00628   fTrueNuComplete.clear();
00629 
00630   return;
00631 }
00632 
00633 //----------------------------------------------------------------------
00634 void NCExtrapolationNS::AddFarEvent(ANtpHeaderInfo *headerInfo,
00635                                     ANtpBeamInfo *beamInfo,
00636                                     ANtpRecoInfo *recoInfo,
00637                                     ANtpAnalysisInfo *analysisInfo,
00638                                     ANtpTruthInfoBeam *truthInfo,
00639                                     NCType::ERunType runType,
00640                                     bool useMCAsData,
00641                                     NCType::EFileType fileType)
00642 {
00643    MAXMSG("NCExtrapolationNS", Msg::kInfo, 1) << "AddFarEvent " << endl;
00644 
00645   //fill histogram for all MC events with a true vertex in the fiducial volume
00646    if(headerInfo->dataType == (int)SimFlag::kMC && !useMCAsData){
00647 
00648       //don't double count and make sure event passed reco because
00649       //reconstruction efficiency is accounted for in the corf function
00650       //(eventually take this out)
00651       if(headerInfo->newSnarl && headerInfo->events!=0)
00652       {
00653          //Far detector efficiency corrections must be wrt all (numu) events in far detector
00654          if(truthInfo->nuFlavor == 14)
00655          {
00656             if(truthInfo->interactionType == 1) fallccmcf->Fill(truthInfo->nuEnergy, recoInfo->weight);
00657 
00659             //Following lines were for fcorf calculation, since we can't do
00660             //this from uDST files, comment out for now - LLH
00662 
00663             //check if this is a new snarl
00664 //      if(headerInfo->newSnarl > 0) AddSingleMCEventsToHistogram(Detector::kFar);
00665 //       fTrueNuEnergy.push_back(truthInfo->nuEnergy);
00666 //       fTrueNuFlavor.push_back(truthInfo->nuFlavor);
00667 //       fTrueNuComplete.push_back(truthInfo->eventCompleteness);
00668          }//end if true vertex is a numu
00669 
00670          //For NC events store all nu types
00671          if(truthInfo->interactionType == 0) fallncmcf->Fill(truthInfo->nuEnergy, recoInfo->weight);
00672       } //end if nonzero events and is new snarl
00673 
00674   }//end if mc
00675 
00676 
00677   NCExtrapolation::AddEvent(headerInfo, beamInfo, recoInfo, analysisInfo,
00678                             truthInfo, runType, useMCAsData, fileType);
00679 
00680   return;
00681 }
00682 
00683 //----------------------------------------------------------------------
00684 //just have to multiply total MC events with true vertex in the
00685 //fiducial volume by the correction factor for split events
00686 void NCExtrapolationNS::PrepareNearDetector()
00687 {
00688   //get the correction histograms for the single events
00689 //  fcor->Divide(fSingleMCNearNuCC, fMCNearNuCC);
00690 //  fcorall->Divide(fSingleMCNearAllCC, fMCNearAllCC);
00691 
00692    cout <<"In PrepareNearDetector" << endl;
00693 
00694   //For now read in fcor and fcorall from Niki's files, since we don't have the
00695   //capability to do this from uDST's
00696 //  TFile *corr      = new TFile("$MINOS_DATA/d84/niki/mad_cedar_daikon/near/correction_plots/near/norm_fix_near.root","READ");
00697    TFile *corr      = new TFile("/afs/fnal.gov/files/data/minos/d79/llhsu/eventFindingEff/norm_fix_near_uDST_strip.root","READ"); //using my own file - LLH with zero event cut
00698   fcor = (TH1F*)(corr->Get("cor_neu"));
00699 
00700   //Correct the true CC events for reconstruction efficiencies and event splitting
00701   fallccmc->Multiply(fcor);
00702   return;
00703 }
00704 
00705 //----------------------------------------------------------------------
00706 //this is where you call the guts of the code to do the prediction and
00707 //fit
00708 void NCExtrapolationNS::PredictSpectrum()
00709 {
00710 
00711    cout <<"In NCExtrapolationNS::PredictSpectrum!" << endl;
00712 
00713   //make the correction histograms used by the method
00714 //LLH deprecating  fcorf->Divide(fSingleMCFarNuCC, fMCFarNuCC);
00715 //LLH deprecating  fcorfall->Divide(fSingleMCFarAllCC, fMCFarAllCC);
00716   TFile *corrf      = new TFile("$MINOS_DATA/d84/niki/mad_cedar_daikon/near/correction_plots/far/norm_fix_far.root","READ");
00717   fcorf = (TH1F*)(corrf->Get("cor_neuf"));
00718 
00719   //find the position in the beam vector
00720   NCType::ERunType runI = NCType::kRunI;
00721   NCType::ERunType runII = runI;//NCType::kRunII; Temp! FIXME - LLH
00722   int farBeamI = NCExtrapolation::FindBeamIndex(Detector::kFar, BeamType::kL010z185i, runI);
00723   int farBeamII = NCExtrapolation::FindBeamIndex(Detector::kFar, BeamType::kL010z185i, runII);
00724   int nearBeamI = NCExtrapolation::FindBeamIndex(Detector::kNear, BeamType::kL010z185i, runI);
00725   int nearBeamII = NCExtrapolation::FindBeamIndex(Detector::kNear, BeamType::kL010z185i, runII);
00726   if(farBeamI > 100 || farBeamII > 100
00727      || nearBeamI > 100 || nearBeamII > 100){
00728     MSG("NCExtrapolationNS", Msg::kWarning) << "could not find "
00729                                             << "FD/ND beam in vector - bail"
00730                                             << " far: " << farBeamI << " " << farBeamII
00731                                             << " near: " << nearBeamI << " " << nearBeamII
00732                                             << endl;
00733     return;
00734   }
00735 
00736   //the next 2 lines just make everything look good without doing any fits
00737   //ie just the nominal mc and data
00738   //return;
00739 
00740   //get near and far histograms
00741   //datsel histograms are near data
00742   TH1F *datsel = new TH1F("datsel", "", fnbins, fbinvect);
00743   TH1F *datselpost = new TH1F("datselpost", "", fnbins, fbinvect);
00744   fdatselFD = new TH1F("datselstoreFD", "", fnbins, fbinvect); //LLH added
00745   fdatselFD_nc = new TH1F("datselstoreFD_nc", "", fnbins, fbinvect); //LLH added
00746   //TH1F *datselnc = fBeams[nearBeamI].GetDataHistogram(NCType::kNC);
00747   //TH1F *datselncpost = fBeams[nearBeamII].GetDataHistogram(NCType::kNC);
00748   FillDataHistogramFromEnergyBins(nearBeamI, NCType::kCC, datsel);
00749   FillDataHistogramFromEnergyBins(nearBeamII, NCType::kCC, datselpost);
00750   FillDataHistogramFromEnergyBins(farBeamI, NCType::kCC, fdatselFD); //LLH added
00751   FillDataHistogramFromEnergyBins(farBeamI, NCType::kNC, fdatselFD_nc); //LLH added
00752 
00753   //put far data histograms into Niki's names for the fitting
00754   FillDataHistogramFromEnergyBins(farBeamI, NCType::kCC, foscreco_rebinned);
00755   FillDataHistogramFromEnergyBins(farBeamII, NCType::kCC, foscreco_rebinned_post);
00756   FillDataHistogramFromEnergyBins(farBeamI, NCType::kNC, foscreco_ncfit_rebinned);
00757   FillDataHistogramFromEnergyBins(farBeamII, NCType::kNC, foscreco_ncfit_rebinned_post);
00758 
00759   //this method needs to know the relative scaling between near and far
00760   //exposures.  the matrix extrapolates based on the ND exposure
00761   //scale near detector exposure to be same as far detector
00762   double potNear = fBeams[nearBeamI].GetDataPOT();
00763   double potFar = fBeams[farBeamI].GetDataPOT();
00764 
00765   fNearExposureScale = potFar/potNear;
00766   cout <<"Scaling datsel!!!!!!! scale factor = " << fNearExposureScale
00767        <<"\npotFar = " << potFar
00768        <<"\npotNear = " << potNear
00769        << endl; //LLH
00770   fdatselND = (TH1F*)datsel->Clone("datselstoreND"); //LLH unmodified datsel
00771   datsel->Scale(fNearExposureScale);
00772   datselpost->Scale(fNearExposureScale);
00773 
00774   double scale_factNC = 1.;
00775   double nc_shift  = 1.;
00776   double enu_shift = 1.;
00777 
00778   // fit loop
00779   float dmsqmin     = 0.;
00780   float s2tmin      = 0.;
00781   float s2smin      = 0.;
00782   float chisqmin    = 9999.;
00783   float norm_min    = 9999.;
00784   float enu_shiftmin= 9999.;
00785   float nc_shiftmin = 9999.;
00786   int best_enu      = 9999;
00787   int best_nc       = 9999;
00788 
00789   int ndf = 39 + 39 - 3;
00790   if(fDoMC < 1) ndf=15+15-3;
00791 
00792   vector <vector <vector <double> > > chisqarraymins;
00793   chisqarraymins.resize(NCType::kNumDeltaMSqrBins);
00794   for(int i = 0; i < NCType::kNumDeltaMSqrBins; i++){
00795     chisqarraymins[i].resize(NCType::kNumUS3SqrBins);
00796     for(int j = 0; j < NCType::kNumUS3SqrBins; j++){
00797       chisqarraymins[i][j].resize(NCType::kNumUMu3SqrBins);
00798       for(int k = 0; k < NCType::kNumUMu3SqrBins; k++){
00799         chisqarraymins[i][j][k] = 1.e12;
00800       }
00801     }
00802   }
00803 
00804   float dmsqfit = 0.;
00805   float s2tfit = 0.;
00806   float s2sfit = 0.;
00807 
00808   TString math;
00809 
00810   TString noosc;
00811 
00812   noosctruearray.resize(n1);
00813   noosctruearraypost.resize(n1);
00814 
00815   //initializing arrays
00816   for(int i = 0; i < n1; ++i){
00817     math = TString::Format("truerecomat%d",i);
00818     matarray.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
00819     math = TString::Format("truerecomatnc%d",i);
00820     matarraync.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
00821     math = TString::Format("truerecomatb%d",i);
00822     matarrayb.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
00823     math = TString::Format("truerecomatncb%d",i);
00824     matarrayncb.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
00825     //TAUS
00826     math = TString::Format("truerecomatt%d",i);
00827     matarrayt.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
00828     math = TString::Format("truerecomatbt%d",i);
00829     matarraybt.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
00830 
00831     math = TString::Format("truerecomat_post%d",i);
00832     matarray_post.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
00833     math = TString::Format("truerecomatnc_post%d",i);
00834     matarraync_post.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
00835     math = TString::Format("truerecomatb_post%d",i);
00836     matarrayb_post.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
00837     math = TString::Format("truerecomatncb_post%d",i);
00838     matarrayncb_post.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
00839     //TAUS
00840     math = TString::Format("truerecomatt_post%d",i);
00841     matarrayt_post.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
00842     math = TString::Format("truerecomatbt_post%d",i);
00843     matarraybt_post.push_back(new TH2F(math,"",fnbins,fbinvect,fnbins,fbinvect));
00844 
00845     for(int j = 0; j < nc1; ++j){
00846 
00847       noosc = TString::Format("noosctrue_%d_%d",i,j);
00848       noosctruearray[i].push_back(new TH1F(noosc,"",fnbins,fbinvect));
00849       noosc = TString::Format("noosctrue_post_%d_%d",i,j);
00850       noosctruearraypost[i].push_back(new TH1F(noosc,"",fnbins,fbinvect));
00851     }
00852   }//end initialization of arrays
00853 
00854   //Getting efficiencies (Far only)
00855   GetEff(feffic, fefficnc, fefficbg, fefficncbg, fratccnc, frattau, farBeamI);
00856   GetEff(feffic_post, fefficnc_post, fefficbg_post, fefficncbg_post, fratccnc_post, frattau_post, farBeamII);
00857   if(fDoMC==0){
00858     GetEfft(feffict, fefficbgt, farBeamI);
00859     GetEfft(feffict_post, fefficbgt_post, farBeamII);
00860   }
00861 
00862   double purbin = 0.;
00863   double newbin = 0.;
00864 
00865   // loop for truereco mat efficiency and purity For Far and Near
00866   int start = n1d/2;
00867   int end = start + 1;
00868   int startnc = ncd/2;
00869   int endnc = startnc + 1;
00870 
00871   if(fSystematicsToUse[NCType::kShowerEnergy]){
00872     start   = n0;
00873     end     = n1;
00874   }
00875   if(fSystematicsToUse[NCType::kNCBackground]){
00876     startnc = nc0;
00877     endnc   = nc1;
00878   }
00879 
00880   //use TH1::Copy function to copy the ftruerecomat* matricies into the
00881   //array of matricies matarray*
00882 
00883   for(int ii = start; ii < end; ii++){
00884     // SHOWER LOOP
00885     enu_shift = les+ii*size_e;
00886     MSG("NCExtrapolationNS", Msg::kInfo) << " enu  shift " << enu_shift << endl;
00887     cout <<"First Shower Loop!" << endl;
00888 
00889     RecoTrueEffPur(matarray[ii], matarraync[ii], matarrayb[ii],
00890                    matarrayncb[ii], enu_shift, farBeamI);
00891     // here TAUS
00892     if(fDoMC==0)
00893       RecoTrueEffPurTau(matarrayt[ii], matarraybt[ii], enu_shift, farBeamI);
00894 
00895     RecoTruePur(frecotruemat, fpurityn, fefficiencyn, enu_shift, nearBeamI);
00896 
00897     for(int gg = 0; gg < fnbins; gg++){
00898       purbin = fpurityn->GetBinContent(gg+1);
00899       newbin = purbin/(purbin+scale_factNC*(1-purbin));
00900       fpurityn->SetBinContent(gg+1,newbin);
00901     }
00902 
00903     // post shutdown
00904 //     RecoTrueEffPur(matarray_post[ii], matarraync_post[ii], matarrayb_post[ii],
00905 //                 matarrayncb_post[ii], enu_shift, farBeamII);
00906 
00907 //     // here TAUS
00908 //     if(fDoMC==0)
00909 //       RecoTrueEffPurTau(matarrayt_post[ii], matarraybt_post[ii], enu_shift, farBeamII);
00910 
00911 //     RecoTruePur(frecotruemat_post, f_purityn_post, fefficiencyn_post, enu_shift, nearBeamII);
00912 
00913 //     for(int gg = 0; gg < fnbins; gg++){
00914 //       purbin = fpurityn_post->GetBinContent(gg+1);
00915 //       newbin = purbin/(purbin+scale_factNC*(1-purbin));
00916 //       fpurityn_post->SetBinContent(gg+1,newbin);
00917 //     }
00918 
00919     // NC LOOP
00920     for(int jj = startnc; jj < endnc; jj++){
00921       nc_shift = lnc+jj*size_nc;
00922       MSG("NCExtrapolationNS", Msg::kInfo) << " nc  shift " << nc_shift << endl;
00923       cout <<"NC loop within first shower loop" << endl;
00924 
00925       // MAKE PREDICTION UNOSCILLATED //
00926       for(int gg = 0; gg < fnbins; gg++){
00927         purbin = fpurityn->GetBinContent(gg+1);
00928         newbin = purbin/(purbin+nc_shift*(1-purbin));
00929         fpurityn->SetBinContent(gg+1,newbin);
00930       }
00931 
00932       //Temp - Code speedup for debugging only!
00933       Int_t bypass = 1;
00934       if(bypass == 1)
00935       {
00936          cout <<"Bypass 3!" << endl;
00937          TFile* f = new TFile("/local/scratch04/llhsu/tempFiles/L010z185i/extrapolation_fn_NS_16kND488FD_ncfvcor.root","READ");
00938          TString noosctruename = TString::Format("noosctrue_%d_%d", ii, jj);
00939          TH1F* tempHisto = (TH1F*)f->Get("plotsBeamMatrix/"+noosctruename);
00940          noosctruename += "_clone";
00941          noosctruearray[ii][jj] = (TH1F*)tempHisto->Clone(noosctruename);
00942          noosctruearray[ii][jj]->SetDirectory(0);
00943          //       fdatselND = (TH1F*)((f->Get("datselstoreND"))->Clone("datselstoreND"));
00944 //         fpurityn = (TH1F*)((f->Get("purityn"))->Clone("purityn"));
00945 //         fefficiencyn = (TH1F*)((f->Get("efficiencyn"))->Clone("efficiencyn"));
00946 //         frecotruemat = (TH2F*)((f->Get("recotruemat"))->Clone("recotruemat"));
00947          //f->Close();
00948          cout <<"Done 3!" << endl;
00949       }
00950       else
00951       {
00952 
00953       //take ND "data" and generate predicted true far spectrum
00954       MakePrediction(datsel, fpurityn, fefficiencyn,
00955                      frecotruemat, noosctruearray[ii][jj], runI);
00956 
00957      } //remove else brackets with bypass
00958 
00959       // MAKE PREDICTION UNOSCILLATED //
00960 //       for(int gg = 0; gg < fnbins; gg++){
00961 //      purbin = fpuritynp->GetBinContent(gg+1);
00962 //      newbin = purbin/(purbin+nc_shift*(1-purbin));
00963 //      fpuritynp->SetBinContent(gg+1,newbin);
00964 //       }
00965 
00966 //       MakePrediction(datselpost, fpurityn_post, fefficiencyn_post,
00967 //                   frecotruemat_post, noosctruearraypost[ii][jj], runII);
00968 
00969       // restore purity
00970       for(int gg = 0; gg < fnbins; gg++){
00971         purbin = fpurityn->GetBinContent(gg+1);
00972         newbin = (nc_shift*purbin)/(1-purbin+nc_shift*purbin);
00973         fpurityn->SetBinContent(gg+1,newbin);
00974 
00975 //      purbin = fpuritynp->GetBinContent(gg+1);
00976 //      newbin = (nc_shift*purbin)/(1-purbin+nc_shift*purbin);
00977 //      fpuritynp->SetBinContent(gg+1,newbin);
00978       }
00979     }//end nc loop
00980 
00981 //Temp LLH, for 1 iteration only!    frecotruemat->Reset();
00982 //Temp LLH, for 1 iteration only!    fpurityn->Reset();
00983 //Temp LLH, for 1 iteration only!    fefficiencyn->Reset();
00984 
00985 //     frecotruemat_post->Reset();
00986 //     fpurityn_post->Reset();
00987 //     fefficiencyn_post->Reset();
00988 
00989   } // end shower loop
00990 
00991   MSG("NCExtrapolationNS", Msg::kInfo) << " DONE EFF PUR RECO->TRUE NEAR FAR " << endl;
00992 
00993   float chisqtemp1 = 1.e12;
00994   //  float chisqtemp2 = 1.e12;
00995   float chisqtemp3 = 1.e12;
00996   //  float chisqtemp4 = 1.e12;
00997   float chisqmin_norm = 1.e12;
00998 
00999   float nc_min   = 9999.;
01000   float sh_min   = 9999.;
01001   double norm = 1.;
01002   int decayMode = 1;
01003   int decayModeNC = 0;
01004   int decayModeTau = 11;
01005   if(fDecay == 1){
01006     decayMode = -1;
01007     decayModeNC = -1;
01008     decayModeTau = -11;
01009   }
01010 
01011   MSG("NCExtrapolationNS", Msg::kInfo) << " outside loop SHOWER Start "
01012                                        << start << " END " << end
01013                                        << " NC Start " << startnc
01014                                        << " END " << endnc << endl;
01015 
01016 
01017   int maxsin = 1;//99;//NCType::kNumUMu3SqrBins;
01018   int maxdm = 1; //119;//NCType::kNumDeltaMSqrBins;
01019   int maxfs = 1;//NCType::kNumUS3SqrBins;
01020   vector<double> sinSqrDeltaMSqr;
01021 
01022   //Fitting and oscillations??
01023   // SHOWER LOOP // ii
01024   for(int ii = start; ii < end; ii++){
01025     enu_shift = les + ii*size_e;
01026     MSG("NCExtrapolationNS", Msg::kDebug) << " Shower Loop " << ii << endl;
01027     cout <<"starting second shower loop, enu_shift = " << enu_shift << endl;
01028     // NC LOOP  // jj
01029     for(int jj = startnc; jj < endnc; jj++){
01030       nc_shift = lnc + jj*size_nc;
01031       MSG("NCExtrapolationNS", Msg::kDebug) << " NC Loop " << jj
01032                                             << " nc_shift " << nc_shift << endl;
01033       cout  <<"nc loop within second shower loop, nc_shift = " << nc_shift << endl;
01034 
01035       //I need four histos here
01036       //true cc sel as cc, true cc sel as nc,
01037       //true nc sel as nc, true nc sel as cc
01038 
01039       fnooscar->Reset();
01040       fnooscar_post->Reset();
01041       fnooscar_ccselcc->Reset();
01042       fnooscar_ccselcc_post->Reset();
01043       fnooscar_ccselnc->Reset();
01044       fnooscar_ccselnc_post->Reset();
01045       fnooscar_ncselnc->Reset();
01046       fnooscar_ncselnc_post->Reset();
01047       fnooscar_ncselcc->Reset();
01048       fnooscar_ncselcc_post->Reset();
01049       // taus
01050       fnooscar_ccselcct->Reset();
01051       fnooscar_ccselcct_post->Reset();
01052       fnooscar_ccselnct->Reset();
01053       fnooscar_ccselnct_post->Reset();
01054 
01055 //Temp disabling to check fitting      fnooscar->Add(noosctruearray[ii][jj]); //Copying prediction into working histo!!
01056       fnooscar->Add(fallccmcf); //Temp using true FD spectrum (before pid efficiency corr) instead of prediction!
01057 
01058       fTrueFDAfterBM->Add(noosctruearray[ii][jj]);  //for storage
01059 //       fnooscar_post->Add(noosctruearraypost[ii][jj]);
01060 
01061       //PRE
01062       fnooscar_ccselcc->Multiply(fnooscar,feffic);  //truecc, cclike, efficiency correction
01063       fTrueFDAfterBM_posteffcorr=(TH1F*)fnooscar_ccselcc->Clone("trueFDAfterBM_posteffcorr"); //LLH adding to compare to fTrueFD after PID'as applied
01064       fnooscar_ccselnc->Multiply(fnooscar,fefficbg); //truecc, nclike, efficiency correction
01065       //Obtain true NC Spectrum from the NC/CC Ratio
01066       fnooscar_ncselnc->Multiply(fnooscar,fratccnc); //truenc, nclike, cross-section correction
01067       fnooscar_ncselcc->Multiply(fnooscar,fratccnc); //truenc, cclike, cross-section correction
01068       fnooscar_ncselnc->Multiply(fefficnc); //efficiency correction
01069       fnooscar_ncselcc->Multiply(fefficncbg); //efficiency correction
01070 
01071       // POST
01072 //       fnooscar_ccselcc_post->Multiply(fnooscar_post,fefficp);
01073 //       fnooscar_ccselnc_post->Multiply(fnooscar_post,fefficbgp);
01074 //       //Obtain true NC Spectrum from the NC/CC Ratio
01075 //       fnooscar_ncselnc_post->Multiply(fnooscar_post,fratccncp);
01076 //       fnooscar_ncselcc_post->Multiply(fnooscar_post,fratccncp);
01077 //       fnooscar_ncselnc_post->Multiply(fefficncp);
01078 //       fnooscar_ncselcc_post->Multiply(fefficncbgp);
01079 
01080       // taus
01081       if(fDoMC==0){
01082         // PRE
01083         fnooscar_ccselcct->Multiply(fnooscar,frattau);
01084         fnooscar_ccselcct->Multiply(feffict);
01085         fnooscar_ccselnct->Multiply(fnooscar,frattau);
01086         fnooscar_ccselnct->Multiply(fefficbgt);
01087 
01088         // POST
01089 //      fnooscar_ccselcct_post->Multiply(fnooscar_post,frattaup);
01090 //      fnooscar_ccselcct_post->Multiply(feffictp);
01091 //      fnooscar_ccselnct_post->Multiply(fnooscar_post,frattaup);
01092 //      fnooscar_ccselnct_post->Multiply(fefficbgtp);
01093       }
01094 
01095       //Loop over deltamsq
01096       for(int i = 0; i < maxdm; i++){
01097         dmsqfit = (1.*i + 0.5)*NCType::kDeltaDeltaMSqr+NCType::kDeltaMSqrStart;
01098         dmsqfit *= 1.e-3;
01099 
01100         MSG("NCExtrapolationNS", Msg::kInfo) << " on Delta m^2 = " << dmsqfit
01101                                              << " chi^2 min = " << chisqmin << endl;
01102         cout <<"Inside deltamsq loop! dmsqfit = " << dmsqfit << endl;
01103 
01104         sinSqrDeltaMSqr.clear();
01105         //really just need the value of the sin^2(1.27*\Delta m^2 L/E) factor
01106         //at the middle of each bin as the bins are of size 0.25 GeV in the
01107         //region of oscillation interest - that is good enough to get the oscillated
01108         //true value to within 0.0004%
01109         for(int dmb = 0; dmb < fnbins; ++dmb){  //Why was this set to 27?? - LLH
01110           sinSqrDeltaMSqr.push_back(TMath::Power(TMath::Sin(NCType::k127*dmsqfit
01111                                                             *NCType::kBaseLineFar/fnooscreco->GetBinCenter(dmb+1)),2.));
01112         }//end loop over energy bins
01113 
01114         for(int j = 0; j < maxsin; j++){
01115           MSG("NCExtrapolationNS", Msg::kDebug) << " In sin loop " << j <<  endl;
01116           s2tfit = (1.*j)*NCType::kDeltaUMu3Sqr+NCType::kUMu3SqrStart; //Temp! (1.*j + 0.5)*NCType::kDeltaUMu3Sqr+NCType::kUMu3SqrStart;
01117           cout <<"Inside sinsq loop! s2tfit = " << s2tfit << endl;
01118 
01119           for(int k = 0; k < maxfs; k++){
01120             s2sfit = (1.*k + 0.5)*NCType::kDeltaUS3Sqr+NCType::kUS3SqrStart;
01121             cout <<"Inside fsterile loop! s2sfit = " << s2sfit << endl;
01122 
01123             //PRE SHUTDOWN
01124 
01125             //Oscillate true cc's selected as cc's after efficiency and turn them to reco
01126             // LAST 1 IS FOR CC
01127             //oscweight(fnooscar_ccselcc,*ftemptrue,dmsqfit,s2tfit,1,decayMode);
01128             oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselcc,ftemptrue,14,NCType::kCC);
01129             fTrueFDAfterOsc = (TH1F*)ftemptrue->Clone("trueFDAfterOsc"); //should be identical to fTrueFDAfterBM_posteffcorr w/o osc
01130             //truetoreco(fTrueFD,*ftempreco,matarray[ii]);  //Temp! To debug unsmearing - LLH
01131             truetoreco(ftemptrue,*ftempreco,matarray[ii]);
01132             fRecoFDAfterBM->Add(ftempreco);
01133             ftemptrue->Reset();
01134             //Oscillate true nc's selected as cc's after efficiency and turn them to reco
01135             //oscweight (fnooscar_ncselcc,*ftemptrue,dmsqfit,s2sfit,1,decayModeNC);
01136             oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ncselcc,ftemptrue,14,NCType::kNC);
01137             truetoreco(ftemptrue,*ftempreconc,matarrayncb[ii]);
01138             ftempreconc->Scale(nc_shift);
01139             ftemptrue->Reset();
01140             //Oscillate true nc's selected as nc's after efficiency and turn them to reco
01141             //oscweight(fnooscar_ncselnc,*ftemptrue,dmsqfit,s2sfit,1,decayModeNC);
01142             oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ncselnc,ftemptrue,14,NCType::kNC);
01143             truetoreco(ftemptrue,*ftempreco_ncfit,matarraync[ii]);
01144             ftempreco_ncfit->Scale(nc_shift);
01145             ftemptrue->Reset();
01146             //Oscillate true cc's selected as nc's after efficiency and turn them to reco
01147             //oscweight(fnooscar_ccselnc,*ftemptrue,dmsqfit,s2tfit,1,decayMode);
01148             oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselnc,ftemptrue,14,NCType::kCC);
01149             truetoreco(ftemptrue,*ftemprecocc_ncfit,matarrayb[ii]);
01150             ftemptrue->Reset();
01151 
01152             // TAUS  HERE
01153             if(fDoMC==0){
01154               //Oscillate true cc's selected as cc's after efficiency and turn them to reco
01155               // LAST 11 IS FOR CC TAU
01156               //oscweight(fnooscar_ccselcct,*ftemptrue,dmsqfit,s2tfit,1,decayModeTau);
01157               oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselcct,ftemptrue,16,NCType::kCC);
01158               truetoreco(ftemptrue,*ftemprecotau,matarrayt[ii]);
01159               ftemptrue->Reset();
01160 
01161               //Oscillate true cc's selected as nc's after efficiency and turn them to reco
01162               //oscweight(fnooscar_ccselnct,*ftemptrue,dmsqfit,s2tfit,1,decayModeTau);
01163               oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselnct,ftemptrue,16,NCType::kNC);
01164               truetoreco(ftemptrue,*ftemprecotau_ncfit,matarraybt[ii]);
01165               ftemptrue->Reset();
01166             }   //end taus
01167 
01168 //          //POST SHUTDOWN
01169 
01170 //          //Oscillate true cc's selected as cc's after efficiency and turn them to reco
01171 //          // LAST 1 IS FOR CC
01172 //          //oscweight(fnooscar_ccselcc_post,*ftemptrue,dmsqfit,s2tfit,1,decayMode);
01173 //          oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselcc_post,ftemptrue,14,NCType::kCC);
01174 //          truetoreco(ftemptrue,*ftempreco_post,matarray_post[ii]);
01175 //          ftemptrue->Reset();
01176 //          //Oscillate true nc's selected as cc's after efficiency and turn them to reco
01177 //          //oscweight(fnooscar_ncselcc_post,*ftemptrue,dmsqfit,s2sfit,1,decayModeNC);
01178 //          oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ncselcc_post,ftemptrue,14,NCType::kNC);
01179 //          truetoreco(ftemptrue,*ftempreconc_post,matarrayncb_post[ii]);
01180 //          ftempreconc_post->Scale(nc_shift);
01181 //          ftemptrue->Reset();
01182 
01183 //          //Oscillate true nc's selected as nc's after efficiency and turn them to reco
01184 //          //oscweight(fnooscar_ncselnc_post,*ftemptrue,dmsqfit,s2sfit,1,decayModeNC);
01185 //          oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ncselnc_post,ftemptrue,14,NCType::kNC);
01186 //          truetoreco(ftemptrue,*ftempreco_ncfit_post,matarraync_post[ii]);
01187 //          ftempreco_ncfit_post->Scale(nc_shift);
01188 //          ftemptrue->Reset();
01189 //          //Oscillate true cc's selected as nc's after efficiency and turn them to reco
01190 //          //oscweight(fnooscar_ccselnc_post,*ftemptrue,dmsqfit,s2tfit,1,decayMode);
01191 //          oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselnc_post,ftemptrue,14,NCType::kCC);
01192 //          truetoreco(ftemptrue,*ftemprecocc_ncfit_post,matarrayb_post[ii]);
01193 //          ftemptrue->Reset();
01194 
01195 //          // TAUS  HERE
01196 //          if(fDoMC==0){
01197 //            //Oscillate true cc's selected as cc's after efficiency and turn them to reco
01198 //            // LAST 11 IS FOR CC TAU
01199 //            //oscweight(fnooscar_ccselcct_post,*ftemptrue,dmsqfit,s2tfit,1,decayModeTau);
01200 //            oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselcct_post,ftemptrue,16,NCType::kCC);
01201 //            truetoreco(ftemptrue,*ftemprecotau_post,matarrayt_post[ii]);
01202 //            ftemptrue->Reset();
01203 
01204 //            //Oscillate true cc's selected as nc's after efficiency and turn them to reco
01205 //            //oscweight(fnooscar_ccselnct_post,*ftemptrue,dmsqfit,s2tfit,1,decayModeTau);
01206 //            oscweight(sinSqrDeltaMSqr,s2tfit,s2sfit,fnooscar_ccselnct_post,ftemptrue,16,NCType::kCC);
01207 //            truetoreco(ftemptrue,*ftemprecotau_ncfit_post,matarraybt_post[ii]);
01208 //            ftemptrue->Reset();
01209 //          }
01210 
01211             chisqmin_norm = 1.e12;
01212             // NORMALIZATION // kk
01213             for(int kk = nr0; kk < nr1; kk++){
01214               norm = lnr + kk*size_nr;
01215               // CC
01216               ftempreconc->Scale(norm);
01217               ftempreco->Scale(norm);
01218 //            ftempreconc_post->Scale(norm);
01219 //            ftempreco_post->Scale(norm);
01220               if(fDoMC==0){
01221                 ftemprecotau->Scale(norm);
01222 //              ftemprecotau_post->Scale(norm);
01223               }
01224               // NC
01225               ftemprecocc_ncfit->Scale(norm);
01226               ftempreco_ncfit->Scale(norm);
01227 //            ftemprecocc_ncfit_post->Scale(norm);
01228 //            ftempreco_ncfit_post->Scale(norm);
01229               if(fDoMC==0){
01230                 ftemprecotau_ncfit->Scale(norm);
01231 //              ftemprecotau_ncfit_post->Scale(norm);
01232               }
01233               //ADD THEM ALL CC
01234               ftempreco->Add(ftempreconc);
01235 //            ftempreco_post->Add(ftempreconcp);
01236               //ADD THEM ALL NC
01237               ftempreco_ncfit->Add(ftemprecocc_ncfit);
01238 //            ftempreco_ncfit_post->Add(ftemprecocc_ncfitp);
01239               if(fDoMC == 0){
01240                 ftempreco->Add(ftemprecotau);
01241 //              ftempreco_post->Add(ftemprecotaup);
01242                 ftempreco_ncfit->Add(ftemprecotau_ncfit);
01243 //              ftempreco_ncfit_post->Add(ftemprecotau_ncfitp);
01244               }
01245 
01246               //REBIN
01247               Rebin(ftempreco, *ftempreco_rebinned);
01248               Rebin(ftempreconc, *ftempreconc_rebinned);
01249 //            Rebin(ftempreco_post, *ftempreco_rebinned_post);
01250 //            Rebin(ftempreconc_post, *ftempreconc_rebinned_post);
01251               Rebin(ftempreco_ncfit, *ftempreco_ncfit_rebinned);
01252               Rebin(ftemprecocc_ncfit, *ftemprecocc_ncfit_rebinned);
01253 //            Rebin(ftempreco_ncfit_post, *ftempreco_ncfit_rebinned_post);
01254 //            Rebin(ftemprecocc_ncfit_post, *ftemprecocc_ncfit_rebinned_post);
01255 
01256               // CHISQUARE
01257               chisqtemp1 = chisq(foscreco_rebinned, ftempreco_rebinned, ftempreconc_rebinned);
01258               chisqtemp1 += ( ((1.-norm)/signorm)*((1.-norm)/signorm)
01259                               +((1.-enu_shift)/siges)*((1.-enu_shift)/siges)
01260                               +((1.-nc_shift)/signc)*((1.-nc_shift)/signc) );
01261 //            chisqtemp2 = chisq(foscreco_rebinned_post, ftempreco_rebinned_post, ftempreconc_rebinned_post);
01262 //            chisqtemp2 += ( ((1.-norm)/signorm)*((1.-norm)/signorm)
01263 //                            +((1.-enu_shift)/siges)*((1.-enu_shift)/siges)
01264 //                            +((1.-nc_shift)/signc)*((1.-nc_shift)/signc) );
01265               chisqtemp3 = chisq(foscreco_ncfit_rebinned, ftempreco_ncfit_rebinned, ftemprecocc_ncfit_rebinned);
01266 //            chisqtemp3 += ( ((1.-norm)/signorm)*((1.-norm)/signorm)
01267 //                            +((1.-enu_shift)/siges)*((1.-enu_shift)/siges)
01268 //                            +((1.-nc_shift)/signc)*((1.-nc_shift)/signc) );
01269 //            chisqtemp4 = chisq(foscreco_ncfit_rebinned_post, ftempreco_ncfit_rebinned_post, ftemprecocc_ncfit_rebinned_post);
01270 //            chisqtemp4 += ( ((1.-norm)/signorm)*((1.-norm)/signorm)
01271 //                            +((1.-enu_shift)/siges)*((1.-enu_shift)/siges)
01272 //                            +((1.-nc_shift)/signc)*((1.-nc_shift)/signc) );
01273 
01274               chisqtemp1 += chisqtemp3;//chisqtemp2 + chisqtemp3 + chisqtemp4;
01275 
01276               if(chisqtemp1 < chisqmin_norm) chisqmin_norm = chisqtemp1;
01277 
01278               if(chisqtemp1 < chisqmin){
01279                 chisqmin = chisqtemp1;
01280                 dmsqmin  = dmsqfit;
01281                 s2tmin   = s2tfit;
01282                 s2smin   = s2sfit;
01283                 norm_min = norm;
01284                 nc_min   = nc_shift;
01285                 sh_min   = enu_shift;
01286                 best_enu = ii;
01287                 best_nc  = jj;
01288               }
01289               // Un-Normalize
01290               if(fDoMC == 0){
01291                 ftempreco->Add(ftemprecotau,-1);
01292                 ftemprecotau->Scale(1./norm);
01293 //              ftempreco_post->Add(ftemprecotau_post,-1);
01294 //              ftemprecotau_post->Scale(1./norm);
01295                 ftempreco_ncfit->Add(ftemprecotau_ncfit,-1);
01296                 ftemprecotau_ncfit->Scale(1./norm);
01297 //              ftempreco_ncfit_post->Add(ftemprecotau_ncfit_post,-1);
01298 //              ftemprecotau_ncfit_post->Scale(1./norm);
01299               }
01300               ftempreco->Add(ftempreconc,-1);
01301               ftempreco->Scale(1./norm);
01302               ftempreconc->Scale(1./norm);
01303 //            ftempreco_post->Add(ftempreconc_post,-1);
01304 //            ftempreco_post->Scale(1./norm);
01305 //            ftempreconc_post->Scale(1./norm);
01306               ftempreco_ncfit->Add(ftemprecocc_ncfit,-1);
01307               ftempreco_ncfit->Scale(1./norm);
01308               ftemprecocc_ncfit->Scale(1./norm);
01309 //            ftempreco_ncfit_post->Add(ftemprecocc_ncfit_post,-1);
01310 //            ftempreco_ncfit_post->Scale(1./norm);
01311 //            ftemprecocc_ncfit_post->Scale(1./norm);
01312 
01313             }// NORMALIZATION
01314 
01315             chisqarraymins[i][k][j] = chisqmin_norm;
01316 
01317             MSG("NCExtrapolationNS", Msg::kInfo) << " (dm^2, sin^2 2theta, fsterile) = ("
01318                                                  << dmsqfit << "," << s2tfit << ","
01319                                                  << s2sfit << ") chi^2 = "
01320                                                  << chisqmin_norm << endl;
01321           }// fs
01322         }// sin2theta
01323       }// dm2
01324 
01325       //Unnormalized NC
01326       ftempreconc->Scale(1./nc_shift);
01327       ftempreco_ncfit->Scale(1./nc_shift);
01328 //       ftempreconc_post->Scale(1./nc_shift);
01329 //       ftempreco_ncfit_post->Scale(1./nc_shift);
01330 
01331     }// NC LOOP
01332   }// SHOWER LOOP
01333 
01335 
01336   //set the best fit values
01337   BestUMu3Sqr() = s2tmin;
01338   BestUS3Sqr() = s2smin;
01339   BestDeltaMSqr() = dmsqmin;
01340   fMinChiSqr = chisqmin;
01341 
01342   cout <<"Bestfit values stored as Oscillation Parameters!: "
01343        <<"\nUMu3SqrTheta = " << BestUMu3Sqr()
01344        <<"\nUS3Sqr = " << BestUS3Sqr()
01345        <<"\nDeltaMSqr = " << BestDeltaMSqr()
01346        <<"\nMinChiSqr = " << fMinChiSqr
01347        << endl;
01348 
01349 //   fBeams[farBeamI].SetOscillationParameters(fBestUMu3Sqr, fBestDeltaMSqr, fBestUS3Sqr);
01350 //   fBeams[farBeamII].SetOscillationParameters(fBestUMu3Sqr, fBestDeltaMSqr, fBestUS3Sqr);
01351 
01352   NCExtrapolation::FillDeltaChiSqrMaps(chisqarraymins);
01353 
01354   // Now fill chisq
01355   for(int i = 0; i < maxdm; i++){// dm2
01356     dmsqfit = (1.*i + 0.5)*NCType::kDeltaDeltaMSqr+NCType::kDeltaMSqrStart;
01357 
01358     for(int j = 0; j < maxsin; j++){ //sin2theta
01359       s2tfit = (1.*j + 0.5)*NCType::kDeltaUMu3Sqr+NCType::kUMu3SqrStart;
01360 
01361       double chismin =1e9;
01362       int    kmin    =-1;
01363       for(int k = 0; k < maxfs; k++){ //fsterile
01364         s2sfit = (1.*k + 0.5)*NCType::kDeltaUS3Sqr+NCType::kUS3SqrStart;
01365 
01366         if(chisqarraymins[i][k][j]<chismin){
01367           chismin = chisqarraymins[i][k][j];
01368           kmin = k;
01369         }
01370       } // fs
01371       fchisqall->Fill(s2tfit, dmsqfit, chismin);
01372     } // sin2theta
01373   } // dm2
01374 
01375   double chismin = 1e9;
01376   int    jmin    = -1;
01377 
01378   // Now fill chisq other projection
01379   for(int i = 0; i < maxdm; i++){// dm2
01380     dmsqfit = (1.*i + 0.5)*NCType::kDeltaDeltaMSqr+NCType::kDeltaMSqrStart;
01381 
01382     for(int k = 0; k < maxfs; k++){ //fsterile
01383       s2sfit = (1.*k + 0.5)*NCType::kDeltaUS3Sqr+NCType::kUS3SqrStart;
01384 
01385       chismin = 1e9;
01386       jmin    = -1;
01387 
01388       for(int j = 0; j < maxsin; j++){ //sin2theta
01389         s2tfit = (1.*j + 0.5)*NCType::kDeltaUMu3Sqr+NCType::kUMu3SqrStart;
01390 
01391         if(chisqarraymins[i][k][j] < chismin){
01392           chismin = chisqarraymins[i][k][j];
01393           jmin = j;
01394         }
01395       } // sin2theta
01396       fchisqallnc->Fill(s2sfit, dmsqfit, chismin);
01397     } // fs
01398   } // dm2
01399 
01400   TString name = "";
01401   TString junk = "";
01402   TString domch = "data";
01403   TString dodecayh = "osc";
01404 
01405   if(fDoMC==1) domch = "mc";
01406   if(fDecay==1) dodecayh = "dec";
01407 
01408   int sele = n1d/2;
01409   int selnc = ncd/2;
01410   enu_shiftmin = 1.;
01411   nc_shiftmin  = 1.;
01412   if(fSystematicsToUse[NCType::kShowerEnergy]){
01413     sele = best_enu;
01414     enu_shiftmin = les+best_enu*size_e;
01415   }
01416   if(fSystematicsToUse[NCType::kNCBackground]){
01417     selnc = best_nc;
01418     nc_shiftmin = lnc+best_nc *size_nc;
01419   }
01420 
01421   if(fDoMC==1){
01422     name = "res_fit_";
01423     name += "_" + junk + "_" + domch + "_" + dodecayh;
01424     name += "new.txt";
01425   }
01426   if(fDoMC==0){
01427     name = "res_fit_" + domch + "_"  + dodecayh;
01428     name += "_test_post_part_new.txt";
01429   }
01430 
01431   //Temp! - removal of systematic loops, setting normalizations to 1
01432   norm_min = 1.0;
01433   enu_shiftmin = 1.0;
01434   nc_shiftmin = 1.0;
01435 
01436   ofstream text_file(name);
01437   text_file << " Best fit at dmsq= " << dmsqmin
01438             << " sin2theta= "       << s2tmin
01439             << " fs " << s2smin
01440             << "  chisq=" << chisqmin/ndf
01441             << " norm_min " << norm_min
01442             << " enu_shiftmin " << enu_shiftmin
01443             << " nc_shiftmin " << nc_shiftmin << endl;
01444   text_file.close();
01445 
01446   MSG("NCExtrapolationNS", Msg::kInfo) << " Best fit at dmsq = " << dmsqmin
01447                                        << " sin2theta = " << s2tmin
01448                                        << " fs = " << s2smin << "  chisq = "
01449                                        << chisqmin << " ndf = "  << ndf << " norm_min "
01450                                        << norm_min << " enu_shiftmin "
01451                                        << enu_shiftmin << "nc_shiftmin " << nc_shiftmin
01452                                        <<  endl;
01453   MSG("NCExtrapolationNS", Msg::kInfo) << " Best enu integer= " << best_enu
01454                                        << " best nc integer "  << best_nc << endl;
01455 
01456   // Make best fit histo//
01457   // For that we need best shower energy, best nc value,
01458   //best normalization value and the best fit points//
01459   fnooscar->Reset();
01460 //   fnooscar_post->Reset();
01461   fnooscar_ccselcc->Reset();
01462 //   fnooscar_ccselcc_post->Reset();
01463   fnooscar_ccselnc->Reset();
01464 //   fnooscar_ccselnc_post->Reset();
01465   fnooscar_ncselnc->Reset();
01466 //   fnooscar_ncselnc_post->Reset();
01467   fnooscar_ncselcc->Reset();
01468 //   fnooscar_ncselcc_post->Reset();
01469   // taus
01470   fnooscar_ccselcct->Reset();
01471 //   fnooscar_ccselcct_post->Reset();
01472   fnooscar_ccselnct->Reset();
01473 //   fnooscar_ccselnct_post->Reset();
01474 
01475   //Copying noosctruearray for best showerenergy and nc values
01476 //Temp disabling to check fitting, no looping over showerenergy and nc values at the moment   noosctruearray[sele][selnc]->Copy(*fnooscar);
01477   fallccmcf->Copy(*fnooscar); //temp using true FD spectrum (before pid efficiency corr) instead of prediction!
01478 //   noosctruearraypost[sele][selnc]->Copy(*fnooscar_post);
01479 
01480   // PRE
01481   fnooscar_ccselcc->Multiply(fnooscar,feffic);
01482   fnooscar_ccselnc->Multiply(fnooscar,fefficbg);
01483   fnooscar_ncselnc->Multiply(fnooscar,fratccnc);
01484   fnooscar_ncselcc->Multiply(fnooscar,fratccnc);
01485   fnooscar_ncselnc->Multiply(fefficnc);
01486   fnooscar_ncselcc->Multiply(fefficncbg);
01487 
01488   // POST
01489 //   fnooscar_ccselcc_post->Multiply(fnooscar_post,fefficp);
01490 //   fnooscar_ccselnc_post->Multiply(fnooscar_post,fefficbgp);
01491 //   fnooscar_ncselnc_post->Multiply(fnooscar_post,fratccncp);
01492 //   fnooscar_ncselcc_post->Multiply(fnooscar_post,fratccncp);
01493 //   fnooscar_ncselnc_post->Multiply(fefficncp);
01494 //   fnooscar_ncselcc_post->Multiply(fefficncbgp);
01495 
01496   // taus
01497   if(fDoMC==0){
01498     // PRE
01499     fnooscar_ccselcct->Multiply(fnooscar,frattau);
01500     fnooscar_ccselcct->Multiply(feffict);
01501     fnooscar_ccselnct->Multiply(fnooscar,frattau);
01502     fnooscar_ccselnct ->Multiply(fefficbgt);
01503 
01504     // POST
01505 //     fnooscar_ccselcct_post->Multiply(fnooscar_post,frattaup);
01506 //     fnooscar_ccselcct_post->Multiply(feffictp);
01507 //     fnooscar_ccselnct_post->Multiply(fnooscar_post,frattaup);
01508 //     fnooscar_ccselnct_post->Multiply(fefficbgtp);
01509 
01510   }
01511 
01512   //PRE SHUTDOWN
01513 
01514   for(int dmb = 0; dmb < fnbins; ++dmb){
01515     sinSqrDeltaMSqr.push_back(TMath::Power(TMath::Sin(NCType::k127*dmsqfit
01516                                                       *NCType::kBaseLineFar/fnooscreco->GetBinCenter(dmb+1)),2.));
01517   }//end loop over energy bins
01518 
01519   //Oscillate true cc's selected as cc's after efficiency and turn them to reco
01520   //oscweight(fnooscar_ccselcc,*ftemptrue,dmsqmin,s2tmin,1,decayMode);
01521   oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselcc,ftemptrue,14,NCType::kCC);
01522   truetoreco(ftemptrue,*fbestfitreco,matarray[sele]);
01523   ftemptrue->Reset();
01524   //Oscillate true nc's selected as cc's after efficiency and turn them to reco
01525   //oscweight(fnooscar_ncselcc,*ftemptrue,dmsqmin,s2smin,1,decayModeNC);
01526   oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ncselcc,ftemptrue,14,NCType::kNC);
01527   truetoreco(ftemptrue,*fbestfitreconc,matarrayncb[sele]);
01528   fbestfitreconc->Scale(nc_shiftmin);
01529   ftemptrue->Reset();
01530   //Oscillate true nc's selected as nc's after efficiency and turn them to reco
01531   //oscweight(fnooscar_ncselnc,*ftemptrue,dmsqmin,s2smin,1,decayModeNC);
01532   oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ncselnc,ftemptrue,14,NCType::kNC);
01533   truetoreco(ftemptrue,*fbestfitreco_ncfit,matarraync[sele]);
01534   fbestfitreco_ncfit->Scale(nc_shiftmin);
01535   ftemptrue->Reset();
01536   //Oscillate true cc's selected as nc's after efficiency and turn them to reco
01537   //oscweight(fnooscar_ccselnc,*ftemptrue,dmsqmin,s2tmin,1,decayMode);
01538   oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselnc,ftemptrue,14,NCType::kCC);
01539   truetoreco(ftemptrue,*fbestfitrecocc_ncfit,matarrayb[sele]);
01540   ftemptrue->Reset();
01541 
01542   // TAUS  HERE
01543   if(fDoMC==0){
01544     //Oscillate true cc's selected as cc's after efficiency and turn them to reco
01545     // LAST 11 IS FOR CC TAU
01546     //oscweight(fnooscar_ccselcct,*ftemptrue,dmsqmin,s2tmin,1,decayModeTau);
01547     oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselcct,ftemptrue,16,NCType::kCC);
01548     truetoreco(ftemptrue,*fbestfitrecotau,matarrayt[sele]);
01549     ftemptrue->Reset();
01550 
01551     //Oscillate true cc's selected as nc's after efficiency and turn them to reco
01552     //oscweight(fnooscar_ccselnct,*ftemptrue,dmsqmin,s2tmin,1,decayModeTau);
01553     oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselnct,ftemptrue,16,NCType::kCC);
01554     truetoreco(ftemptrue,*fbestfitrecotau_ncfit,matarraybt[sele]);
01555     ftemptrue->Reset();
01556   }
01557 
01558 //   //POST SHUTDOWN
01559 //   //Oscillate true cc's selected as cc's after efficiency and turn them to reco
01560 //   //oscweight(fnooscar_ccselcc_post,*ftemptrue,dmsqmin,s2tmin,1,decayMode);
01561 //   oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselcc_post,ftemptrue,14,NCType::kCC);
01562 //   truetoreco(ftemptrue,*fbestfitreco_post,matarray_post[sele]);
01563 //   ftemptrue->Reset();
01564 //   //Oscillate true nc's selected as cc's after efficiency and turn them to reco
01565 //   //oscweight(fnooscar_ncselcc_post,*ftemptrue,dmsqmin,s2smin,1,decayModeNC);
01566 //   oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ncselcc_post,ftemptrue,14,NCType::kNC);
01567 //   truetoreco(ftemptrue,*fbestfitreconc_post,matarrayncb_post[sele]);
01568 //   fbestfitreconc_post->Scale(nc_shiftmin);
01569 //   ftemptrue->Reset();
01570 //   //Oscillate true nc's selected as nc's after efficiency and turn them to reco
01571 //   //oscweight(fnooscar_ncselnc_post,*ftemptrue,dmsqmin,s2smin,1,decayModeNC);
01572 //   oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ncselnc_post,ftemptrue,14,NCType::kNC);
01573 //   truetoreco(ftemptrue,*fbestfitreco_ncfit_post,matarraync_post[sele]);
01574 //   fbestfitreco_ncfit_post->Scale(nc_shiftmin);
01575 //   ftemptrue->Reset();
01576 //   //Oscillate true cc's selected as nc's after efficiency and turn them to reco
01577 //   //oscweight(fnooscar_ccselnc_post,*ftemptrue,dmsqmin,s2tmin,1,decayMode);
01578 //   oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselnc_post,ftemptrue,14,NCType::kCC);
01579 //   truetoreco(ftemptrue,*fbestfitrecocc_ncfit_post,matarrayb_post[sele]);
01580 //   ftemptrue->Reset();
01581 
01582 //   // TAUS  HERE
01583 //   if(fDoMC==0){
01584 //     //Oscillate true cc's selected as cc's after efficiency and turn them to reco
01585 //     //oscweight(fnooscar_ccselcct_post,*ftemptrue,dmsqmin,s2tmin,1,decayModeTau);
01586 //     oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselcct_post,ftemptrue,16,NCType::kCC);
01587 //     truetoreco(ftemptrue,*fbestfitrecotau_post,matarrayt_post[sele]);
01588 //     ftemptrue->Reset();
01589 
01590 //     //Oscillate true cc's selected as nc's after efficiency and turn them to reco
01591 //     //oscweight(fnooscar_ccselnct_post,*ftemptrue,dmsqmin,s2tmin,1,decayModeTau);
01592 //     oscweight(sinSqrDeltaMSqr,s2tmin,s2smin,fnooscar_ccselnct_post,ftemptrue,16,NCType::kCC);
01593 //     truetoreco(ftemptrue,*fbestfitrecotau_ncfit_post,matarraybt_post[sele]);
01594 //     ftemptrue->Reset();
01595 //   }
01596 
01597   // CC
01598   fbestfitreconc->Scale(norm_min);
01599   fbestfitreco->Scale(norm_min);
01600 //   fbestfitreconc_post->Scale(norm_min);
01601 //   fbestfitreco_post->Scale(norm_min);
01602   fbestfitrecocc_ncfit->Scale(norm_min);
01603   fbestfitreco_ncfit->Scale(norm_min);
01604   fbestfitrecocc_ncfit->Scale(norm_min);
01605 //   fbestfitreco_ncfit_post->Scale(norm_min);
01606   fbestfitreco->Add(fbestfitreconc);
01607 //   fbestfitreco_post->Add(fbestfitreconcp);
01608   fbestfitreco_ncfit->Add(fbestfitrecocc_ncfit);
01609 //   fbestfitreco_ncfit_post->Add(fbestfitrecocc_ncfitp);
01610   if(fDoMC==0){
01611     fbestfitrecotau->Scale(norm_min);
01612 //     fbestfitrecotau_post->Scale(norm_min);
01613     fbestfitrecotau_ncfit->Scale(norm_min);
01614 //     fbestfitrecotau_ncfit_post->Scale(norm_min);
01615     fbestfitreco->Add(fbestfitrecotau);
01616 //     fbestfitreco_post->Add(fbestfitrecotaup);
01617     fbestfitreco_ncfit->Add(fbestfitrecotau_ncfit);
01618 //     fbestfitreco_ncfit_post->Add(fbestfitrecotau_ncfitp);
01619   }
01620 
01621   //REBIN
01622   Rebin(fbestfitreco, *fbestfitreco_rebinned);
01623   Rebin(fbestfitreconc, *fbestfitreconc_rebinned);
01624 //   Rebin(fbestfitreco_post, *fbestfitreco_rebinned_post);
01625 //   Rebin(fbestfitreconc_post, *fbestfitreconc_rebinned_post);
01626   Rebin(fbestfitreco_ncfit, *fbestfitreco_ncfit_rebinned);
01627   Rebin(fbestfitrecocc_ncfit, *fbestfitrecocc_ncfit_rebinned);
01628 //   Rebin(fbestfitreco_ncfit_post, *fbestfitreco_ncfit_rebinned_post);
01629 //   Rebin(fbestfitrecocc_ncfit_post, *fbestfitrecocc_ncfit_rebinned_post);
01630   if(fDoMC==0){
01631     Rebin(fbestfitrecotau, *fbestfitrecotau_rebinned);
01632     Rebin(fbestfitrecotau_ncfit, *fbestfitrecotau_ncfit_rebinned);
01633 //     Rebin(fbestfitrecotau_post, *fbestfitrecotau_rebinned_post);
01634 //     Rebin(fbestfitrecotau_ncfit_post, *fbestfitrecotau_ncfit_rebinned_post);
01635   }
01636 
01637   // PREDICTION UNOSCILLATED //
01638   fnooscar->Copy(*fnoosctrue);
01639 //   fnooscar_post->Copy(*fnoosctruep);
01640   fnoosctrue_ncfit->Multiply(fnooscar,fratccnc);
01641 //   fnoosctrue_ncfit_post->Multiply(fnooscar_post,fratccncp);
01642 
01643   //PRE SHUTDOWN
01644   truetoreco(fnooscar_ccselcc,*fnooscreco,  matarray[sele]);
01645   truetoreco(fnooscar_ncselcc,*fnooscreconc,matarrayncb[sele]);
01646   fnooscreconc->Scale(nc_shiftmin);
01647   truetoreco(fnooscar_ncselnc,*fnooscreco_ncfit,matarraync[sele]);
01648   fnooscreconc->Scale(nc_shiftmin);
01649   truetoreco(fnooscar_ccselnc,*fnooscrecocc_ncfit, matarrayb[sele]);
01650 
01651   //POST SHUTDOWN
01652   truetoreco(fnooscar_ccselcc_post,*fnooscreco_post,matarray_post[sele]);
01653   truetoreco(fnooscar_ncselcc_post,*fnooscreconc_post,matarrayncb_post[sele]);
01654 //   fnooscreconc_post->Scale(nc_shiftmin);
01655   truetoreco(fnooscar_ncselnc_post,*fnooscreco_ncfit_post,matarraync_post[sele]);
01656 //   fnooscreconc_post->Scale(nc_shiftmin);
01657   truetoreco(fnooscar_ccselnc_post,*fnooscrecocc_ncfit_post,matarrayb_post[sele]);
01658 
01659   // CC
01660   fnooscreconc->Scale(norm_min);
01661   fnooscreco->Scale(norm_min);
01662 //   fnooscreconc_post->Scale(norm_min);
01663   fnooscreco_post->Scale(norm_min);
01664   fnooscrecocc_ncfit->Scale(norm_min);
01665   fnooscreco_ncfit->Scale(norm_min);
01666 //   fnooscrecocc_ncfit_post->Scale(norm_min);
01667 //   fnooscreco_ncfit_post->Scale(norm_min);
01668   fnooscreco->Add(fnooscreconc);
01669 //   fnooscreco_post->Add(fnooscreconcp);
01670   fnooscreco_ncfit->Add(fnooscrecocc_ncfit);
01671 //   fnooscreconc_post->Add(fnooscrecocc_ncfitp);
01672   if(fDoMC==0){
01673     fnooscrecotau->Scale(norm_min);
01674 //     fnooscrecotau_post->Scale(norm_min);
01675     fnooscrecotau_ncfit->Scale(norm_min);
01676 //     fnooscrecotau_ncfit_post->Scale(norm_min);
01677     fnooscreco->Add(fnooscrecotau);
01678 //     fnooscreco_post->Add(fnooscrecotaup);
01679     fnooscreco_ncfit->Add(fnooscrecotau_ncfit);
01680 //     fnooscreconc_post->Add(fnooscrecotau_ncfitp);
01681   }
01682 
01683   //REBIN
01684   Rebin(fnooscreco,          *fnooscreco_rebinned);
01685   Rebin(fnooscreconc,        *fnooscreconc_rebinned);
01686 //   Rebin(fnooscreco_post,         *fnooscreco_rebinned_post);
01687 //   Rebin(fnooscreconc_post,       *fnooscreconc_rebinned_post);
01688   Rebin(fnooscreco_ncfit,       *fnooscreco_ncfit_rebinned);
01689   Rebin(fnooscrecocc_ncfit,        *fnooscrecocc_ncfit_rebinned);
01690 //   Rebin(fnooscreco_ncfit_post,      *fnooscreco_ncfit_rebinned_post);
01691 //   Rebin(fnooscrecocc_ncfit_post,       *fnooscrecocc_ncfit_rebinned_post);
01692   if(fDoMC==0){
01693     Rebin(fnooscrecotau, *fnooscrecotau_rebinned);
01694     Rebin(fnooscrecotau_ncfit, *fnooscrecotau_ncfit_rebinned);
01695 //     Rebin(fnooscrecotau_post, *fnooscrecotau_rebinned_post);
01696 //     Rebin(fnooscrecotau_ncfit_post, *fnooscrecotau_ncfit_rebinned_post);
01697   }
01698 
01699 //   //clear all the far beam MC histograms and then add the
01700 //   //results from the extrapolations to the appropriate ones
01701 // LLH uncommented
01702   int farBeams[] = {farBeamI, farBeamII};
01703   for(int b = 0; b < 2; ++b){
01704      for(int i = 0; i < 2; ++i){
01705         fBeams[farBeams[b]].GetMCFitHistogram(i)->Reset();
01706         fBeams[farBeams[b]].GetMCBackgroundHistogram(i)->Reset();
01707         fBeams[farBeams[b]].GetMCFitNuMuToNuTauHistogram(i)->Reset();
01708         fBeams[farBeams[b]].GetMCHistogram(i)->Reset();
01709         fBeams[farBeams[b]].GetMCBackgroundHistogram(i)->Reset();
01710         fBeams[farBeams[b]].GetMCNoFitNuMuToNuTauHistogram(i)->Reset();
01711      }
01712   }
01713 
01714   float scaleFactor = 1.; //what is this?  LLH
01715 
01716 //LLH uncommented
01717    //put the fit spectra into the beam objects
01718    fBeams[farBeamI].GetMCFitHistogram(NCType::kCC)->Add(fbestfitreco_rebinned, 1./scaleFactor);
01719    fBeams[farBeamI].GetMCFitHistogram(NCType::kNC)->Add(fbestfitreco_ncfit_rebinned, 1./scaleFactor);
01720    fBeams[farBeamI].GetMCBackgroundHistogram(NCType::kCC)->Add(fbestfitreconc_rebinned, 1./scaleFactor);
01721    fBeams[farBeamI].GetMCBackgroundHistogram(NCType::kNC)->Add(fbestfitrecocc_ncfit_rebinned, 1./scaleFactor);
01722    fBeams[farBeamI].GetMCFitNuMuToNuTauHistogram(NCType::kCC)->Add(fbestfitrecotau_rebinned, 1./scaleFactor);
01723    fBeams[farBeamI].GetMCFitNuMuToNuTauHistogram(NCType::kNC)->Add(fbestfitrecotau_ncfit_rebinned, 1./scaleFactor);
01724 
01725    //and the nooscillation histograms
01726    fBeams[farBeamI].GetMCHistogram(NCType::kCC)->Add(fnooscreco_rebinned, 1./scaleFactor);
01727    fBeams[farBeamI].GetMCHistogram(NCType::kNC)->Add(fnooscreco_ncfit_rebinned, 1./scaleFactor);
01728    fBeams[farBeamI].GetMCBackgroundHistogram(NCType::kCC)->Add(fnooscreconc_rebinned, 1./scaleFactor);
01729    fBeams[farBeamI].GetMCBackgroundHistogram(NCType::kNC)->Add(fnooscrecocc_ncfit_rebinned, 1./scaleFactor);
01730    fBeams[farBeamI].GetMCNoFitNuMuToNuTauHistogram(NCType::kCC)->Add(fnooscrecotau_rebinned, 1./scaleFactor);
01731    fBeams[farBeamI].GetMCNoFitNuMuToNuTauHistogram(NCType::kNC)->Add(fnooscrecotau_ncfit_rebinned, 1./scaleFactor);
01732 
01733    //put the fit spectra into the beam objects
01734    fBeams[farBeamII].GetMCFitHistogram(NCType::kCC)->Add(fbestfitreco_rebinned_post, 1./scaleFactor);
01735    fBeams[farBeamII].GetMCFitHistogram(NCType::kNC)->Add(fbestfitreco_ncfit_rebinned_post, 1./scaleFactor);
01736    fBeams[farBeamII].GetMCBackgroundHistogram(NCType::kCC)->Add(fbestfitreconc_rebinned_post, 1./scaleFactor);
01737    fBeams[farBeamII].GetMCBackgroundHistogram(NCType::kNC)->Add(fbestfitrecocc_ncfit_rebinned_post, 1./scaleFactor);
01738    fBeams[farBeamII].GetMCFitNuMuToNuTauHistogram(NCType::kCC)->Add(fbestfitrecotau_rebinned_post, 1./scaleFactor);
01739    fBeams[farBeamII].GetMCFitNuMuToNuTauHistogram(NCType::kNC)->Add(fbestfitrecotau_ncfit_rebinned_post, 1./scaleFactor);
01740 
01741    //and the nooscillation histograms
01742    fBeams[farBeamII].GetMCHistogram(NCType::kCC)->Add(fnooscreco_rebinned_post, 1./scaleFactor);
01743    fBeams[farBeamII].GetMCHistogram(NCType::kNC)->Add(fnooscreco_ncfit_rebinned_post, 1./scaleFactor);
01744    fBeams[farBeamII].GetMCBackgroundHistogram(NCType::kCC)->Add(fnooscreconc_rebinned_post, 1./scaleFactor);
01745    fBeams[farBeamII].GetMCBackgroundHistogram(NCType::kNC)->Add(fnooscrecocc_ncfit_rebinned_post, 1./scaleFactor);
01746    fBeams[farBeamII].GetMCNoFitNuMuToNuTauHistogram(NCType::kCC)->Add(fnooscrecotau_rebinned_post, 1./scaleFactor);
01747    fBeams[farBeamII].GetMCNoFitNuMuToNuTauHistogram(NCType::kNC)->Add(fnooscrecotau_ncfit_rebinned_post, 1./scaleFactor);
01748 
01749   return;
01750 }
01751 
01752 //----------------------------------------------------------------------
01753 //write out histograms etc used in the fit
01754 void NCExtrapolationNS::WriteResources()
01755 {
01756 
01757   //fill the histograms in the NCBeam objects in fBeams vector
01758   //before calling WriteResources in the base class
01759 
01760    cout <<"Writing!!!!!!" << endl;
01761 
01762   NCExtrapolation::WriteResources();
01763 
01764   // get a pointer to the current directory
01765   // this is one of the output files
01766   TDirectory* fileDir = gDirectory;
01767   TDirectory* saveDir = gDirectory;
01768   fileDir->cd("plots"+GetShortName());
01769 
01770   fTrueND->Write();
01771   fTrueNDAfterRT->Write();
01772   fTrueNDAfterRT_noeffcorr->Write(); //LLH adding
01773   fRecoND->Write();
01774   fTrueFD->Write();
01775   fTrueFD_nc->Write();
01776   fRecoFD->Write();
01777   fRecoFD_nc->Write();
01778   fTrueFDAfterBM->Write();
01779   fTrueFDAfterBM_noacccorr->Write(); //LLH addng
01780   fTrueFDAfterBM_posteffcorr->Write(); //LLH addng
01781   fTrueFDAfterOsc->Write(); //LLH addng
01782   fRecoFDAfterBM->Write();
01783 
01784   //LLH adding
01785   frecotruemat->Write();
01786   fefficiencyn->Write();
01787   fdatselND->Write();
01788   fdatselFD->Write();
01789   fdatselFD_nc->Write();
01790 
01791   //true to reco matricies
01792   for(uint i = 0; i < matarray.size(); ++i){
01793     matarray[i]->Write();
01794     matarraync[i]->Write();
01795     matarrayb[i]->Write();
01796     matarrayncb[i]->Write();
01797     //TAUS
01798     matarrayt[i]->Write();
01799     matarraybt[i]->Write();
01800 
01801     matarray_post[i]->Write();
01802     matarraync_post[i]->Write();
01803     matarrayb_post[i]->Write();
01804     matarrayncb_post[i]->Write();
01805     //TAUS
01806     matarrayt_post[i]->Write();
01807     matarraybt_post[i]->Write();
01808 
01809     for(uint j = 0; j < noosctruearray[i].size(); ++j){
01810       noosctruearray[i][j]->Write();
01811       noosctruearraypost[i][j]->Write();
01812     }
01813   }//end writing true to reco matricies
01814 
01815 
01816   fnoosctrue->Write();
01817   fnooscreco->Write();
01818   fbestfitreco->Write();
01819   fbestfitreconc->Write();
01820   fbestfitrecotau->Write();
01821 
01822   //LLH adding
01823   fnooscar->Write();
01824   fnooscar_ccselcc->Write();
01825   fnooscar_ccselnc->Write();
01826   fnooscar_ncselnc->Write();
01827   fnooscar_ncselcc->Write();
01828 
01829   fnooscar_ccselcct->Write();
01830   fnooscar_ccselnct->Write();
01831 
01832   foscreco_rebinned->Write();
01833   fnooscreco_rebinned->Write();
01834   fbestfitreco_rebinned->Write();
01835   fbestfitreconc_rebinned->Write();
01836   fbestfitrecotau_rebinned->Write();
01837   //
01838   fnoosctrue_ncfit->Write();
01839   fnooscreco_ncfit->Write();
01840   fnooscrecocc_ncfit->Write();
01841   fnooscrecotau_ncfit->Write();
01842   fbestfitreco_ncfit->Write();
01843   fbestfitrecocc_ncfit->Write();
01844   fbestfitrecotau_ncfit->Write();
01845 
01846   foscreco_ncfit_rebinned->Write();
01847 
01848   fnooscreco_ncfit_rebinned->Write();
01849   fbestfitreco_ncfit_rebinned->Write();
01850   fbestfitrecocc_ncfit_rebinned->Write();
01851   fbestfitrecotau_ncfit_rebinned->Write();
01852 
01853   // POst
01854 
01855   fnoosctrue_post->Write();
01856   fnooscreco_post->Write();
01857   fbestfitreco_post->Write();
01858   fbestfitreconc_post->Write();
01859   fbestfitrecotau_post->Write();
01860 
01861   foscreco_rebinned_post->Write();
01862 
01863   fnooscreco_rebinned_post->Write();
01864   fbestfitreco_rebinned_post->Write();
01865   fbestfitreconc_rebinned_post->Write();
01866   fbestfitrecotau_rebinned_post->Write();
01867 
01868   fnoosctrue_ncfit_post->Write();
01869   fnooscreconc_post->Write();
01870   fbestfitreco_ncfit_post->Write();
01871   fbestfitrecocc_ncfit_post->Write();
01872   fbestfitrecotau_ncfit_post->Write();
01873 
01874   foscreco_ncfit_rebinned_post->Write();
01875 
01876   fnooscreco_ncfit_rebinned_post->Write();
01877   fbestfitreco_ncfit_rebinned_post->Write();
01878   fbestfitrecocc_ncfit_rebinned_post->Write();
01879   fbestfitrecotau_ncfit_rebinned_post->Write();
01880 
01881   //
01882   fchisqall->Write();
01883   fchisqallnc->Write();
01884   fratccnc->Write();
01885 
01886   fcorf->Write();
01887   fcor->Write();
01888   fcorfall->Write();
01889   fcorall->Write();
01890   fSingleMCNearNuCC->Write();
01891   fSingleMCNearAllCC->Write();
01892   fSingleMCFarNuCC->Write();
01893   fSingleMCFarAllCC->Write();
01894   fMCNearNuCC->Write();
01895   fMCNearAllCC->Write();
01896   fMCFarNuCC->Write();
01897   fMCFarAllCC->Write();
01898   fallccmc->Write();
01899   fallccmcf->Write();
01900   fallncmcf->Write();
01901 
01902   fpurityn->Write();
01903   feffic->Write();
01904   fefficnc->Write();
01905   fefficbg->Write();
01906   fefficncbg->Write();
01907   //fratccnc->Write(); //LLH oops in here twice!
01908   frattau->Write();
01909   feffic_numcheck->Write(); //LLH adding
01910 
01911   saveDir->cd();
01912 
01913   return;
01914 }
01915 
01916 //***************************************************************************************
01917 void NCExtrapolationNS::MakePrediction(TH1F *datsel,
01918                                        TH1F* purityn,
01919                                        TH1F* efficiencyn,
01920                                        TH2F *matn_norm,
01921                                        TH1F* prediction,
01922                                        NCType::ERunType runType){
01923 
01924   MSG("NCExtrapolationNS", Msg::kInfo) << " in MakePrediction" << endl;
01925 
01926   prediction->Reset();
01927 
01928   TH2D *bmat_norm_r2 = fbmat_norm_r2;
01929   if(runType == NCType::kRunI) bmat_norm_r2 = fbmat_norm_r2_post; //LLH should this be runType == NCType::kRunII?? - fixme
01930 
01931   double ntrue[fnbins] = {0.};
01932   double ftrue[fnbins] = {0.};
01933   double eff[fnbins] = {0.};
01934 
01935   // FIDUCIAL MASS DEFINITIONS
01936   double scalef = kscalef;
01937   double mass_f = kmass_f;
01938   double mass_n = kmass_n;
01939 
01940 //this was old fv for PRL (doMC = 0 is data, doMC=1 is mc)
01941 //   if(fDoMC==-1){ // !!!!!!!!! here change
01942 //     mass_f    = 5.23*1e3+53*484*0.01*1.032; // tons
01943 //     mass_n    = 3.141592654*(1.)*(1.)*((84-17)*0.0254*7.874+(84-17)*0.01*1.032);
01944 //     scalef    = 1.00;
01945 //   }
01946 //commented out, this is the old cc fv volume - LLH
01947 
01948   double  scalemass = mass_f/mass_n;
01949   cout <<"scalemass = " << scalemass
01950        <<"\nmass_f = " << mass_f
01951        <<"\nmass_n = " << mass_n
01952        << endl;
01953 
01954   TH1F *datsel_cor = new TH1F("datsel_cor","",fnbins,fbinvect);
01955   // NEAR TRUE SPECTRA
01956   TH1F *allccmc_obt_2d = new TH1F("allccmc_obt_2d","",fnbins,fbinvect);
01957 
01958   // STEP 1 : CORRECT FOR PURITY THE "DATA SELECTED SPECTGRUM"
01959   datsel_cor->Multiply(datsel,purityn);
01960 
01961   // STEP 2 : TRANSFORM RECO = > TRUE & CORRECT FOR EFFICIENCY
01962   float k = 0.;
01963   for(int i = 0; i < fnbins; i++){
01964     ntrue[i] = 0.;
01965     eff[i] = efficiencyn->GetBinContent(i+1);
01966     k = datsel_cor->GetBinCenter(i+1);
01967 
01968     // TRANSFORMATION USING 2D MATRIX
01969     for(int j = 0; j < fnbins; j++){
01970       ntrue[i] += matn_norm->GetBinContent(j+1,i+1)*datsel_cor->GetBinContent(j+1);
01971     }
01972     if(eff[i] > 0.){
01973        fTrueNDAfterRT_noeffcorr->Fill(k, ntrue[i]); //LLH adding
01974        ntrue[i] /= eff[i];   //efficiency correction!
01975        allccmc_obt_2d->Fill(k, ntrue[i]);
01976        fTrueNDAfterRT->Fill(k, ntrue[i]);
01977     }
01978 
01979 //LLH commenting out
01980 //     MAXMSG("NCExtrapolationNS", Msg::kInfo, fnbins)
01981 //       << "near reco to true: bin center = " << k
01982 //       << ", datsel = " << datsel->GetBinContent(i+1)
01983 //       << ", datsel_cor (purity corr) = " << datsel_cor->GetBinContent(i+1)
01984 //       << ", ntrue (eff and reco->true corr) = " << ntrue[i] << endl;
01985 
01986   }//end loop over bins
01987 
01988   // STEP 3 EXTRAPOLATE & CORRECT FOR ACCEPTANCE
01989   double cor_factf = 0.;
01990 
01991   for(int i = 0; i < fnbins; i++){
01992     // Correction factor for acceptance
01993     cor_factf = fcorf->GetBinContent(i+1);
01994     ftrue[i] = 0.;
01995 
01996     // TRANSFORMATION USING 2D GNUMI MATRIX
01997     for(int j = 0; j < fnbins; j++) {
01998       ftrue[i] += bmat_norm_r2->GetBinContent(j+1,i+1)*allccmc_obt_2d->GetBinContent(j+1);
01999     } //LLH added braces for readability
02000 
02001     //LLH adding
02002     k = fcorf->GetBinCenter(i+1);
02003     fTrueFDAfterBM_noacccorr->Fill(k, ftrue[i]);
02004     //end LLH added
02005 
02006 //LLH    cout <<"cor_factf = " << cor_factf << endl;
02007 
02008     prediction->SetBinContent(i+1,ftrue[i]*cor_factf);  //applying correction for FD acceptance and event recon.
02009   }
02010 
02011   prediction->Scale(scalemass*scalef);
02012 
02013   datsel_cor->Delete();
02014   allccmc_obt_2d->Delete();
02015 }
02016 
02017 //----------------------------------------------------------------------
02018 void NCExtrapolationNS::RecoTruePur(TH2F*recotruemat1,
02019                                     TH1F *purityn1,
02020                                     TH1F *efficiencyn1,
02021                                     double enu_shift,
02022                                     int beam)
02023 {
02024    cout <<"In RecoTruePur! " << endl;
02025 
02026   recotruemat1->Reset();
02027   purityn1->Reset();
02028   efficiencyn1->Reset();
02029 
02030   TH1F selcc_true("selcc_true","",fnbins,fbinvect);
02031   TH1F selcc_reco("selcc_reco","",fnbins,fbinvect);
02032   TH1F sel_reco("sel_reco","",fnbins,fbinvect);
02033 
02034   double mc_weif=1.;
02035   double shower = 0.;
02036   double mumom = 0.;
02037   double reco_ene = 0.;
02038   int    flavour = 0;
02039   double true_enu = 0.;
02040   //double annpid = 0.;
02041   //int cc_nc = 0;
02042   //int initial_state = 0;
02043 
02044   //size of the MC arrays
02045   int sizeSignal =  0;
02046   int sizeBG = 0;
02047 
02048   //loop over bins in the beam
02049   //only use CC events in this stage
02050   NCEnergyBin *bin = 0;
02051   int numBins = fBeams[beam].GetNumberEnergyBins(NCType::kCC);
02052   for(int i = 0; i < numBins; ++i){
02053     bin = fBeams[beam].GetEnergyBin(i, NCType::kCC);
02054     sizeSignal = bin->GetMCSignalVectorSize();
02055     sizeBG = bin->GetMCBackgroundVectorSize();
02056     for(int j = 0; j < sizeSignal; ++j){
02057       bin->GetMCInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02058       shower *= enu_shift;
02059       reco_ene = shower + mumom;
02060       //only use nu_mu, no nu_mubar
02061       if(flavour == 14){
02062         selcc_true.Fill(true_enu,mc_weif);
02063         recotruemat1->Fill(reco_ene,true_enu,mc_weif);
02064         selcc_reco.Fill(reco_ene,mc_weif);
02065         fTrueND->Fill(true_enu, mc_weif); } //LLH temp to check datsel against fRecoND
02066         fRecoND->Fill(reco_ene, mc_weif);
02067 //LLH temp default brace location     }
02068         sel_reco.Fill(reco_ene,mc_weif); //LLH sel_reco includes numu and numubar
02069     }//end loop over signal events
02070 
02071     for(int j = 0; j < sizeBG; ++j){
02072       bin->GetMCBackgroundInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02073       shower *= enu_shift;
02074       reco_ene = shower + mumom;
02075       sel_reco.Fill(reco_ene,mc_weif);
02076       fRecoND->Fill(reco_ene, mc_weif); //LLH temp to check datsel against fRecoND
02077     }//end loop over background events
02078 
02079     //LLH adding
02080     int size_elec = bin->GetMCBeamNuEVectorSize();
02081     int size_tau = bin->GetMCNuTauVectorSize(); //there should be no taus in ND!
02082     for(int j = 0; j < size_tau; ++j){
02083        bin->GetMCNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02084        shower *= enu_shift;
02085        reco_ene = shower + mumom;
02086        sel_reco.Fill(reco_ene,mc_weif);
02087        fRecoND->Fill(reco_ene, mc_weif);
02088        cout <<"Found a tau in the ND MC!!!!!!!!!!" << reco_ene << endl;
02089     }
02090     for(int j = 0; j < size_elec; ++j){
02091        bin->GetMCBeamNuEInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02092        shower *= enu_shift;
02093        reco_ene = shower + mumom;
02094        sel_reco.Fill(reco_ene,mc_weif);
02095        fRecoND->Fill(reco_ene, mc_weif);
02096        cout <<"Found a CC-LIKE nue in the ND MC!!!!!!!!!!" << reco_ene << endl;
02097     }
02098     //end LLH added
02099 
02100   }//end loop over energy bins
02101 
02102 
02103   //for(uint i = 0; i < fweind.size(); ++i){
02104 
02105 //     mc_weif       = fweind[i];
02106 //     shower        = feshnd[i];
02107 //     shower        = shower*enu_shift;
02108 //     mumom         = femund[i];
02109 //     reco_ene      = shower+mumom;
02110 //     annpid        = fpidpnd[i];
02111 //     cc_nc         = fccnd[i];
02112 //     flavour       = fflavnd[i];
02113 //     initial_state = finisnd[i];
02114 //     true_enu      = fenutnd[i];
02115 
02116 //     if(cc_nc == NCType::kCC
02117 //        && flavour == 14
02118 //        && initial_state <= 2){
02119 
02120 //       selcc_true->Fill(true_enu,mc_weif);
02121 //       recotruemat1->Fill(reco_ene,true_enu,mc_weif);
02122 //       selcc_reco->Fill(reco_ene,mc_weif);
02123 //     }
02124 
02125 //     sel_reco->Fill(reco_ene,mc_weif);
02126 
02127 //   }//end loop over near events
02128 
02129   double pur = 1.;
02130   double eff = 1.;
02131   double all = 0.;
02132   double selected = 0.;
02133   for(int i = 0; i < fnbins; i++){
02134     selected = selcc_reco.GetBinContent(i+1);
02135     all = sel_reco.GetBinContent(i+1);
02136     if(all > 0){
02137       pur = selected/all;
02138       purityn1->SetBinContent(i+1,pur);
02139     }
02140 
02141     selected = selcc_true.GetBinContent(i+1);
02142     all = fallccmc->GetBinContent(i+1);
02143     if(all > 0){
02144       eff = selected/all;
02145       efficiencyn1->SetBinContent(i+1,eff);
02146     }
02147   }
02148   //Efficiency check
02149   //---------------WHAT DOES THIS DO?
02150   double check = 1.;
02151   for(int i = 0; i < fnbins; i++){
02152     check=efficiencyn1->GetBinContent(i+1);
02153     if(check == 0 && i < fnbins)
02154       efficiencyn1->SetBinContent(i+1,efficiencyn1->GetBinContent(i+2));
02155   }
02156 
02157   Convert2DMatrixToProbability(recotruemat1);
02158 
02159   return;
02160 }
02161 
02162 //----------------------------------------------------------------------
02163 //Why is this separate from GetEff??
02164 void NCExtrapolationNS::RecoTrueEffPur(TH2F *truerecomat1, //cc sel cc
02165                                        TH2F *truerecomat2, //nc sel nc
02166                                        TH2F *truerecomat1b, //nc sel cc
02167                                        TH2F *truerecomat2b, //cc sel nc
02168                                        double enu_shift,
02169                                        int beam)
02170 {
02171    cout <<"In RecoTrueEffPur! " << endl;
02172 
02173   truerecomat1->Reset();
02174   truerecomat2->Reset();
02175 
02176   truerecomat1b->Reset();
02177   truerecomat2b->Reset();
02178 
02179   double  mc_weif=1.;
02180   double  shower = 1.;
02181   double  mumom = 1.;
02182   double  reco_ene = 0.;
02183   int     flavour = 0;
02184   double  true_enu = 0.;
02185 //   double   showernc = 1.;
02186 //   int     annpid_f1p1 = 1;
02187 //   int     cc_nc = 0;
02188 //   int     initial_state = 0;
02189 
02190 //   bool trcc   = false;
02191 //   bool trcca  = false;
02192 //   bool slcc   = false;
02193 //   bool nslcc  = false;
02194 
02195 //   bool evnc =false;
02196 //   bool slnc =false;
02197 //   bool nslnc=false;
02198 
02199   //size of the MC arrays
02200   int sizeSignal =  0;
02201   int sizeBG = 0;
02202   int sizeElec =0;
02203 
02204   //loop over bins in the beam
02205   NCEnergyBin *bin = 0;
02206   int numBins = fBeams[beam].GetNumberEnergyBins(NCType::kCC);
02207   for(int i = 0; i < numBins; ++i){
02208     bin = fBeams[beam].GetEnergyBin(i, NCType::kCC);
02209     sizeSignal = bin->GetMCSignalVectorSize();
02210     sizeBG = bin->GetMCBackgroundVectorSize();
02211     sizeElec = bin->GetMCBeamNuEVectorSize();
02212 
02213     //true cc selected as cc
02214     for(int j = 0; j < sizeSignal; ++j){
02215       bin->GetMCInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02216       shower *= enu_shift;
02217       reco_ene = shower + mumom;
02218       //only use nu_mu, no nu_mubar <= original case, LLH is including numubar!
02219       if(abs(flavour) == 14){
02220         truerecomat1->Fill(true_enu, reco_ene, mc_weif);
02221         fTrueFD->Fill(true_enu, mc_weif);
02222         fRecoFD->Fill(reco_ene, mc_weif);
02223         }
02224     }//end loop over signal events
02225 
02226     //true nc selected as cc
02227     for(int j = 0; j < sizeBG; ++j){
02228       bin->GetMCBackgroundInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02229       shower *= enu_shift;
02230       reco_ene = shower + mumom;
02231       truerecomat2b->Fill(true_enu, reco_ene, mc_weif);
02232     }//end loop over background events
02233 
02234     //LLH adding
02235 //    int size_tau = bin->GetMCNuTauVectorSize(); //taus must be oscillated before adding, no osc on for debugging, therefore no taus!
02236 //     for(int j = 0; j < size_tau; ++j){
02237 //        bin->GetMCNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02238 //        shower *= enu_shift;
02239 //        reco_ene = shower + mumom;
02240 //       fRecoFD->Fill(reco_ene, mc_weif);
02241 //       cout <<"Finding FD tau!!!!" << endl;
02242 //    }
02243 
02244     //true cc elec selected as cc
02245     for(int j = 0; j < sizeElec; ++j){
02246        bin->GetMCBeamNuEInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02247        shower *= enu_shift;
02248        reco_ene = shower + mumom;
02249        truerecomat2b->Fill(true_enu, reco_ene, mc_weif);
02250 //        fRecoFD->Fill(reco_ene, mc_weif);
02251 //        cout <<"Finding FD nue!!!!" << endl;
02252     }
02253     //end LLH added
02254 
02255   }//end loop over cc energy bins
02256 
02257   numBins = fBeams[beam].GetNumberEnergyBins(NCType::kNC);
02258   for(int i = 0; i < numBins; ++i){
02259     bin = fBeams[beam].GetEnergyBin(i, NCType::kNC);
02260     sizeSignal = bin->GetMCSignalVectorSize();
02261     sizeBG = bin->GetMCBackgroundVectorSize();
02262     sizeElec = bin->GetMCBeamNuEVectorSize();
02263 
02264     //true nc selected as nc
02265     for(int j = 0; j < sizeSignal; ++j){
02266       bin->GetMCInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02267       shower *= enu_shift;
02268       reco_ene = shower + mumom;
02269       truerecomat2->Fill(true_enu,shower,mc_weif);
02270       fTrueFD_nc->Fill(true_enu, mc_weif);
02271       fRecoFD_nc->Fill(reco_ene, mc_weif);
02272     }//end loop over signal events
02273 
02274     //true cc selected as nc
02275     for(int j = 0; j < sizeBG; ++j){
02276       bin->GetMCBackgroundInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02277       shower *= enu_shift;
02278       reco_ene = shower + mumom;
02279       truerecomat1b->Fill(true_enu,reco_ene,mc_weif);
02280     }//end loop over background events
02281 
02282     //true elec selected as nc, true cc only
02283     for(int j = 0; j < sizeElec; ++j){
02284        bin->GetMCBeamNuEInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02285        shower *= enu_shift;
02286        reco_ene = shower + mumom;
02287        truerecomat1b->Fill(true_enu, reco_ene, mc_weif);
02288     }
02289 
02290   }//end loop over nc energy bins
02291 
02292 //   for(uint i = 0; i < fwei.size(); i++){
02293 //     mc_weif       = fwei[i];
02294 //     shower        = fesh[i];
02295 //     shower        = shower*enu_shift;
02296 //     showernc      = feshnc[i];
02297 //     showernc      = showernc*enu_shift;
02298 //     mumom         = femu[i];
02299 //     reco_ene      = shower+mumom;
02300 //     annpid_f1p1   = fpid_post[i];
02301 //     cc_nc         = fccfd[i];
02302 //     flavour       = fflavfd[i];
02303 //     initial_state = finis[i];
02304 //     true_enu      = fenut[i];
02305 
02306 //     // CC Selection
02307 //     trcc   = false;
02308 //     trcca  = false;
02309 //     slcc   = false;
02310 //     nslcc  = false;
02311 
02312 //     if(cc_nc == NCType::kCC
02313 //        && initial_state <= 2
02314 //        && flavour == 14)
02315 //       trcc  = true;
02316 //     if(cc_nc == NCType::kCC)
02317 //       trcca = true;
02318 //     if(annpid_f1p1 == NCType::kCC) slcc  = true;
02319 //     if(annpid_f1p1 == NCType::kNC) nslcc = true;
02320 
02321 //     // NC Selection
02322 //     evnc =false;
02323 //     slnc =false;
02324 //     nslnc=false;
02325 
02326 //     if(cc_nc == NCType::kNC) evnc = true;
02327 //     if(annpid_f1p1 == NCType::kNC) slnc = true;
02328 //     if(annpid_f1p1 == NCType::kCC) nslnc= true;
02329 
02330 //     // CC selected as CC
02331 //     if(trcc == true && slcc==true)
02332 //       truerecomat1->Fill(true_enu,reco_ene,mc_weif);
02333 
02334 //     // NC Selected as NC
02335 //     if(evnc == true && slnc==true)
02336 //       truerecomat2->Fill(true_enu,showernc,mc_weif);
02337 
02338 //     //CC selected as NC
02339 //     if(trcca==true && nslcc==true)
02340 //       truerecomat1b->Fill(true_enu,showernc,mc_weif);
02341 
02342 //     //NC selected as CC
02343 //     if(evnc==true && nslnc==true)
02344 //       truerecomat2b->Fill(true_enu,reco_ene,mc_weif);
02345 
02346 //   }
02347 
02348   Convert2DMatrixToProbability(truerecomat1);
02349   Convert2DMatrixToProbability(truerecomat1b);
02350   Convert2DMatrixToProbability(truerecomat2);
02351   Convert2DMatrixToProbability(truerecomat2b);
02352 
02353   return;
02354 }
02355 //**********************************************************************
02356 // TAUS RECO TO TRUE FOR BOTH CC AND ND NAMELY 2 HISTOS
02357 void NCExtrapolationNS::RecoTrueEffPurTau(TH2F *truerecomat1,
02358                                           TH2F *truerecomat1b,
02359                                           double enu_shift,
02360                                           int beam)
02361 {
02362 
02363   truerecomat1->Reset();
02364   truerecomat1b->Reset();
02365 
02366   double  mc_weif=1.;
02367   double   shower = 1.;
02368   double   mumom = 1.;
02369   double   reco_ene = 0.;
02370   int     flavour = 0;
02371   double   true_enu = 0.;
02372 //   double   showernc = 1.;
02373 //   int     annpid_f1p1 = 1;
02374 //   int     cc_nc = 0;
02375 //   int     initial_state = 0;
02376 //   bool    selfl = false;
02377 //   bool    selflnc = false;
02378 
02379   //size of the MC arrays
02380   int sizeTau =  0;
02381 
02382   //loop over bins in the beam
02383   //only true cc taus are included in the tau backgrounds as the nc are just nc events
02384   NCEnergyBin *bin = 0;
02385   int numBins = fBeams[beam].GetNumberEnergyBins(NCType::kCC);
02386   for(int i = 0; i < numBins; ++i){
02387     bin = fBeams[beam].GetEnergyBin(i, NCType::kCC);
02388     sizeTau = bin->GetMCNuTauVectorSize();
02389 
02390     //true cc tau selected as cc
02391     for(int j = 0; j < sizeTau; ++j){
02392       bin->GetMCNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02393       shower *= enu_shift;
02394       reco_ene = shower + mumom;
02395       truerecomat1->Fill(true_enu,reco_ene,mc_weif);
02396     }//end loop over tau events
02397   }//end loop over cc energy bins
02398 
02399   numBins = fBeams[beam].GetNumberEnergyBins(NCType::kNC);
02400   for(int i = 0; i < numBins; ++i){
02401     bin = fBeams[beam].GetEnergyBin(i, NCType::kNC);
02402     sizeTau = bin->GetMCNuTauVectorSize();
02403 
02404     //true cc tau selected as nc
02405     for(int j = 0; j < sizeTau; ++j){
02406       bin->GetMCNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02407       shower *= enu_shift;
02408       truerecomat1b->Fill(true_enu,shower,mc_weif);
02409     }//end loop over tau events
02410   }//end loop over nc energy bins
02411 
02412 //   for(uint i = 0; i < fweit.size(); i++){
02413 
02414 //     mc_weif       = fweit[i];
02415 //     shower        = fesht[i];
02416 //     shower        = shower*enu_shift;
02417 //     showernc      = feshnct[i];
02418 //     showernc      = showernc*enu_shift;
02419 //     mumom         = femut[i];
02420 //     reco_ene      = shower+mumom;
02421 //     annpid_f1p1   = fpidpt[i];
02422 //     cc_nc         = fccfdt[i];
02423 //     flavour       = fflavfdt[i];
02424 //     initial_state = finist[i];
02425 //     true_enu      = fenutt[i];
02426 
02427 //     selfl = false;
02428 //     selflnc = false;
02429 
02430 //     // CC Selection
02431 //     if(annpid_f1p1 == NCType::kCC) selfl =true;
02432 //     // NC Selection
02433 //     if(annpid_f1p1 == NCType::kNC) selflnc =true;
02434 //     // CCTAUS selected as CC
02435 //     if(cc_nc == NCType::kCC && flavour == 16){
02436 //       if(selfl) truerecomat1->Fill(true_enu,reco_ene,mc_weif);
02437 //     }
02438 //     //CCTAUS selected as NC
02439 //     if(cc_nc == NCType::kCC && flavour == 16){
02440 //       if(selflnc)  truerecomat1b->Fill(true_enu,showernc,mc_weif);
02441 //     }
02442 //   }
02443 
02444   Convert2DMatrixToProbability(truerecomat1);
02445   Convert2DMatrixToProbability(truerecomat1b);
02446 
02447   return;
02448 }
02449 
02450 //********************************************************************************
02451 void NCExtrapolationNS::GetEff(TH1F *effic_noosctrue, //cclike, cc (sig)
02452                                TH1F *effic_noosctrue_ncfit, //nclike, nc (sig)
02453                                TH1F *neffic_noosctrue, //nclike, cc (bg)
02454                                TH1F *neffic_noosctruenc, //cclike, nc (bg)
02455                                TH1F *ratccnc,
02456                                TH1F *rattau,
02457                                int beam)
02458 {
02459   ratccnc->Reset();
02460   rattau->Reset();
02461 
02462   effic_noosctrue->Reset();
02463   effic_noosctrue_ncfit->Reset();
02464 
02465   neffic_noosctrue->Reset();
02466   neffic_noosctruenc->Reset();
02467 
02468 //LLH deprecating  TH1F all("all",   "",fnbins,fbinvect);
02469   TH1F alltau("alltau","",fnbins,fbinvect);
02470 //LLH deprecating  TH1F allnc("allnc", "",fnbins,fbinvect);
02471 
02472   double mc_weif = 1.;
02473   double shower = 0.;
02474   double mumom = 0.;
02475   int    flavour = 0;
02476   double true_enu = 0.;
02477 //   double reco_ene = 0.;
02478 //   int   annpid_f1p1 = 0;
02479 //   int cc_nc = 0;
02480 //   int initial_state = 0;
02481 
02482 //   bool trcc =false;
02483 //   bool trcca=false;
02484 //   bool slcc =false;
02485 //   bool nslcc=false;
02486 //   bool evnc =false;
02487 //   bool slnc =false;
02488 //   bool nslnc=false;
02489 
02490   //size of the MC arrays
02491   int sizeTau =  0;
02492   int sizeSignal =  0;
02493   int sizeBG =  0;
02494   int sizeElec = 0;
02495 
02496   //loop over bins in the beam
02497   //only true cc taus are included in the tau backgrounds as the nc are just nc events
02498   NCEnergyBin *bin = 0;
02499   int numBins = fBeams[beam].GetNumberEnergyBins(NCType::kCC);  //Looping over CC selected
02500   for(int i = 0; i < numBins; ++i){
02501     bin = fBeams[beam].GetEnergyBin(i, NCType::kCC);
02502     sizeSignal = bin->GetMCSignalVectorSize();
02503     sizeBG = bin->GetMCBackgroundVectorSize();
02504     sizeTau = bin->GetMCNuTauVectorSize();
02505     sizeElec = bin->GetMCBeamNuEVectorSize();
02506 
02507     //true cc selected as cc
02508     for(int j = 0; j < sizeSignal; ++j){
02509       bin->GetMCInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02510       if(abs(flavour) == 14) effic_noosctrue->Fill(true_enu, mc_weif); //originally only numu, LLH including numubar!
02511     }//end loop over cc events
02512 
02513     //true nc selected as cc
02514     for(int j = 0; j < sizeBG; ++j){
02515       bin->GetMCBackgroundInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02516       neffic_noosctruenc->Fill(true_enu, mc_weif); //all 3 flavours!! (NC true only)
02517     }//end loop over nc events
02518 
02519     //true cc tau selected as cc
02520     for(int j = 0; j < sizeTau; ++j){
02521       bin->GetMCNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02522       alltau.Fill(true_enu, mc_weif); //fixme! incorrect for fv efficiency? - LLH
02523     }//end loop over tau events
02524 
02525     //true cc elec selected as cc
02526     for(int j = 0; j < sizeElec; ++j){
02527        bin->GetMCBeamNuEInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02528        neffic_noosctruenc->Fill(true_enu, mc_weif); //fixme! incorrect for fv efficiency? - LLH
02529     }//end loop over electron events
02530   }//end loop over cc energy bins
02531 
02532   numBins = fBeams[beam].GetNumberEnergyBins(NCType::kNC);  //Looping over NC selected
02533   for(int i = 0; i < numBins; ++i){
02534     bin = fBeams[beam].GetEnergyBin(i, NCType::kNC);
02535     sizeSignal = bin->GetMCSignalVectorSize();
02536     sizeBG = bin->GetMCBackgroundVectorSize();
02537     sizeTau = bin->GetMCNuTauVectorSize();
02538     sizeElec = bin->GetMCBeamNuEVectorSize();
02539 
02540     //true nc selected as nc
02541     for(int j = 0; j < sizeSignal; ++j){
02542       bin->GetMCInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02543       effic_noosctrue_ncfit->Fill(true_enu, mc_weif); //numu and numubar (and neu neubar? - LLH)
02544     }//end loop over nc events
02545 
02546     //true cc selected as nc
02547     for(int j = 0; j < sizeBG; ++j){
02548       bin->GetMCBackgroundInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02549       neffic_noosctrue->Fill(true_enu, mc_weif);  //numu and numubar
02550     }//end loop over cc events
02551 
02552     //true tau selected as nc
02553     for(int j = 0; j < sizeTau; ++j){
02554       bin->GetMCNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02555       alltau.Fill(true_enu,mc_weif);  //fixme
02556     }//end loop over tau events
02557 
02558     //true elec selected as nc, true cc only
02559     for(int j = 0; j < sizeElec; ++j){
02560        bin->GetMCBeamNuEInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02561        neffic_noosctrue->Fill(true_enu, mc_weif); //fixme! incorrect for fv efficiency? - LLH
02562     }//end loop over electron events
02563   }//end loop over nc energy bins
02564 
02565   rattau->Divide(&alltau,fallccmcf);
02566   ratccnc->Divide(fallncmcf,fallccmcf); //nc/cc - takes pred from allccmcf->allccncf, before eff corr
02567 
02568   feffic_numcheck = (TH1F*)effic_noosctrue->Clone("effic_numcheck"); //LLH adding, at this point should be identical to fTrueFD
02569 
02570   MSG("NCExtrapolationNS", Msg::kInfo) << " DONE WITH EFFICIENCIES " << endl;
02571   effic_noosctrue->Divide(fallccmcf);
02572   effic_noosctrue_ncfit->Divide(fallncmcf);
02573 
02574   neffic_noosctrue->Divide(fallccmcf);
02575   neffic_noosctruenc->Divide(fallncmcf);
02576 
02577 }
02578 
02579 //********************************************************
02580 void NCExtrapolationNS::GetEfft(TH1F *effic_noosctrue,
02581                                 TH1F *neffic_noosctrue,
02582                                 int beam)
02583 {
02584   effic_noosctrue->Reset();
02585   neffic_noosctrue->Reset();
02586 
02587   TH1F all("all",  "",fnbins,fbinvect);
02588   //TH1F *allnc   = new TH1F("allnc","",fnbins,fbinvect);
02589 
02590   double  mc_weif=1.;
02591   double  shower = 0.;
02592   double  mumom = 0.;
02593   int     flavour = 0;
02594   double  true_enu = 0.;
02595 //   double  reco_ene = 0.;
02596 //   int    annpid_f1p1 = 0;
02597 //   int    cc_nc = 0;
02598 //   int    initial_state = 0;
02599 
02600 //   bool trcc = false;
02601 //   bool slcc = false;
02602 //   bool slnc = false;
02603 
02604   //size of the MC arrays
02605   int sizeTau =  0;
02606 
02607   //loop over bins in the beam
02608   //only true cc taus are included in the tau backgrounds as the nc are just nc events
02609   NCEnergyBin *bin = 0;
02610   int numBins = fBeams[beam].GetNumberEnergyBins(NCType::kCC);
02611   for(int i = 0; i < numBins; ++i){
02612     bin = fBeams[beam].GetEnergyBin(i, NCType::kCC);
02613     sizeTau = bin->GetMCNuTauVectorSize();
02614 
02615     //true cc tau selected as cc
02616     for(int j = 0; j < sizeTau; ++j){
02617       bin->GetMCNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02618       all.Fill(true_enu,mc_weif);
02619       effic_noosctrue->Fill(true_enu,mc_weif);
02620     }//end loop over tau events
02621   }//end loop over cc energy bins
02622 
02623   numBins = fBeams[beam].GetNumberEnergyBins(NCType::kNC);
02624   for(int i = 0; i < numBins; ++i){
02625     bin = fBeams[beam].GetEnergyBin(i, NCType::kNC);
02626     sizeTau = bin->GetMCNuTauVectorSize();
02627 
02628     //true cc tau selected as cc
02629     for(int j = 0; j < sizeTau; ++j){
02630       bin->GetMCNuTauInformation(true_enu, shower, mumom, mc_weif, flavour, j);
02631       all.Fill(true_enu,mc_weif);
02632       neffic_noosctrue->Fill(true_enu,mc_weif);
02633     }//end loop over tau events
02634   }//end loop over cc energy bins
02635 
02636 //   for(uint i = 0; i < fweit.size(); i++){
02637 
02638 //     mc_weif       = fweit[i];
02639 //     mc_weif       = fweit[i];
02640 //     shower        = fesht[i];
02641 //     mumom         = femut[i];
02642 //     reco_ene      = shower+mumom;
02643 
02644 //     annpid_f1p1   = fpidpt[i];
02645 //     cc_nc         = fccfdt[i];
02646 //     flavour       = fflavfdt[i];
02647 //     initial_state = finist[i];
02648 //     true_enu      = fenutt[i];
02649 
02650 //     trcc = false;
02651 //     slcc = false;
02652 //     slnc = false;
02653 
02654 //     if(cc_nc == NCType::kCC && flavour == 16) trcc = true;
02655 //     if(annpid_f1p1 == NCType::kCC) slcc = true;
02656 //     if(annpid_f1p1 == NCType::kNC) slnc = true;
02657 
02658 //     if(trcc==true)                all->Fill(true_enu,mc_weif);
02659 //     if(slcc==true  && trcc==true) effic_noosctrue->Fill(true_enu,mc_weif);
02660 //     if(slnc==true  && trcc==true) neffic_noosctrue->Fill(true_enu,mc_weif);
02661 //   }
02662 
02663   MSG("NCExtrapolationNS", Msg::kInfo) << " DONE WITH EFFICIENCIES TAU " << endl;
02664   effic_noosctrue->Divide(&all);
02665   neffic_noosctrue->Divide(&all);
02666 
02667   return;
02668 }
02669 
02670 //***************************************************************
02671 
02672 //LLH Note: this only works if inihist has smaller bins that fit exactly inside of rebinhist bins
02673 void NCExtrapolationNS::Rebindata(TH1F *inihist,
02674                                   TH1F &rebinhist)
02675 {
02676   //loop over the output histogram and find the edges of the
02677   //bins.  loop over the input histogram and put all bins that land
02678   //within those edges into the output histogram
02679   double binCenter = 0.;
02680   double binContent = 0.;
02681   double outBinLow = 0.;
02682   double outBinHigh = 0.;
02683   double inBinLow = 0.;
02684   double inBinHigh = 0.;
02685   rebinhist.Reset();
02686   for(int i = 0; i < rebinhist.GetXaxis()->GetNbins(); ++i){
02687     outBinLow = rebinhist.GetXaxis()->GetBinLowEdge(i+1);
02688     outBinHigh = rebinhist.GetXaxis()->GetBinUpEdge(i+1);
02689     for(int j = 0; j < inihist->GetXaxis()->GetNbins(); ++j){
02690       binContent = inihist->GetBinContent(j+1);
02691       binCenter = inihist->GetBinCenter(j+1);
02692       inBinLow = inihist->GetXaxis()->GetBinLowEdge(j+1);
02693       inBinHigh = inihist->GetXaxis()->GetBinUpEdge(j+1);
02694 
02695       //if the input histogram bin is withint the output bin,
02696       //add the content to the output histogram
02697       if(binCenter >= outBinLow && binCenter <= outBinHigh){
02698         if(inBinLow >= outBinLow && inBinHigh <= outBinHigh){
02699           rebinhist.Fill(binCenter, binContent);
02700           MSG("NCExtrapolationNS", Msg::kDebug)
02701             << rebinhist.GetName() << " " << outBinLow
02702             << " " << outBinHigh << " " << inBinLow
02703             << " " << inBinHigh << " " << rebinhist.GetBinContent(i+1)
02704             << " " << binContent << endl;
02705         }
02706         else if(inBinLow >= outBinLow && inBinHigh >= outBinHigh)
02707           MSG("NCExtrapolationNS", Msg::kWarning) << "input bin goes over output bin high edge "
02708                                                   << binCenter << " "
02709                                                   << outBinLow << "/" << inBinLow << " "
02710                                                   << inBinHigh << "/" << outBinHigh << endl;
02711         else if(inBinLow <= outBinLow && inBinHigh <= outBinHigh)
02712           MSG("NCExtrapolationNS", Msg::kWarning) << "input bin goes under output bin low edge "
02713                                                   << binCenter << " "
02714                                                   << outBinLow << "/" << inBinLow << " "
02715                                                   << inBinHigh << "/" << outBinHigh << endl;
02716       }//end if input bin center is within output bin
02717 
02718     }//end loop over input histogram
02719   }//end loop over output histogram
02720 
02721   return;
02722 }
02723 
02724 //----------------------------------------------------------------------
02725 //Rebindata now works for any input histogram as it uses the information
02726 //from the rebinhist and inihist to determine how to combine the bins in
02727 //inihist
02728 void NCExtrapolationNS::Rebin(TH1F *inihist,
02729                               TH1F &rebinhist)
02730 {
02731   Rebindata(inihist, rebinhist);
02732 
02733 //   int kkk=1;
02734 
02735 //   for(int kk=1;kk<200;kk++){
02736 //     if(kk==3){
02737 //       double reco2 =inihist->GetBinContent(kk)+inihist ->GetBinContent(kk-1)+inihist ->GetBinContent(kk-2);
02738 //       rebinhist.SetBinContent(kkk,reco2);
02739 //     }
02740 //     if(kk>3 && kk<=39 && kk%2==1){
02741 //       kkk=kkk+1;
02742 //       double reco2 =inihist->GetBinContent(kk)+inihist ->GetBinContent(kk-1);
02743 //       rebinhist.SetBinContent(kkk,reco2);
02744 //     }
02745 //     if(kk>39 && kk%4==3 && kk<=119){ // 75
02746 //       kkk=kkk+1;
02747 //       double reco2 =inihist->GetBinContent(kk)+inihist ->GetBinContent(kk-1)+inihist->GetBinContent(kk-2)+inihist ->GetBinContent(kk-3);
02748 //       rebinhist.SetBinContent(kkk,reco2);
02749 //     }
02750 
02751 //   }
02752 
02753   return;
02754 }
02755 
02756 //----------------------------------------------------------------------
02757 float NCExtrapolationNS::chisq(TH1F *obs,
02758                                TH1F *exp,
02759                                TH1F *expnc)
02760 {
02761 
02762   // calculate chisq between 2 histograms - use Poisson form
02763 
02764   float chisquare=0;
02765   float nobs = 0.;
02766   float nexp = 0.;
02767   float nexpnc = 0.;
02768 
02769   for(int i = 0; i < obs->GetNbinsX(); i++){
02770     nobs = obs->GetBinContent(i+1);
02771     nexp = exp->GetBinContent(i+1);
02772 
02773     nexpnc = expnc->GetBinContent(i+1);
02774 
02775     if(nexp > 0){
02776       // LOG LIKELIHOOD
02777       if(nobs > 0)  chisquare += 2*(nexp-nobs)+2*nobs*log(nobs/nexp);
02778       if(nobs == 0) chisquare += 2*(nexp-nobs);
02779     }
02780 
02781 // Temp - LLH
02782 //     MSG("NCExtrapolationNS", Msg::kInfo) << "find chi^2 " << obs->GetBinCenter(i+1)
02783 //                                       << " " << nobs << " / " << nexp
02784 //                                       << " " << nexp/nobs << " "
02785 //                                       << " " << chisquare << endl;
02786   }
02787 
02788   return chisquare;
02789 }
02790 
02791 //----------------------------------------------------------------------
02792 void NCExtrapolationNS::truetoreco(TH1F *truth,
02793                                    TH1F& reco,
02794                                    TH2F *matrix)
02795 {
02796   //Adding this back (got lost in translation for some reason) - LLH
02797    reco.Reset();
02798 
02799   // use matrix to convert from true to reco
02800   int nbins = truth->GetNbinsX();
02801 
02802   float count = 0.;
02803   for(int i = 0; i < nbins; i++){
02804     count = 0;
02805     for(int j = 0; j < nbins; j++){
02806       count += truth->GetBinContent(j+1)*matrix->GetBinContent(j+1, i+1);
02807     }
02808     reco.SetBinContent(i+1,count);
02809 
02810 // Temp
02811 //     MAXMSG("NCExtrapolationNS", Msg::kInfo, nbins) << "far true to reco "
02812 //                                                 << ":  bin_center = " << reco.GetBinCenter(i+1) << " "
02813 //                                                 << ", reco = " << count << ", truth = " << truth->GetBinContent(i+1)
02814 //                                                 << endl;
02815   }
02816 
02817   return;
02818 }
02819 
02820 //----------------------------------------------------------------------
02821 //this method uses the value of sin^2(1.27Delta m^2 L/E) appropriate for each bin
02822 //in the provided vector.
02823 void NCExtrapolationNS::oscweight(std::vector<double> &sinSqrDeltaMSqr,
02824                                   double sinSqr2Theta, double fSterile,
02825                                   TH1F *noosc, TH1F *osc, int flavor,
02826                                   int interactionType)
02827 {
02828   double prob = 1.;
02829   double noOscVal = 0.;
02830 
02831   for(uint i = 0; i < sinSqrDeltaMSqr.size(); ++i){
02832     noOscVal = noosc->GetBinContent(i+1);
02833 
02834     if(TMath::Abs(flavor) == 14){
02835       if(interactionType == NCType::kCC)
02836         prob = 1. - sinSqr2Theta*sinSqrDeltaMSqr[i];
02837       else if(interactionType == NCType::kNC)
02838         prob = 1. - sinSqr2Theta*fSterile*sinSqrDeltaMSqr[i];
02839     }
02840     else if(TMath::Abs(flavor) == 16){
02841       if(interactionType == NCType::kCC)
02842         prob = (1. - fSterile)*sinSqr2Theta*sinSqrDeltaMSqr[i];
02843       else if(interactionType == NCType::kNC)
02844         prob = 0.;
02845     }
02846 
02847     osc->SetBinContent(i+1, prob*noOscVal);
02848   }//end loop over energy bins
02849 
02850   return;
02851 }
02852 
02853 //----------------------------------------------------------------------
02854 void NCExtrapolationNS::oscweightOld(TH1F *noosc,
02855                                      TH1F &osc,
02856                                      float dm,
02857                                      float st,
02858                                      int method,
02859                                      int what)
02860 {
02861 
02862   // what = 1 cc disappearance
02863   // what = 2 nc disappearance
02864   // what = 11 taus appearance
02865   // what = -1 decay
02866 
02867   // reweight spectrum with oscillation parameters
02868   // interpolation method:  0 - none,  1 - quadratic,  anything else - linear
02869   //osc.Reset();
02870   float distance = NCType::kBaseLineFar;
02871 
02872   int nbins = noosc->GetNbinsX();
02873   float prob = 0.;
02874   float binsize = 0.25;
02875   float wgt1 = 0.;
02876   float wgt2 = 0.;
02877   float wgt3 = 0.;
02878   float de = 0.;
02879   float de2 = 0.;
02880 
02881   float avgosc = 0.;
02882   float avgosc2 = 0.;
02883   float avgwgt1 = 0.;
02884   float avgwgt2 = 0.;
02885   float enu = 0.;
02886   float delta = 0.;
02887   float linwgt = 0.;
02888   float quadwgt = 0.;
02889 
02890   for (int i = 0; i < nbins; i++) {
02891 
02892     if(what == 1 || what == 0)
02893       prob = 1-st*TMath::Power(TMath::Sin(NCType::k127*dm*distance/noosc->GetBinCenter(i+1)),2);
02894     if(what == 11)
02895       prob = st*TMath::Power(TMath::Sin(NCType::k127*dm*distance/noosc->GetBinCenter(i+1)),2);
02896     if(what == -1)
02897       prob=TMath::Power(st+(1-st)*TMath::Exp(-dm*distance/(2*noosc->GetBinCenter(i+1))),2);
02898     if(what == -11)
02899       prob = 0.;
02900 
02901     binsize = 0.25;
02902     if (i == 0)    binsize = 0.5;
02903     if (i == 199)  binsize = 200.-49.75;
02904 
02905     if(i > nbins-2 || method == 0){
02906       osc.SetBinContent(i+1, noosc->GetBinContent(i+1)*prob);
02907     }
02908     else{
02909 
02910       wgt1 = noosc->GetBinContent(i+1);
02911       wgt2 = noosc->GetBinContent(i+2);
02912       wgt3 = noosc->GetBinContent(i+3);
02913 
02914       if(i == 0){
02915         wgt2 = noosc->GetBinContent(i+2)+noosc->GetBinContent(i+3);
02916         wgt3 = noosc->GetBinContent(i+4)+noosc->GetBinContent(i+5);
02917 
02918       }
02919 
02920       de = wgt2-wgt1;
02921       de2 = wgt3-2*wgt2+wgt1;
02922 
02923       avgosc = 0.;
02924       avgosc2 = 0.;
02925       avgwgt1 = 0.;
02926       avgwgt2 = 0.;
02927 
02928       for(int k = 0; k < 10; k++){
02929 
02930         enu = noosc->GetBinLowEdge(i+1)+k*binsize/10.+binsize/20.;
02931 
02932         delta = enu-noosc->GetBinCenter(i+1);
02933 
02934         linwgt = wgt1+de*delta;
02935         quadwgt = wgt1+de*delta+de2/2*TMath::Power(delta,2);
02936 
02937         if(linwgt < 0) linwgt = 0.;
02938         if(quadwgt < 0) quadwgt = 0.;
02939 
02940         if(what == 1 || what == 0){
02941           avgosc += quadwgt*(1-st*TMath::Power(TMath::Sin(NCType::k127*dm*distance/enu),2));
02942           avgosc2 += linwgt*(1-st*TMath::Power(TMath::Sin(NCType::k127*dm*distance/enu),2));
02943         }
02944         if(what==11){
02945           avgosc += quadwgt*(st*TMath::Power(TMath::Sin(NCType::k127*dm*distance/enu),2));
02946           avgosc2 += linwgt*(st*TMath::Power(TMath::Sin(NCType::k127*dm*distance/enu),2));
02947         }
02948         if(what==-1){
02949           avgosc += quadwgt*TMath::Power(st+(1-st)*exp(-dm*distance/(2*enu)),2);
02950           avgosc2 += linwgt*TMath::Power(st+(1-st)*exp(-dm*distance/(2*enu)),2);
02951         }
02952         if(what==-11){
02953           avgosc += quadwgt*0.;
02954           avgosc2 += linwgt*0.;
02955         }
02956         avgwgt1 += quadwgt;
02957         avgwgt2 += linwgt;
02958       }
02959 
02960       if(method == 1 && avgosc > 0){
02961         if(avgwgt1 > 0) avgosc /= avgwgt1;
02962         else{
02963           MSG("NCExtrapolationNS", Msg::kWarning) << " attempted divide by 0 to get average "
02964                                                   << " oscillation for bin -  not averaged "
02965                                                   << avgosc << " " << avgwgt1 << " "
02966                                                   << osc.GetBinCenter(i+1) << " "
02967                                                   << noosc->GetBinContent(i+1) << endl;
02968         }
02969         osc.SetBinContent(i+1, avgosc*noosc->GetBinContent(i+1));
02970       }
02971       else if(avgwgt2 > 0){
02972         if(avgwgt2 > 0) avgosc2 /= avgwgt2;
02973         else{
02974           MSG("NCExtrapolationNS", Msg::kWarning) << " attempted divide by 0 to get average "
02975                                                   << " oscillation for bin -  not averaged "
02976                                                   << avgosc2 << " " << avgwgt2 << endl;
02977         }
02978         osc.SetBinContent(i+1, avgosc2*noosc->GetBinContent(i+1));
02979       }
02980 
02981     }
02982   }//end loop over energy bins
02983 
02984   if((st == 0 || dm == 0) && (what == 1 || what == 0)){
02985     for(int i = 0; i < nbins; i++) osc.SetBinContent(i+1,noosc->GetBinContent(i+1));
02986   }
02987 
02988   if((st == 0 || dm == 0) && (what == 11 || what == -11)){
02989     for(int i = 0; i < nbins; i++) osc.SetBinContent(i+1,0.);
02990   }
02991 
02992   if(dm == 0 && what == -1){
02993     for(int i = 0; i < nbins; i++) osc.SetBinContent(i+1,noosc->GetBinContent(i+1));
02994   }
02995 
02996   return;
02997 }
02998 
02999 //----------------------------------------------------------------------
03000 //matricies are either reco to true or true to reco.  if reco to true,
03001 //the x axis is reconstructed energy, the y axis is true energy.
03002 //if true to reco, the x axis is true energy, the y axis is reco'ed
03003 //energy
03004 void NCExtrapolationNS::Convert2DMatrixToProbability(TH2F *matrix)
03005 {
03006   double sum[fnbins];
03007   double maxf[fnbins][fnbins];
03008 
03009   //initialize
03010   for(int j = 0; j < fnbins; j++){
03011     sum[j] = 0.;
03012     for(int i = 0; i < fnbins; i++) maxf[i][j] = 0.;
03013   }
03014 
03015   //Make True to reco or reco to true
03016   //this step takes the y bins and puts them in the
03017   //x bins and vice versa
03018   for(int i = 0; i < fnbins; i++){
03019     for(int j = 0; j < fnbins; j++){
03020       maxf[j][i] = matrix->GetBinContent(i+1, j+1);
03021     }
03022   }
03023 
03024   // MAKE 2D MATRIX A PROBABILITY ONE
03025 //LLH in here twice  for(int i = 0; i < fnbins; i++) sum[i] = 0.;
03026 
03027   //sum over the new x bins for each y value
03028   for(int j = 0; j < fnbins; j++){
03029     for(int i = 0; i < fnbins; i++){
03030       sum[j] += maxf[i][j];
03031     }
03032   }
03033 
03034   //each x,y bin is the fraction of total for
03035   //given row in y.  set bin content in matrix to
03036   //be true,reco (reco,true)
03037   for(int i = 0; i < fnbins; i++){
03038     for(int j = 0; j < fnbins; j++){
03039       if(sum[j] > 0.){
03040         maxf[i][j] /= sum[j];
03041         matrix->SetBinContent(j+1, i+1, maxf[i][j]);
03042       }
03043     }
03044   }
03045 
03046   return;
03047 }
03048 
03049 //----------------------------------------------------------------------
03050 void NCExtrapolationNS::FillDataHistogramFromEnergyBins(int beam,
03051                                                         int nccc,
03052                                                         TH1F *hist)
03053 {
03054   hist->Reset();
03055 
03056   NCEnergyBin *bin = 0;
03057   int sizeData = 0;
03058   int numBins = fBeams[beam].GetNumberEnergyBins(nccc);
03059   double energy = 0.;
03060   double weight = 0.;
03061   for(int i = 0; i < numBins; ++i){
03062     bin = fBeams[beam].GetEnergyBin(i, nccc);
03063     sizeData = bin->GetDataVectorSize();
03064     for(int j = 0; j < sizeData; ++j){
03065       bin->GetDataInformation(energy, weight, j);
03066       hist->Fill(energy, weight);
03067     }
03068   }//end loop over bins
03069 
03070   return;
03071 }