loon: NueExtrapolationJB.cxx Source File

00001 
00002 // THis will be the actual engine that handles a Full Extrapolation
00003 #include <vector>
00004 #include "NueData.h"
00005 #include "NueAna/NueAnaTools/NueConvention.h"
00006 #include "NueAna/Extrapolation/NueExtrapolationJB.h"
00007 #include "NueAna/NueStandard.h"
00008 
00009 #include "TGraphAsymmErrors.h"
00010 
00011 using namespace std;
00012 
00013 NueExtrapolationJB::NueExtrapolationJB():
00014  fNZBins(0), fZBins(0)
00015 {
00016    fCurrentSys = 0;
00017    outFileName = "DefaultOut.root";
00018    fTargetPOT = 3.25e20;
00019    fUseMCAsData = true;
00020    fUseMCAsCCData = true;
00021    fNearCCPOT = 0.0;
00022    fFarCCPOT = 0.0;
00023 
00024    fNearDataPOT = 0.0;
00025    fNearCCDataPOT = 0.0;
00026 
00027    fExtrapMethod = 1;
00028 }
00029 
00030 void NueExtrapolationJB::SetNueRecoBins(int nx,Double_t lx,Double_t ux)
00031 {
00032   fNXBins = nx;
00033   fXBins = new Double_t[fNXBins+1];
00034   Float_t bwidth = (ux-lx)/float(fNXBins);
00035   for(int i=0;i<fNXBins+1;i++) fXBins[i] = lx + float(i)*bwidth;
00036 }
00037                                                                                 
00038 void NueExtrapolationJB::SetTrueBins(int ny,Double_t ly,Double_t uy)
00039 {
00040   fNYBins = ny;
00041   fYBins = new Double_t[fNYBins+1];
00042   Float_t bwidth = (uy-ly)/float(fNYBins);
00043   for(int i=0;i<fNYBins+1;i++) fYBins[i] = ly + float(i)*bwidth;
00044 }
00045                                                                                 
00046 void NueExtrapolationJB::SetCCRecoBins(int nz,Double_t lz,Double_t uz)
00047 {
00048   fNZBins = nz;
00049   fZBins = new Double_t[fNZBins+1];
00050   Float_t bwidth = (uz-lz)/float(fNZBins);
00051   for(int i=0;i<fNZBins+1;i++) fZBins[i] = lz + float(i)*bwidth;
00052 }
00053 
00054 void NueExtrapolationJB::Initialize()
00055 {
00056    this->Init();
00057 }
00058 
00059 bool NueExtrapolationJB::LoadFiles()
00060 {
00061    std::map<TChain*, int>::iterator mapBeg =  fChainDataMap.begin();
00062    std::map<TChain*, int>::iterator mapEnd =  fChainDataMap.end();
00063 
00064    int chainCount = 0;
00065    int events = 0;
00066 
00067    while(mapBeg != mapEnd)
00068    {
00069       if(!mapBeg->first) cout<<"Error invalid chain"<<endl;
00070       TChain *chain = mapBeg->first;
00071       if(strstr(chain->GetName(), "ana_nue")!=0)  this->SetUpNueAnaChain(chain);
00072       if(strstr(chain->GetName(), "NueMini")!=0)  this->SetUpNueMiniChain(chain);
00073 
00074       int nEntries = chain->GetEntries();
00075       for(int i=0;i<nEntries;i++){
00076         if(i%10000==0) std::cout << "Processed "
00077                            << 100*Float_t(i)/Float_t(nEntries)
00078                            << "% of Set " << chainCount << std::endl;
00079          chain->GetEntry(i);
00080 //         if(i > 10000) break;// && 
00081 //            (fData[mapBeg->second]->GetDetector() != 2
00082 //            || !fData[mapBeg->second]->IsNueData())) break; 
00083  
00084          if(strstr(chain->GetName(), "ana_nue")!=0){
00085             NueConvention::NueEnergyCorrection(fRecord);
00086             fData[mapBeg->second]->AddEvent(fRecord);
00087          }
00088          if(strstr(chain->GetName(),"NueMini")!=0)
00089             fData[mapBeg->second]->AddEvent(fMini);
00090 
00091       }
00092       mapBeg++;
00093       events += chain->GetEntries();
00094       chainCount++;
00095    }
00096 
00097    for(unsigned int i = 0; i <  fData.size(); i++) {
00098       if(fData[i]->GetDetector() == Detector::kFar) continue;
00099       if(fData[i]->IsData() || fUseMCAsData){
00100          if(fData[i]->IsNueData()) fNearDataPOT   = fData[i]->GetPOT();
00101          else                      fNearCCDataPOT = fData[i]->GetPOT();
00102       }
00103       if(!fUseMCAsData){ 
00104          if(fData[i]->IsNueData()) fNearDataPOT   = 1e19;
00105          else                      fNearCCDataPOT = 1e19;
00106       } 
00107    }
00108 
00109 //   MSG("NueExtrapolationJB", Msg::kInfo) << "Loaded Files "
00110 //      << chainCount<<" Chains yieled "<<events<<" events."<<endl;
00111    std::cout<<"Finished Loading all files"<<std::endl;
00112 
00113    TFile Input("NueAna/data/xsec_minos_modbyrs4_v3_5_0_mk.root", "READ");
00114 
00115    string names[5] = {"tot", "qe", "res", "dis", "coh"};
00116 
00117    for(int i = 0; i < 5; i++){
00118      string id = "h_numu_cc_" + names[i];
00119      Input.GetObject(id.c_str(), fNuMuCCXSec[i]);
00120      id = "h_nutau_cc_" + names[i];
00121      Input.GetObject(id.c_str(), fNuTauCCXSec[i]);
00122      id = "h_nue_cc_" + names[i];
00123      Input.GetObject(id.c_str(),fNueCCXSec[i]);
00124 
00125      id = "h_numubar_cc_" + names[i];
00126      Input.GetObject(id.c_str(), fNuMuBarCCXSec[i]);
00127      id = "h_nutaubar_cc_" + names[i];
00128      Input.GetObject(id.c_str(), fNuTauBarCCXSec[i]);
00129      id = "h_nuebar_cc_" + names[i];
00130      Input.GetObject(id.c_str(),fNueBarCCXSec[i]);
00131 
00132      fNueBarCCXSec[i]->SetDirectory(0);
00133      fNueCCXSec[i]->SetDirectory(0);
00134      fNuMuBarCCXSec[i]->SetDirectory(0);
00135      fNuMuCCXSec[i]->SetDirectory(0);
00136      fNuTauBarCCXSec[i]->SetDirectory(0);
00137      fNuTauCCXSec[i]->SetDirectory(0);
00138    }
00139 
00140    for(int k = 0; k < fNuMuCCXSec[0]->GetNbinsX()+1; k++){
00141       fXSecPos[k] = fNuMuCCXSec[0]->GetBinLowEdge(k);
00142    }
00143 
00144    return true;
00145 }
00146 
00147 void NueExtrapolationJB::AddChain(TChain *chain, double POT, bool isNue)
00148 {
00149    if(strstr(chain->GetName(), "ana_nue")!=0)  this->SetUpNueAnaChain(chain);
00150    if(strstr(chain->GetName(), "NueMini")!=0)  this->SetUpNueMiniChain(chain);
00151 
00152    chain->GetEntry(0);
00153 
00154    ReleaseType::Release_t rel = ReleaseType::kUnknown;
00155    Detector::Detector_t det = Detector::kUnknown;
00156    BeamType::BeamType_t beam = BeamType::kUnknown;
00157  
00158    if(strstr(chain->GetName(),"ana_nue") != 0){
00159       rel = fRecord->GetHeader().GetRelease();
00160       det = fRecord->GetHeader().GetVldContext().GetDetector();
00161       beam = fRecord->GetHeader().GetBeamType();
00162    }
00163    if(strstr(chain->GetName(),"NueMini") !=0){
00164       rel = fMini->fRelease;
00165       det = fMini->fDet;
00166       beam = fMini->fBeam;
00167    }
00168 
00169    int pos = -1;
00170 
00171    for(unsigned int i = 0; i <  fData.size(); i++) {
00172       if(fData[i]->GetBeamType() == beam && fData[i]->GetRelease() == rel
00173            && fData[i]->GetDetector() == det && fData[i]->IsNueData() == isNue)
00174       {     pos = i; break;  }
00175    }
00176    if(pos < 0){
00177      NueData* dataSet = new NueData(beam, det, rel, isNue);
00178      fData.push_back(dataSet);
00179      pos = fData.size() - 1;
00180    }
00181 
00182    if(det == Detector::kNear && isNue) fNearPOT += POT;
00183    if(det == Detector::kFar && isNue) fFarPOT += POT;
00184    if(det == Detector::kNear && !isNue) fNearCCPOT += POT;
00185    if(det == Detector::kFar && !isNue) fFarCCPOT += POT;   
00186    
00187    fData[pos]->AddPOT(POT);
00188 
00189    fChainDataMap[chain] = pos;
00190 }
00191 
00192 void NueExtrapolationJB::PrepareExtrapHistograms(Selection::Selection_t sel)
00193 {
00194    ResetExtrapHistograms();
00195 
00196    Double_t totWeight = 1.0;
00197    Selection::Selection_t tSel = sel;
00198 
00199    for(unsigned int i = 0; i < fData.size(); i++){
00200      if(fData[i]->IsData() ) continue;   //Don't need to remake Data Hists
00201      Detector::Detector_t det = fData[i]->GetDetector();
00202 
00203      NueHeader nh;
00204      fData[i]->SetupNueHeader(nh); 
00205      NueRecord* nr = new NueRecord(nh);
00206 
00207      bool isNue = fData[i]->IsNueData();
00208      tSel = sel;
00209      if(!isNue)  tSel = Selection::kCC;     
00210 
00211      for(int j = 0; j < fData[i]->NumberOfEntries(); j++)
00212      {
00213        fData[i]->FillRecord(nr, j);
00214 
00215        Background::Background_t bg =
00216           Background::TranslateFromNueClass(nr->mctrue.fNueClass);
00217 
00218        Double_t OscSysVal = fCurrentSys->GetSysValue(Systematic::kOscProb);
00219 
00220        if(fExtrapMethod != 1){         
00221          //If we are going for speed we still oscillate the NC since it doesn't matter
00222          // we need the full BNue and NuMu without oscillations to get true vs reco spectrum
00223          // for the nue and the nutau the oscillations only show up in the selected CC spectrum
00224          //  or in the true to reco which we should oscillate since it changes the map slightly
00225          if(bg == Background::kNuMuCC || bg == Background::kBNueCC
00226             || tSel == Selection::kCC)
00227          fCurrentSys->SetSysValue(Systematic::kOscProb, 0);
00228        }                                                                                      
00229        totWeight = fCurrentSys->UpdateRecord(nr, tSel);
00230 
00231        if(fExtrapMethod != 1){
00232          // reset osc state
00233          if(bg == Background::kNuMuCC || bg == Background::kBNueCC
00234             || tSel == Selection::kCC)
00235          fCurrentSys->SetSysValue(Systematic::kOscProb, OscSysVal);
00236        }
00237        
00238        double trueEnergy = nr->mctrue.nuEnergy;
00239        double recoEnergy = this->GetNueEnergy(nr, tSel);
00240 
00241        if(!isNue){
00242          bg = Background::kSelCC;
00243 
00244          if(nr->mctrue.fNueClass == 1 && det == Detector::kFar){
00245            fHistMap[bg]->fFD_numu_TrueEnergyBase->Fill(trueEnergy,totWeight);
00246          }
00247        }
00248       
00249        bool PassCuts = NueStandard::PassesPreSelection(nr) && NueStandard::PassesPIDSelection(nr, tSel);
00250        if(!isNue) PassCuts = NueStandard::PassesCCSelection(nr);
00251 
00252        if(det == Detector::kFar){
00253          //nutau and nue osc need the full fiducial voume true to reco spectrum
00254          //the others want true to reco after they have been selected
00255          if(bg == Background::kNuTauCC || bg == Background::kNueCC ){
00256            fHistMap[bg]->fFD_RecoVsTrue->Fill(recoEnergy,trueEnergy,totWeight);
00257          }else{
00258            // If we are using the accurate method we don't need to do anything special
00259            // to make things fast though we separately need to track the numu's that survive
00260            // and the numu's that appear from bnue -> numu
00261 
00262            if(PassCuts && fExtrapMethod == 1)  
00263                 fHistMap[bg]->fFD_RecoVsTrue->Fill(recoEnergy,trueEnergy,totWeight);
00264            if(PassCuts && fExtrapMethod != 1)
00265            {
00266              if(bg == Background::kNuMuCC)
00267              {
00268                bool wasNumu = TMath::Abs(nr->mctrue.nonOscNuFlavor) == 14;
00269                if(wasNumu)  fHistMap[bg]->fFD_RecoVsTrue->Fill(recoEnergy,trueEnergy,totWeight);
00270                else fHistMap[bg]->fFD_RecoVsTrue_BNue->Fill(recoEnergy,trueEnergy,totWeight);
00271              }else
00272                fHistMap[bg]->fFD_RecoVsTrue->Fill(recoEnergy,trueEnergy,totWeight);
00273            }
00274          }
00275        }
00276        
00277        if(PassCuts)
00278        {
00279           if(det == Detector::kNear){ 
00280              fHistMap[bg]->fND_RecoEnergy->Fill(recoEnergy,totWeight);
00281              fHistMap[bg]->fND_TrueEnergy->Fill(trueEnergy,totWeight);
00282           }
00283           if(det == Detector::kFar){
00284              fHistMap[bg]->fFD_RecoEnergy->Fill(recoEnergy,totWeight);
00285              fHistMap[bg]->fFD_TrueEnergy->Fill(trueEnergy,totWeight);
00286           }
00287        }
00288      }
00289      delete nr;
00290    }
00291 
00292    return;
00293 }
00294 
00295 void NueExtrapolationJB::MakeDataHistograms(Selection::Selection_t sel)
00296 {
00297    // So this isn't correct at least for Near Det
00298    // There I will be given Histograms for Data
00299      Double_t totWeight = 1.0;
00300      Selection::Selection_t tSel = sel;
00301 
00302      NueSystematic *nueSys = new NueSystematic("Standard");
00303      nueSys->AddSystematic(Systematic::kOscProb,Systematic::GetDefaultValue(Systematic::kOscProb));
00304      nueSys->SetOscParams(0.554,0.7854,0.,8.2e-5,2.4e-3,0,1);
00305      nueSys->AddSystematic(Systematic::kSKZP,0);
00306 
00307      for(unsigned int i = 0; i < fData.size(); i++){
00308        if(fData[i]->IsData()) continue;   //Don't need to deal with Data
00309   
00310        Detector::Detector_t det = fData[i]->GetDetector();
00311 //       if(det != Detector::kNear) continue;
00312 
00313        bool isNue = fData[i]->IsNueData();
00314        tSel = sel;
00315        if(!isNue)  tSel = Selection::kCC;
00316 
00317        NueHeader nh;
00318        fData[i]->SetupNueHeader(nh);
00319        NueRecord* nr = new NueRecord(nh);
00320 
00321        for(int j = 0; j < fData[i]->NumberOfEntries(); j++)
00322        {
00323           fData[i]->FillRecord(nr, j);                                                                                                      
00324           totWeight = nueSys->UpdateRecord(nr, tSel);
00325           double recoEnergy = this->GetNueEnergy(nr, tSel);           
00326 
00327           Background::Background_t bg =
00328              Background::TranslateFromNueClass(nr->mctrue.fNueClass);
00329           if(!isNue)  bg = Background::kSelCC;
00330 
00331           bool PassCuts = NueStandard::PassesPreSelection(nr) && NueStandard::PassesPIDSelection(nr, tSel);
00332           if(!isNue) PassCuts = NueStandard::PassesCCSelection(nr);
00333                                                                                                          
00334           if(PassCuts){
00335            if(det == Detector::kNear)
00336               fDataMap[bg]->fND_RecoEnergy->Fill(recoEnergy,totWeight);  
00337            if(det == Detector::kFar)
00338               fDataMap[bg]->fFD_RecoEnergy->Fill(recoEnergy,totWeight);
00339           }
00340        }
00341        delete nr;
00342      }
00343    
00344 }
00345 
00346 void  NueExtrapolationJB::LoadDataHistograms(string /*fname*/, Selection::Selection_t sel)
00347 {
00348   if(fUseMCAsData)  MakeDataHistograms(sel);
00349   else{
00350     TFile *fCC = new TFile("DataHists/NearDetCCData.root", "READ");
00351     TH1D* numuCC = 0;
00352 
00353     fCC->GetObject("ccdata", numuCC);
00354     numuCC->SetDirectory(0);
00355 
00356     
00357 
00358     fDataMap[Background::kSelCC]->fND_RecoEnergy = new TH1D("ND_Data_RecoEnergy_SelCC", "ND_Data_RecoEnergy_SelCC", fNZBins, fZBins);
00359     fDataMap[Background::kSelCC]->fND_RecoEnergy->SetDirectory(0);
00360  
00361     int k = 0;
00362     for(int i = 0; i < fNZBins; i++)
00363     {
00364        double pos = fDataMap[Background::kSelCC]->fND_RecoEnergy->GetBinCenter(i);
00365        while(pos > numuCC->GetBinCenter(k)) k++;
00366        
00367        if(pos == numuCC->GetBinCenter(k)){
00368          fDataMap[Background::kSelCC]->fND_RecoEnergy->SetBinContent(i, numuCC->GetBinContent(k));
00369          fDataMap[Background::kSelCC]->fND_RecoEnergy->SetBinError(i, numuCC->GetBinError(k));
00370        }       
00371     }
00372 
00373     string method = fSeparation;
00374     TFile *f1 = 0;
00375     TH1D* nc = 0;
00376     TH1D* cc = 0;
00377     TH1D* bnue = 0;
00378 
00379     nc = new TH1D("ND_RecoEnergy","ND Reco Energy",fNXBins, fXBins);
00380     nc->Sumw2();
00381     cc = new TH1D("FD_RecoEnergy","FD Reco Energy",fNXBins, fXBins);
00382     cc->Sumw2();
00383     bnue = new TH1D("FD_RecoEnergy_BNue","FD Reco Energy",fNXBins, fXBins);
00384     bnue->Sumw2();
00385                                                                                                   
00386     nc->SetDirectory(0);
00387     cc->SetDirectory(0);
00388     bnue->SetDirectory(0);
00389 
00390     f1 = new TFile(inDataFileName.c_str());
00391 
00392     if(f1 == 0){
00393        cout<<"Error loading files"<<endl;
00394        return;
00395     }
00396 
00397     TGraphAsymmErrors* INcc;
00398     TGraphAsymmErrors* INnc;
00399     TH1D* INbnue;
00400                                                                                 
00401     string idstring = "TG_NC_" + string(Selection::AsString(sel));
00402     f1->GetObject(idstring.c_str(), INnc);
00403     idstring = "TG_CC_" + string(Selection::AsString(sel));
00404     f1->GetObject(idstring.c_str(), INcc);
00405     idstring = "MC_BNue_" + string(Selection::AsString(sel));
00406     f1->GetObject(idstring.c_str(), INbnue);
00407     double x,y;
00408     if(INcc == 0 || INnc == 0 || INbnue == 0){
00409          cout<<"Error loading histograms"<<endl;
00410     }
00411 
00412     for(int i = 0; i < INnc->GetN(); i++)
00413     {
00414        INnc->GetPoint(i,x,y);
00415        nc->Fill(x,y);
00416        nc->SetBinError(nc->FindBin(x), INnc->GetErrorY(i));
00417        INcc->GetPoint(i,x,y);
00418        cc->Fill(x,y);
00419        cc->SetBinError(cc->FindBin(x), INcc->GetErrorY(i));
00420     }
00421                                                                                 
00422     for(int i = 1; i < 20; i++){
00423        double pos = INbnue->GetBinCenter(i);
00424        double val = INbnue->GetBinContent(i);
00425                                                                                 
00426        bnue->Fill(pos, val);
00427        bnue->SetBinError(bnue->FindBin(pos), INbnue->GetBinError(i));
00428     }
00429                                                                                 
00430     if(nc == 0 || cc == 0 || bnue == 0) cout<<"Failure to load files "<<nc<<"  "<<cc<<"  "<<bnue<<endl;
00431     nc->SetDirectory(0);
00432     cc->SetDirectory(0);
00433     bnue->SetDirectory(0);
00434 
00435     fDataMap[Background::kNuMuCC]->fND_RecoEnergy = (TH1D*) cc->Clone("ND_Data_RecoEnergy_NuMuCC");
00436     fDataMap[Background::kNuMuCC]->fND_RecoEnergy->SetDirectory(0);
00437 
00438     fDataMap[Background::kNC]->fND_RecoEnergy = (TH1D*) nc->Clone("ND_Data_RecoEnergy_NC");
00439     fDataMap[Background::kNC]->fND_RecoEnergy->SetDirectory(0);
00440 
00441     fDataMap[Background::kBNueCC]->fND_RecoEnergy = (TH1D*) bnue->Clone("ND_Data_RecoEnergy_BNueCC");
00442     fDataMap[Background::kBNueCC]->fND_RecoEnergy->SetDirectory(0);
00443 
00444      delete nc;
00445      delete cc;
00446      delete bnue;
00447      delete f1;
00448      delete fCC;
00449      delete numuCC;
00450   }
00451 
00452    TFile* EffInput = new TFile("NueAna/data/tau.root", "READ");;
00453 
00454    string name="tau_eff_"+ string(Selection::AsString(sel));  
00455 
00456    EffInput->GetObject(name.c_str(), fTauEff);
00457    if(fTauEff == 0) cout<<"error loading efficiency "<<name<<endl; 
00458    fTauEff->SetDirectory(0);
00459 
00460    delete EffInput;
00461 
00462    EffInput = new TFile("NueAna/data/MREEff.root", "READ");;
00463                                                                                 
00464    name="eff_"+ string(Selection::AsString(sel));
00465                                                                                 
00466    EffInput->GetObject(name.c_str(), fMREEff);
00467    if(fMREEff == 0) cout<<"error loading efficiency "<<name<<endl;
00468    fMREEff->SetDirectory(0);
00469                                                                                 
00470    delete EffInput;
00471 
00472    cout<<"Finished Loading the Data"<<endl;
00473 }
00474 
00475 void NueExtrapolationJB::MakePrediction()
00476 {
00477   std::map<Background::Background_t, TH1D*>::iterator mapBeg =  fPredMap.begin();
00478   std::map<Background::Background_t, TH1D*>::iterator mapEnd =  fPredMap.end();
00479                                                                                                                                  
00480   while(mapBeg != mapEnd){
00481      if(mapBeg->second){
00482             delete (mapBeg->second);
00483             mapBeg->second = NULL;
00484      }
00485      mapBeg->second = (TH1D*) GetPrediction(mapBeg->first);
00486      mapBeg++;
00487   }    
00488 }
00489 
00490 
00491 TH1* NueExtrapolationJB::GetPrediction(Background::Background_t bg)
00492 {
00493   TH1D* extrapHist = NULL;
00494   TString name = "FarPrediction_" + string(Background::AsString(bg));
00495                                                                                 
00496   //NuMu and NC get a straight F/N extrapolation
00497   if(bg == Background::kNuMuCC || bg == Background::kNC){
00498     extrapHist = (TH1D*) fDataMap[bg]->fND_RecoEnergy->Clone(name);
00499     
00500     TH1D* fHelperRatio = NULL;
00501     fHelperRatio = (TH1D*) this->GetFNRatio(bg);
00502 
00503     for(int i=0;i<extrapHist->GetNbinsX()+1;i++){
00504       Double_t ratio = fHelperRatio->GetBinContent(i);
00505       extrapHist->SetBinContent(i,extrapHist->GetBinContent(i)*ratio);
00506       extrapHist->SetBinError(i,extrapHist->GetBinError(i)*ratio);
00507       if(i > 10) i += 1000;
00508     }
00509 
00510     delete fHelperRatio;
00511     if(fNearDataPOT == 0) cout<<"No pot"<<endl;
00512     extrapHist->Scale(fTargetPOT/fNearDataPOT);
00513     return extrapHist;
00514   }
00515 
00516   if(bg == Background::kBNueCC)  // Just return the FarMC
00517   {
00518      if(fExtrapMethod == 1){
00519        extrapHist = (TH1D*) fHistMap[bg]->fFD_RecoEnergy->Clone(name);
00520      }else{
00521        TH2D* farHist = (TH2D*) fHistMap[bg]->fFD_RecoVsTrue;
00522        extrapHist = (TH1D*) CreateOscHist(farHist, 12, 12);
00523      }                                                                                                               
00524      extrapHist->Scale(fTargetPOT/fFarPOT);
00525      return extrapHist;
00526   }
00527 
00528    //If nue, nutau use more complicated approach:
00529   if(bg ==Background::kNuTauCC ||  bg ==Background::kNueCC){
00530      extrapHist = (TH1D*) fHistMap[bg]->fFD_RecoEnergy->Clone(name);
00531      extrapHist->Reset("ICE");
00532                                                                                                       
00533      TH1D* trueHist   = (TH1D*) fHistMap[bg]->fFD_TrueEnergy->Clone("TrueE");
00534      trueHist->Reset("ICE");
00535                                                                                                         
00536      TH2D* TrueToReco = (TH2D*) fHistMap[bg]->fFD_RecoVsTrue->Clone("TtoR");
00537      for(int j = 0; j < TrueToReco->GetNbinsY()+1; j++){
00538        double total = 0;
00539        for(int i = 0; i < TrueToReco->GetNbinsX()+1; i++){
00540           total += TrueToReco->GetBinContent(i,j);
00541        }
00542        if(total > 0){
00543          for(int i = 0; i < TrueToReco->GetNbinsX()+1; i++){
00544           TrueToReco->SetBinContent(i,j, TrueToReco->GetBinContent(i,j)/total);
00545          }
00546        }
00547      }
00548                                                                                                         
00549      if(fExtrapMethod == 1) this->BuildAppTrueHistExact(bg, trueHist);
00550      else                   this->BuildAppTrueHistFast(bg, trueHist);
00551 
00552      TH1D* eff = fTauEff;
00553      if(bg ==Background::kNueCC) eff = fMREEff;
00554                                                                                                         
00555      int nBin = eff->GetNbinsX();
00556      double* pos = new double[nBin];
00557      double* Eff = new double[nBin];
00558      for(int k = 0; k < nBin; k++){
00559         pos[k] = eff->GetBinCenter(k);
00560         Eff[k] = eff->GetBinContent(k);
00561      }
00562                                                                                                         
00563      for(int k = 0; k < TrueToReco->GetNbinsY(); k++)  //Bins in TrueGeV
00564      {
00565         double trueEnt = trueHist->GetBinContent(k);
00566         for(int m = 0; m < TrueToReco->GetNbinsX()+1; m++){  //bins in Reco E NueGeV
00567            double recoE = TrueToReco->GetXaxis()->GetBinCenter(m);
00568            if(recoE > 9){ m += 1000; continue; }
00569            double TtoR = TrueToReco->GetBinContent(m,k);
00570            extrapHist->Fill(recoE, trueEnt*TtoR*Eff[m]);
00571            if(recoE != pos[m]) cout<<"Skew on efficiency "<<m<<"  "<<pos[m]<<"  "<<recoE<<endl;
00572         }
00573      } //End of True To Reco
00574 
00575      extrapHist->Scale(fTargetPOT/fFarCCPOT);
00576      delete [] pos;
00577      delete [] Eff;
00578      delete TrueToReco;
00579      delete trueHist;                                                                                                     
00580      return extrapHist;
00581   }
00582      
00583   return NULL;
00584   
00585 }
00586 
00587 void NueExtrapolationJB::AddNueSystematic(NueSystematic *nueSys)
00588 {
00589   fSystematics.push_back(nueSys);
00590 }
00591 
00592 void NueExtrapolationJB::RunSystematicStudy(Selection::Selection_t sel)
00593 {
00594   vector<Selection::Selection_t> temp;
00595   temp.push_back(sel);
00596  
00597   return RunSystematicStudy(temp);
00598 }
00599 
00600 
00601 void NueExtrapolationJB::RunSystematicStudy(vector<Selection::Selection_t> &sel)
00602 {
00603    string fname;
00604    InitializeExtrapHistograms();
00605    InitializePredictionHistograms();
00606    LoadFiles();
00607 
00608    bool isNeugen = false;
00609    for(unsigned int j = 0; j < fSystematics.size(); j++){
00610      Systematic::Systematic_t sys = Systematic::kKNO;
00611      if(fSystematics[j]->GetSysValue(sys) != Systematic::GetDefaultValue(sys))
00612         isNeugen = true;
00613      sys = Systematic::kMA_QE;
00614      if(fSystematics[j]->GetSysValue(sys) != Systematic::GetDefaultValue(sys))
00615         isNeugen = true;
00616      sys = Systematic::kMA_RES;
00617      if(fSystematics[j]->GetSysValue(sys) != Systematic::GetDefaultValue(sys))
00618         isNeugen = true;
00619    }
00620 
00621    if(isNeugen)  InitializeNeugen();
00622 
00623    for(unsigned int i = 0; i < sel.size(); i++){
00624      LoadDataHistograms(fname, sel[i]);
00625 
00626      for(unsigned int j = 0; j < fSystematics.size(); j++){
00627         fCurrentSys = fSystematics[j];
00628         std::cout<<"systematic "<<j<<endl;
00629         PrepareExtrapHistograms(sel[i]);
00630         MakePrediction();
00631         WriteToFile(sel[i]);
00632      }
00633    }
00634 }
00635 
00636 void NueExtrapolationJB::RunExtrapStudy(Selection::Selection_t sel)
00637 {
00638    string fname;
00639    InitializeExtrapHistograms();
00640    InitializePredictionHistograms();
00641    LoadFiles();
00642    
00643    LoadDataHistograms(fname, sel);
00644 
00645    for(unsigned int j = 0; j < fSystematics.size(); j++){
00646       fCurrentSys = fSystematics[j];
00647       if(fExtrapMethod != 3 || j == 0) PrepareExtrapHistograms(sel);
00648       if(fExtrapMethod == 3 && j == 0){
00649         for(unsigned int k = 0; k < fData.size(); k++)   fData[k]->Clear();
00650       }
00651 
00652       MakePrediction();
00653 
00654       double th13, delta, dm23, dum;
00655       int mH;
00656 
00657       fCurrentSys->GetOscParams(dum,dum,th13,dum,dm23,delta,mH);
00658       std::map<Background::Background_t, TH1D*>::iterator mapBeg =  fPredMap.begin();
00659       std::map<Background::Background_t, TH1D*>::iterator mapEnd =  fPredMap.end();
00660 
00661       cout<<th13<<"  "<<delta/3.1415926<<"  "<<dm23<<"  "<<mH<<"  ";
00662 
00663       double nc, numu, nue, bnue, tau;
00664       nc = numu= nue = bnue = tau = 0.0;
00665       double total = 0.0;
00666       while(mapBeg != mapEnd){
00667         if(mapBeg->second){
00668           double sum = mapBeg->second->GetSum();
00669           if(mapBeg->first == Background::kNC) nc = sum;
00670           if(mapBeg->first == Background::kNuMuCC) numu = sum;
00671           if(mapBeg->first == Background::kNuTauCC) tau = sum;
00672           if(mapBeg->first == Background::kNueCC) nue = sum;
00673           if(mapBeg->first == Background::kBNueCC) bnue = sum;
00674           total += sum;
00675         }
00676         mapBeg++;
00677       }
00678      cout<<nc<<"  "<<numu<<"  "<<bnue<<"  "<<tau<<"  "<<nue<<"  "<<total<<endl;
00679    }
00680 }
00681 
00682 
00683 void NueExtrapolationJB::SetOutputFile(string name)
00684 {
00685    outFileName = name;
00686 }
00687 
00688 
00689 void NueExtrapolationJB::WriteToFile(Selection::Selection_t sel)
00690 {
00691   static TFile* file = 0;  
00692   static TTree* tree = 0;
00693   static char selection[256];
00694 
00695   if(file == 0){
00696      file = new TFile(outFileName.c_str(),"RECREATE");
00697      tree = new TTree("energytree","energytree");
00698      tree->Branch("Selection",selection,"Selection/C");
00699      tree->Branch("nearPOT",&fNearPOT,"nearPOT/D");
00700      tree->Branch("farPOT",&fFarPOT,"farPOT/D");
00701      tree->Branch("nearCCPOT",&fNearCCPOT,"nearCCPOT/D");
00702      tree->Branch("farCCPOT",&fFarCCPOT,"farCCPOT/D");
00703   }
00704   file->cd();
00705                                                                                 
00706   sprintf(selection,"%s",Selection::AsString(sel));
00707 
00708   std::map<Background::Background_t,FNHistsM2*>::iterator FNbeg = fHistMap.begin();
00709   std::map<Background::Background_t,FNHistsM2*>::iterator FNend = fHistMap.end();
00710 
00711   while(FNbeg!=FNend) {
00712     string fnh_name = FNbeg->second->fDirectory->GetName();
00713     fnh_name += "_" + string(fCurrentSys->GetName()) + "_" + string(Selection::AsString(sel));
00714 
00715     TDirectory *filedir = file->mkdir(fnh_name.c_str());
00716     filedir->cd();
00717     TList *list = FNbeg->second->fDirectory->GetList();
00718     TIter iter(list->MakeIterator());
00719     TObject *ob = 0;
00720     while((ob = iter())) ob->Write();
00721     file->cd();
00722     FNbeg++;
00723   }
00724 
00725   //Now the predictions
00726   string fnh_name = "Predictions_"+ string(fCurrentSys->GetName()) + "_" + string(Selection::AsString(sel));
00727 
00728   TDirectory* filedir = file->mkdir(fnh_name.c_str());
00729   filedir->cd();
00730   std::map<Background::Background_t, TH1D*>::iterator mapBeg =  fPredMap.begin();
00731   std::map<Background::Background_t, TH1D*>::iterator mapEnd =  fPredMap.end();
00732 
00733   std::map<Background::Background_t, FNHistsM2*>::iterator dataBeg =  fDataMap.begin();
00734   std::map<Background::Background_t, FNHistsM2*>::iterator dataEnd =  fDataMap.end();
00735                                                                                 
00736                                                                                 
00737   while(mapBeg != mapEnd){
00738      if(mapBeg->second) mapBeg->second->Write();
00739      mapBeg++;
00740      dataBeg->second->fND_RecoEnergy->Write();
00741      dataBeg->second->fFD_RecoEnergy->Write();
00742      dataBeg++;
00743   }
00744   file->cd();
00745 
00746   fCurrentSys->MakeBranches(tree);                                                                                
00747   tree->Fill();
00748 
00749   if(fCurrentSys == fSystematics[fSystematics.size()-1])
00750      tree->Write();
00751 }
00752 
00753 void NueExtrapolationJB::InitializeExtrapHistograms()
00754 {
00755   Int_t max_bg_index = 0;
00756 
00757   int nRecoBins;
00758   Double_t* RecoBins;
00759 
00760   while(strcmp(Background::
00761                AsString(Background::EBackground(max_bg_index)),
00762                "?Unknown?")!=0) {
00763     gDirectory->cd("/");
00764     string fnh_name = string(Background::
00765                              AsString(Background::EBackground(max_bg_index)));
00766     FNHistsM2 *fnh = new FNHistsM2(fnh_name.c_str());
00767  
00768     nRecoBins = fNXBins;
00769     RecoBins = fXBins;
00770 
00771     if(Background::EBackground(max_bg_index) == Background::kSelCC){
00772        nRecoBins = fNZBins;
00773        RecoBins  = fZBins;
00774     }
00775 
00776     fnh->fDirectory->cd();
00777     fnh->fND_RecoEnergy = new TH1D("ND_RecoEnergy","ND Reco Energy",nRecoBins, RecoBins);
00778     fnh->fND_RecoEnergy->Sumw2();
00779     fnh->fFD_RecoEnergy = new TH1D("FD_RecoEnergy","FD Reco Energy",nRecoBins, RecoBins);
00780     fnh->fFD_RecoEnergy->Sumw2();
00781     fnh->fND_TrueEnergy = new TH1D("ND_TrueEnergy","ND True Energy",fNYBins,fYBins);
00782     fnh->fND_TrueEnergy->Sumw2();
00783     fnh->fFD_TrueEnergy = new TH1D("FD_TrueEnergy","FD True Energy",fNYBins,fYBins);
00784     fnh->fFD_TrueEnergy->Sumw2();
00785                                                                                 
00786     fnh->fFD_RecoVsTrue = new TH2D("FD_RecoVsTrue", "FD Reco v True E", nRecoBins, RecoBins,
00787                                        fNYBins, fYBins);
00788     fnh->fFD_RecoVsTrue->Sumw2();
00789  
00790     if(Background::EBackground(max_bg_index) == Background::kSelCC){
00791       fnh->fFD_numu_TrueEnergyBase = new TH1D("fFD_numu_TrueEnergyBase","FD True Numu in Fid vs. True Energy",fNYBins,fYBins);
00792       fnh->fFD_numu_TrueEnergyBase->Sumw2();
00793     }
00794 
00795     if(fExtrapMethod != 1 && Background::EBackground(max_bg_index) == Background::kNuMuCC){
00796       fnh->fFD_RecoVsTrue_BNue = new TH2D("FD_RecoVSTrue_BNue","FD BNue to Numu Reco vs. True Energy",nRecoBins, RecoBins,fNYBins,fYBins);
00797       fnh->fFD_RecoVsTrue_BNue->Sumw2();
00798     }
00799 
00800                                                                                 
00801     (fHistMap)[Background::EBackground(max_bg_index)] = fnh;
00802     max_bg_index++;
00803   }
00804 }
00805 
00806 void NueExtrapolationJB::InitializePredictionHistograms()
00807 {
00808   vector<Background::Background_t> bgs;
00809   bgs.push_back(Background::kNC);
00810   bgs.push_back(Background::kNuMuCC);
00811   bgs.push_back(Background::kNuTauCC);
00812   bgs.push_back(Background::kNueCC);
00813   bgs.push_back(Background::kBNueCC);
00814   bgs.push_back(Background::kSelCC);
00815                                                                                                                                  
00816   int nRecoBins;
00817   Double_t* RecoBins;
00818                                                                                                                                  
00819   for(unsigned int i = 0; i < bgs.size(); i++)
00820   {
00821     string fnh_name = "Data Sample " + string(Background::AsString(bgs[i]));
00822 
00823     nRecoBins = fNXBins;
00824     RecoBins = fXBins;
00825      
00826     TH1D* temp = NULL; //new TH1D(fnh_name.c_str(),"Reco Energy",nRecoBins, RecoBins);
00827 //    temp->Sumw2();
00828 
00829     (fPredMap)[bgs[i]] = temp;
00830 
00831     gDirectory->cd("/");
00832 
00833     FNHistsM2 *fnh = new FNHistsM2(fnh_name.c_str());
00834                                                                                 
00835     if(bgs[i] == Background::kSelCC){
00836        nRecoBins = fNZBins;
00837        RecoBins  = fZBins;
00838     }
00839 
00840     TString name = "ND_Data_RecoEnergy_" + string(Background::AsString(bgs[i]));
00841 
00842     fnh->fND_RecoEnergy = new TH1D(name,"ND Reco Energy",nRecoBins, RecoBins);
00843     fnh->fND_RecoEnergy->Sumw2();
00844 
00845     name = "FD_Data_RecoEnergy_" + string(Background::AsString(bgs[i]));
00846     fnh->fFD_RecoEnergy = new TH1D(name,"FD Reco Energy",nRecoBins, RecoBins);
00847     fnh->fFD_RecoEnergy->Sumw2();
00848 
00849     fDataMap[bgs[i]] = fnh;
00850   }
00851 }
00852 
00853 
00854 void NueExtrapolationJB::ResetExtrapHistograms()
00855 {
00856 
00857   if(fHistMap.size() == 0) InitializeExtrapHistograms();
00858 
00859    //Really I should change this to iterate over the fHistMap
00860 
00861   Int_t max_bg_index = 0;
00862   while(strcmp(Background::
00863                AsString(Background::EBackground(max_bg_index)),
00864                "?Unknown?")!=0) {
00865 
00866     Background::Background_t bg = Background::EBackground(max_bg_index);
00867     FNHistsM2* fnh = fHistMap[bg];
00868 
00869     fnh->fDirectory->cd();
00870     fnh->fND_RecoEnergy->Reset("ICE");
00871     fnh->fFD_RecoEnergy->Reset("ICE");
00872     fnh->fND_TrueEnergy->Reset("ICE");
00873     fnh->fFD_TrueEnergy->Reset("ICE");
00874     fnh->fFD_RecoVsTrue->Reset("ICE");
00875 
00876     if(bg == Background::kSelCC) fnh->fFD_numu_TrueEnergyBase->Reset("ICE");
00877     if(bg == Background::kNuMuCC && fExtrapMethod != 1) fnh->fFD_RecoVsTrue_BNue->Reset("ICE");                                                                                                        
00878     max_bg_index++;
00879   }
00880 
00881 }
00882 
00883 TH1 *NueExtrapolationJB::GetFNRatio(Background::Background_t bg)
00884 {
00885   Double_t NearPOT = 0.0;
00886   Double_t FarPOT = 0.0;
00887 
00888   for(unsigned int i = 0; i < fData.size(); i++){
00889      if(fData[i]->IsData() ) continue;  
00890      bool isNue = fData[i]->IsNueData();
00891 
00892      if(bg == Background::kSelCC && isNue) continue;
00893      else if(!isNue) continue;
00894 
00895      Detector::Detector_t det = fData[i]->GetDetector();
00896          
00897      if(det == Detector::kNear) NearPOT = fData[i]->GetPOT();
00898      if(det == Detector::kFar) FarPOT = fData[i]->GetPOT();
00899   }
00900 
00901   TH1 *helperHistFD = (TH1*) fHistMap[bg]->fFD_RecoEnergy;
00902   TH1* farRecoHist = 0;
00903                                                                                                         
00904   if(bg == Background::kNuMuCC && fExtrapMethod != 1)
00905   {
00906      TH2D* farHist = (TH2D*) fHistMap[bg]->fFD_RecoVsTrue;
00907      farRecoHist = this->CreateOscHist(farHist, 14, 14);
00908 
00909      farHist = (TH2D*) fHistMap[bg]->fFD_RecoVsTrue_BNue;
00910      TH1* bnue = this->CreateOscHist(farHist, 12, 14);
00911      farRecoHist->Add(bnue);
00912      delete bnue;
00913                                                                                                         
00914      helperHistFD = farRecoHist;
00915   }
00916                                                                                                         
00917   TH1 *helperHistND = (TH1*) fHistMap[bg]->fND_RecoEnergy;
00918   TH1D* fHelperRatio = (TH1D*) helperHistFD->Clone("helperRatio");
00919   fHelperRatio->SetDirectory(0);
00920   fHelperRatio->Divide(helperHistND);
00921   if(farRecoHist != 0) delete farRecoHist;
00922 
00923   if(FarPOT!=0) fHelperRatio->Scale(NearPOT/FarPOT);
00924   return (TH1*) fHelperRatio;
00925 }
00926 
00927 TH1 *NueExtrapolationJB::GetCCRatio(string hist, Background::Background_t /*bg*/)
00928 {
00929   TH1D *helperHistFD = 0;
00930   TH1D *helperHistND = 0;
00931 
00932   Double_t NearPOT = 1.0;
00933   Double_t FarPOT = 1.0;
00934                                                                                                          
00935   if(hist == "PE"){  // P/E 
00936     helperHistFD = (TH1D*) fHistMap[Background::kSelCC]->fFD_numu_TrueEnergyBase;
00937     helperHistND = (TH1D*) fHistMap[Background::kSelCC]->fFD_TrueEnergy;
00938   }
00939  
00940   if(hist == "NN"){   // NearData / NearMC
00941     helperHistFD = (TH1D*) fDataMap[Background::kSelCC]->fND_RecoEnergy;
00942     helperHistND = (TH1D*) fHistMap[Background::kSelCC]->fND_RecoEnergy;
00943                                                                                                          
00944     for(unsigned int i = 0; i < fData.size(); i++){
00945       if(fData[i]->GetDetector() != Detector::kNear) continue;
00946       if(fData[i]->IsNueData()) continue;
00947                                                                                                          
00948       if(fData[i]->IsData() ) continue;   //Don't need to remake Data Hists
00949                                                                                                          
00950       FarPOT  = fNearCCDataPOT;
00951       if(!fData[i]->IsData())                NearPOT = fData[i]->GetPOT();
00952     }
00953   }
00954 
00955   TH1D* fHelperRatio = (TH1D*) helperHistFD->Clone("helperRatio");
00956   fHelperRatio->SetDirectory(0);
00957   fHelperRatio->Divide(helperHistND);
00958 
00959   if(FarPOT!=0) fHelperRatio->Scale(NearPOT/FarPOT);
00960   return (TH1*) fHelperRatio;
00961 }
00962 
00963 void NueExtrapolationJB::InitializeNeugen()
00964 {
00965   for(unsigned int i = 0; i < fData.size(); i++){
00966     if(fData[i]->IsData() ) continue;   //Don't need to remake Data Hists
00967   
00968     NueHeader nh;
00969     fData[i]->SetupNueHeader(nh);
00970     NueRecord* nr = new NueRecord(nh);
00971 
00972     Systematic::Systematic_t sys = Systematic::kMA_QE;
00973     double defV = Systematic::GetDefaultValue(sys);
00974 
00975     NueSystematic nuesys("temp");
00976 
00977     for(int j = 0; j < fData[i]->NumberOfEntries(); j++)
00978     {
00979       fData[i]->FillRecord(nr, j);
00980   
00981       double xsec = nuesys.DoNeugenCalc(nr, sys, defV);
00982       fData[i]->SetNeugenStdXsec(xsec, j);
00983     }
00984     delete nr;
00985   }
00986 } 
00987 
00988 
00989 TH1* NueExtrapolationJB::CreateOscHist(TH2* hist, int nonOsc, int nuFlavor)
00990 {
00991  //create a NueRecord for later use
00992   NueRecord* nr = 0;
00993   Background::Background_t bg = Background::TranslateFromMC(1, nuFlavor, nonOsc);
00994 
00995   static int TrueEndBin[10] = {fNYBins,fNYBins,fNYBins,fNYBins,fNYBins, 
00996                                fNYBins,fNYBins,fNYBins,fNYBins,fNYBins};
00997   static int RecoEndBin[10] = {10,10,10,10,10, 10,10,10,10,10};
00998 
00999   for(unsigned int i = 0; i < fData.size(); i++){
01000      if(fData[i]->IsData() ) continue;
01001      bool isNue = fData[i]->IsNueData();
01002      if(!isNue) continue;
01003      Detector::Detector_t det = fData[i]->GetDetector();
01004      if(det == Detector::kFar){
01005        NueHeader nh;
01006        fData[i]->SetupNueHeader(nh);
01007        nr = new NueRecord(nh);
01008      }
01009   }
01010  
01011   TH2D* farHist = (TH2D*) hist;
01012   TH1D* extrapHist = (TH1D*) fHistMap[bg]->fFD_RecoEnergy->Clone("farreco");
01013   extrapHist->Reset();                                              
01014 
01015   int tempTrue = 0;
01016   int pos = bg;
01017   if(nonOsc == 12 && nuFlavor == 14) pos = 8;
01018  
01019   for(int i = 0; i < farHist->GetNbinsY()+1; i++){
01020     double E = farHist->GetYaxis()->GetBinCenter(i);
01021     double oscWeight = 1.0;
01022     nr->mctrue.nuEnergy = E;
01023     
01024     nr->mctrue.nonOscNuFlavor = nonOsc;
01025     nr->mctrue.nuFlavor = nuFlavor; 
01026     oscWeight *= fCurrentSys->GetAppearanceWeight(nr, bg);
01027 
01028     if(i > TrueEndBin[pos]) i = farHist->GetNbinsY()+1;
01029 
01030     for(int j = 0; j < farHist->GetNbinsX()+1; j++){
01031        double events = farHist->GetBinContent(j,i);
01032        events *= oscWeight;
01033        double recoE = farHist->GetXaxis()->GetBinCenter(j);
01034        extrapHist->Fill(recoE,events);
01035 
01036        if(events > 0) {tempTrue = i;}
01037        if(j > RecoEndBin[pos]) j = farHist->GetNbinsX()+1;
01038     }
01039   }
01040 
01041   if(TrueEndBin[pos] != tempTrue){
01042      TrueEndBin[pos] = tempTrue; 
01043 //     cout<<"Setting max pos to "<<tempTrue<<" for "<<pos<<endl;
01044   }
01045 
01046   if(nr != 0) delete nr;
01047   return (TH1*) extrapHist;
01048 }
01049 
01050 void NueExtrapolationJB::BuildAppTrueHistExact(Background::Background_t bg, TH1D* trueHist)
01051 {
01052   TH1D* NearONear  = (TH1D*) this->GetCCRatio("NN", Background::kSelCC);
01053   TH1D* POverE     = (TH1D*) this->GetCCRatio("PE", Background::kSelCC);
01054 
01055   trueHist->Reset();
01056   for(unsigned int i = 0; i < fData.size(); i++){
01057     //Only Need the CC Far MC spectrum
01058     if(fData[i]->IsData() ) continue;   //Don't need to remake Data Hists
01059     if(fData[i]->GetDetector() != Detector::kFar) continue;
01060     if(fData[i]->IsNueData()) continue;   //Only loop over the Far CC events
01061 
01062     NueHeader nh;
01063     fData[i]->SetupNueHeader(nh);
01064     NueRecord* nr = new NueRecord(nh);
01065 
01066     for(int j = 0; j < fData[i]->NumberOfEntries(); j++)
01067     {
01068        fData[i]->FillRecord(nr, j);
01069        double totWeight = 1.0;
01070 
01071        int resCode = nr->mctrue.resonanceCode;
01072        if(resCode < 1001 || resCode > 1004)  resCode = 0;
01073 
01074        resCode = 0;  //just using the totals for now
01075 
01076        //Reweight immediately for the x-sec
01077        TH1F *denom =  fNuMuCCXSec[resCode%1000];
01078        TH1F *num =  fNuTauCCXSec[resCode%1000];
01079 
01080        if(bg == Background::kNueCC)  num =  fNueCCXSec[resCode%1000];
01081 
01082        if(nr->mctrue.nonOscNuFlavor < 0){
01083           denom =  fNuMuBarCCXSec[resCode%1000];
01084           num = fNuTauBarCCXSec[resCode%1000];
01085           if(bg == Background::kNueCC) num =  fNueBarCCXSec[resCode%1000];
01086        }
01087 
01088        double trueEnergy = nr->mctrue.nuEnergy;
01089 
01090        for(int k = 0; k < num->GetNbinsX(); k++){
01091           double start = fXSecPos[k];
01092           double end   = fXSecPos[k+1];
01093           if(trueEnergy > start && trueEnergy < end)
01094           {
01095               float low = denom->GetBinContent(k);
01096               float high = num->GetBinContent(k);
01097               if(low > 0) totWeight *= high/low;
01098               k = num->GetNbinsX();
01099           }
01100        }
01101 
01102        // below threshold just move on
01103        if(totWeight == 0) continue;
01104        totWeight *= fCurrentSys->UpdateRecord(nr, Selection::kCC, bg);
01105        if(!NueStandard::PassesCCSelection(nr)) continue;  //Only things passing CC sel
01106        double recoEnergy = this->GetNueEnergy(nr, Selection::kCC);
01107        for(int k = 1; k < NearONear->GetNbinsX(); k++){
01108           if(recoEnergy > fXBins[k-1] && recoEnergy < fXBins[k]){
01109               totWeight *= NearONear->GetBinContent(k);  k += 1000; 
01110           }
01111        }
01112 
01113        for(int k = 1; k < POverE->GetNbinsX(); k++)  //Bins in TrueGeV
01114        {
01115           double start = fYBins[k-1];
01116           double end   = fYBins[k];
01117 
01118           if(trueEnergy > start && trueEnergy < end)
01119           {
01120                totWeight *= POverE->GetBinContent(k);
01121                trueHist->Fill(trueEnergy, totWeight);
01122                k = POverE->GetNbinsX();
01123           }
01124        }
01125     }// End of loop over entries in chain
01126     delete nr;
01127   }// End of loop over all data
01128 
01129   delete NearONear;
01130   delete POverE;
01131 }
01132 
01133 void NueExtrapolationJB::BuildAppTrueHistFast(Background::Background_t bg, TH1D* trueHist)
01134 {
01135    TH1D* NearONear  = (TH1D*) this->GetCCRatio("NN", Background::kSelCC);
01136    TH1D* POverE     = (TH1D*) this->GetCCRatio("PE", Background::kSelCC);
01137 
01138    trueHist->Reset();
01139 
01140    TH2D* RecoToTrue = (TH2D*) fHistMap[Background::kSelCC]->fFD_RecoVsTrue->Clone("RtoT");
01141 
01142    NueRecord* nr = 0;
01143    for(unsigned int i = 0; i < fData.size(); i++){
01144      //Only Need the CC Far MC spectrum
01145      if(fData[i]->IsData() ) continue;   //Don't need to remake Data Hists
01146      if(fData[i]->GetDetector() != Detector::kFar) continue;
01147      if(fData[i]->IsNueData()) continue;   //Only loop over the Far CC events
01148                                                                                                         
01149      NueHeader nh;
01150      fData[i]->SetupNueHeader(nh);
01151      nr = new NueRecord(nh);
01152    }
01153 
01154    for(int i = 0; i < RecoToTrue->GetNbinsX()+1; i++){
01155       double nOverNWeight = NearONear->GetBinContent(i);
01156       if(nOverNWeight < 0) cout<<"negative weight"<<endl;
01157       for(int j = 0; j < RecoToTrue->GetNbinsY()+1; j++){
01158          double trueE = RecoToTrue->GetYaxis()->GetBinCenter(j);
01159          double entries = RecoToTrue->GetBinContent(i,j)*nOverNWeight;
01160          if(entries < 0) cout<<"entries less than zero"<<i<<"  "<<j<<entries<<"  "<<nOverNWeight<<endl;
01161          trueHist->Fill(trueE,entries);
01162       }
01163    }
01164 
01165    Double_t totWeight = 1.0;
01166 
01167    //Reweight immediately for the x-sec
01168    TH1F *denom =  fNuMuCCXSec[0];
01169    TH1F *num =  fNuTauCCXSec[0];
01170    if(bg == Background::kNueCC)  num =  fNueCCXSec[0];
01171 
01172    int xSecInd = 0;
01173 
01174    for(int i = 0; i < trueHist->GetNbinsX()+1; i++){
01175      double trueE = trueHist->GetBinCenter(i);
01176      if(trueE < 0) continue;
01177      totWeight = 1.0;
01178 
01179      //First do the xsec reweighting (which may not be synchronized)
01180      if(trueE > fXSecPos[fNuMuCCXSec[0]->GetNbinsX()]) continue;
01181      while(trueE > fXSecPos[xSecInd+1]) xSecInd++;
01182 
01183      double start = fXSecPos[xSecInd];
01184      double end   = fXSecPos[xSecInd+1];
01185      if(trueE > start && trueE < end)
01186      {
01187         float low = denom->GetBinContent(xSecInd);
01188         float high = num->GetBinContent(xSecInd);
01189         if(low > 0) totWeight *= high/low;
01190      }
01191      //Next add in the oscillation and tau prod weight
01192      nr->mctrue.nuEnergy = trueE;
01193      totWeight *= fCurrentSys->GetAppearanceWeight(nr, bg);
01194      //Finally add the POverE weight
01195      totWeight *= POverE->GetBinContent(i);
01196      trueHist->SetBinContent(i, trueHist->GetBinContent(i)*totWeight);
01197    }// End of loop over all data
01198 
01199                                                                                                        
01200   delete NearONear;
01201   delete POverE;
01202   delete RecoToTrue;
01203   delete nr;
01204 
01205 }