NueExtrapolationJB Class Reference

#include <NueExtrapolationJB.h>

Inheritance diagram for NueExtrapolationJB:
NueExtrapHelper

List of all members.

Public Member Functions

 NueExtrapolationJB ()
virtual ~NueExtrapolationJB ()
void Reset ()
void Initialize ()
void AddChain (TChain *ch, double POT, bool isNue=true)
bool LoadFiles ()
void AddNueSystematic (NueSystematic *)
void SetSeparation (string)
void PrepareExtrapHistograms (Selection::Selection_t)
TH1 * GetPrediction (Background::Background_t bg)
void RunSystematicStudy (Selection::Selection_t sel)
void RunSystematicStudy (vector< Selection::Selection_t > &sel)
void RunExtrapStudy (Selection::Selection_t sel)
void SetOutputFile (string name)
void WriteToFile (Selection::Selection_t sel)
void UseMCAsData (bool in)
void UseMCAsCCData (bool in)
void SetNueRecoBins (int, double, double)
void SetNueRecoBins (int, double *)
void SetTrueBins (int, double, double)
void SetTrueBins (int, double *)
void SetCCRecoBins (int, double, double)
void SetCCRecoBins (int, double *)
void SetDataFileName (string name)
void SetExtrapMethod (int method)

Private Member Functions

void ResetExtrapHistograms ()
void InitializeExtrapHistograms ()
void InitializePredictionHistograms ()
void InitializeNeugen ()
void MakePrediction ()
void MakeDataHistograms (Selection::Selection_t sel)
void LoadDataHistograms (string fname, Selection::Selection_t sel)
void BuildAppTrueHistExact (Background::Background_t bg, TH1D *trueHist)
void BuildAppTrueHistFast (Background::Background_t bg, TH1D *trueHist)
TH1 * CreateOscHist (TH2 *hist, int nonOsc, int nuFlavor)
TH1 * GetCCRatio (string hist, Background::Background_t bg)
TH1 * GetFNRatio (Background::Background_t bg)

Private Attributes

vector< NueData * > fData
std::map< TChain *, int > fChainDataMap
NueSystematicfCurrentSys
vector< NueSystematic * > fSystematics
std::map
< Background::Background_t,
FNHistsM2 * > 
fHistMap
std::map
< Background::Background_t,
FNHistsM2 * > 
fDataMap
std::map
< Background::Background_t,
TH1D * > 
fPredMap
std::string inDataFileName
std::string outFileName
double fTargetPOT
bool fUseMCAsData
bool fUseMCAsCCData
Int_t fNZBins
Double_t * fZBins
double fFarCCPOT
double fNearCCPOT
double fNearDataPOT
double fNearCCDataPOT
TH1F * fNuMuXSec
TH1F * fNuTauXSec
TH1F * fNueXSec
TH1F * fNuMuBarXSec
TH1F * fNuTauBarXSec
TH1F * fNueBarXSec
TH1F * fNuMuCCXSec [5]
TH1F * fNuTauCCXSec [5]
TH1F * fNueCCXSec [5]
TH1F * fNuMuBarCCXSec [5]
TH1F * fNuTauBarCCXSec [5]
TH1F * fNueBarCCXSec [5]
double fXSecPos [2000]
TH1D * fTauEff
TH1D * fMREEff
std::string fSeparation
int fExtrapMethod

Detailed Description

Definition at line 10 of file NueExtrapolationJB.h.


Constructor & Destructor Documentation

NueExtrapolationJB::NueExtrapolationJB (  ) 

Definition at line 13 of file NueExtrapolationJB.cxx.

References fCurrentSys, fExtrapMethod, fFarCCPOT, fNearCCDataPOT, fNearCCPOT, fNearDataPOT, fTargetPOT, fUseMCAsCCData, fUseMCAsData, and outFileName.

00013                                       :
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 }

virtual NueExtrapolationJB::~NueExtrapolationJB (  )  [inline, virtual]

Definition at line 15 of file NueExtrapolationJB.h.

00015 {};


Member Function Documentation

void NueExtrapolationJB::AddChain ( TChain *  ch,
double  POT,
bool  isNue = true 
)

Definition at line 147 of file NueExtrapolationJB.cxx.

References det, NueMini::fBeam, fChainDataMap, fData, NueMini::fDet, fFarCCPOT, NueExtrapHelper::fFarPOT, NueExtrapHelper::fMini, fNearCCPOT, NueExtrapHelper::fNearPOT, NueExtrapHelper::fRecord, NueMini::fRelease, DataUtil::GetBeamType(), NueHeader::GetBeamType(), VldContext::GetDetector(), RecRecordImp< T >::GetHeader(), NueHeader::GetRelease(), RecHeader::GetVldContext(), Detector::kFar, Detector::kNear, ReleaseType::kUnknown, Detector::kUnknown, BeamType::kUnknown, NueExtrapHelper::SetUpNueAnaChain(), and NueExtrapHelper::SetUpNueMiniChain().

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 }

void NueExtrapolationJB::AddNueSystematic ( NueSystematic nueSys  )  [virtual]

Reimplemented from NueExtrapHelper.

Definition at line 587 of file NueExtrapolationJB.cxx.

References fSystematics.

00588 {
00589   fSystematics.push_back(nueSys);
00590 }

void NueExtrapolationJB::BuildAppTrueHistExact ( Background::Background_t  bg,
TH1D *  trueHist 
) [private]

Definition at line 1050 of file NueExtrapolationJB.cxx.

References fCurrentSys, fData, fNueBarCCXSec, fNueCCXSec, fNuMuBarCCXSec, fNuMuCCXSec, fNuTauBarCCXSec, fNuTauCCXSec, NueExtrapHelper::fXBins, fXSecPos, NueExtrapHelper::fYBins, GetCCRatio(), DataUtil::GetDetector(), NueExtrapHelper::GetNueEnergy(), ReleaseType::IsData(), Selection::kCC, Detector::kFar, Background::kNueCC, Background::kSelCC, NueRecord::mctrue, ANtpTruthInfoBeam::nonOscNuFlavor, ANtpTruthInfo::nuEnergy, NueStandard::PassesCCSelection(), ANtpTruthInfoBeam::resonanceCode, and NueSystematic::UpdateRecord().

Referenced by GetPrediction().

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 }

void NueExtrapolationJB::BuildAppTrueHistFast ( Background::Background_t  bg,
TH1D *  trueHist 
) [private]

Definition at line 1133 of file NueExtrapolationJB.cxx.

References fCurrentSys, fData, fHistMap, fNueCCXSec, fNuMuCCXSec, fNuTauCCXSec, fXSecPos, NueSystematic::GetAppearanceWeight(), GetCCRatio(), DataUtil::GetDetector(), ReleaseType::IsData(), Detector::kFar, Background::kNueCC, Background::kSelCC, NueRecord::mctrue, and ANtpTruthInfo::nuEnergy.

Referenced by GetPrediction().

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 }

TH1 * NueExtrapolationJB::CreateOscHist ( TH2 *  hist,
int  nonOsc,
int  nuFlavor 
) [private]

Definition at line 989 of file NueExtrapolationJB.cxx.

References det, fCurrentSys, fData, fHistMap, NueSystematic::GetAppearanceWeight(), ReleaseType::IsData(), Detector::kFar, NueRecord::mctrue, ANtpTruthInfoBeam::nonOscNuFlavor, ANtpTruthInfo::nuEnergy, ANtpTruthInfo::nuFlavor, and Background::TranslateFromMC().

Referenced by GetFNRatio(), and GetPrediction().

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 }

TH1 * NueExtrapolationJB::GetCCRatio ( string  hist,
Background::Background_t  bg 
) [private]

Definition at line 927 of file NueExtrapolationJB.cxx.

References fData, fDataMap, fHistMap, fNearCCDataPOT, DataUtil::GetDetector(), ReleaseType::IsData(), Detector::kNear, and Background::kSelCC.

Referenced by BuildAppTrueHistExact(), and BuildAppTrueHistFast().

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 }

TH1 * NueExtrapolationJB::GetFNRatio ( Background::Background_t  bg  )  [private]

Definition at line 883 of file NueExtrapolationJB.cxx.

References NueConvention::bnue, CreateOscHist(), det, fData, fExtrapMethod, fHistMap, ReleaseType::IsData(), Detector::kFar, Detector::kNear, Background::kNuMuCC, and Background::kSelCC.

Referenced by GetPrediction().

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 }

TH1 * NueExtrapolationJB::GetPrediction ( Background::Background_t  bg  ) 

Definition at line 491 of file NueExtrapolationJB.cxx.

References Background::AsString(), BuildAppTrueHistExact(), BuildAppTrueHistFast(), CreateOscHist(), fDataMap, fExtrapMethod, fFarCCPOT, NueExtrapHelper::fFarPOT, fHistMap, fMREEff, fNearDataPOT, fTargetPOT, fTauEff, GetFNRatio(), Background::kBNueCC, Background::kNC, Background::kNueCC, Background::kNuMuCC, Background::kNuTauCC, Munits::m, and total().

Referenced by MakePrediction().

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 }

void NueExtrapolationJB::Initialize (  ) 

Definition at line 54 of file NueExtrapolationJB.cxx.

References NueExtrapHelper::Init().

00055 {
00056    this->Init();
00057 }

void NueExtrapolationJB::InitializeExtrapHistograms (  )  [private]

Definition at line 753 of file NueExtrapolationJB.cxx.

References bfld::AsString(), FNHistsM2::fDirectory, fExtrapMethod, FNHistsM2::fFD_numu_TrueEnergyBase, FNHistsM2::fFD_RecoEnergy, FNHistsM2::fFD_RecoVsTrue, FNHistsM2::fFD_RecoVsTrue_BNue, FNHistsM2::fFD_TrueEnergy, fHistMap, FNHistsM2::fND_RecoEnergy, FNHistsM2::fND_TrueEnergy, NueExtrapHelper::fXBins, NueExtrapHelper::fYBins, fZBins, Background::kNuMuCC, and Background::kSelCC.

Referenced by ResetExtrapHistograms(), RunExtrapStudy(), and RunSystematicStudy().

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 }

void NueExtrapolationJB::InitializeNeugen (  )  [private]

Definition at line 963 of file NueExtrapolationJB.cxx.

References NueSystematic::DoNeugenCalc(), fData, Systematic::GetDefaultValue(), ReleaseType::IsData(), and Systematic::kMA_QE.

Referenced by RunSystematicStudy().

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 } 

void NueExtrapolationJB::InitializePredictionHistograms (  )  [private]

Definition at line 806 of file NueExtrapolationJB.cxx.

References Background::AsString(), fDataMap, FNHistsM2::fFD_RecoEnergy, FNHistsM2::fND_RecoEnergy, fPredMap, NueExtrapHelper::fXBins, fZBins, Background::kBNueCC, Background::kNC, Background::kNueCC, Background::kNuMuCC, Background::kNuTauCC, and Background::kSelCC.

Referenced by RunExtrapStudy(), and RunSystematicStudy().

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 }

void NueExtrapolationJB::LoadDataHistograms ( string  fname,
Selection::Selection_t  sel 
) [private]

Definition at line 346 of file NueExtrapolationJB.cxx.

References Selection::AsString(), NueConvention::bnue, fDataMap, fMREEff, fSeparation, fTauEff, fUseMCAsData, NueExtrapHelper::fXBins, fZBins, inDataFileName, Background::kBNueCC, Background::kNC, Background::kNuMuCC, Background::kSelCC, and MakeDataHistograms().

Referenced by RunExtrapStudy(), and RunSystematicStudy().

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 }

bool NueExtrapolationJB::LoadFiles (  ) 

Definition at line 59 of file NueExtrapolationJB.cxx.

References fChainDataMap, fData, NueExtrapHelper::fMini, fNearCCDataPOT, fNearDataPOT, fNueBarCCXSec, fNueCCXSec, fNuMuBarCCXSec, fNuMuCCXSec, fNuTauBarCCXSec, fNuTauCCXSec, NueExtrapHelper::fRecord, fUseMCAsData, fXSecPos, DataUtil::GetDetector(), ReleaseType::IsData(), Detector::kFar, NueConvention::NueEnergyCorrection(), NueExtrapHelper::SetUpNueAnaChain(), and NueExtrapHelper::SetUpNueMiniChain().

Referenced by RunExtrapStudy(), and RunSystematicStudy().

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 }

void NueExtrapolationJB::MakeDataHistograms ( Selection::Selection_t  sel  )  [private]

Definition at line 295 of file NueExtrapolationJB.cxx.

References NueSystematic::AddSystematic(), det, MuELoss::e, fData, fDataMap, ANtpTruthInfoBeamNue::fNueClass, Systematic::GetDefaultValue(), NueExtrapHelper::GetNueEnergy(), ReleaseType::IsData(), Selection::kCC, Detector::kFar, Detector::kNear, Systematic::kOscProb, Background::kSelCC, Systematic::kSKZP, NueRecord::mctrue, NueExtrapHelper::PassCuts(), NueStandard::PassesCCSelection(), NueStandard::PassesPIDSelection(), NueStandard::PassesPreSelection(), NueSystematic::SetOscParams(), Background::TranslateFromNueClass(), and NueSystematic::UpdateRecord().

Referenced by LoadDataHistograms().

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 }

void NueExtrapolationJB::MakePrediction (  )  [private]

Definition at line 475 of file NueExtrapolationJB.cxx.

References fPredMap, and GetPrediction().

Referenced by RunExtrapStudy(), and RunSystematicStudy().

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 }

void NueExtrapolationJB::PrepareExtrapHistograms ( Selection::Selection_t  sel  ) 

Definition at line 192 of file NueExtrapolationJB.cxx.

References det, fCurrentSys, fData, fExtrapMethod, fHistMap, ANtpTruthInfoBeamNue::fNueClass, NueExtrapHelper::GetNueEnergy(), NueSystematic::GetSysValue(), ReleaseType::IsData(), Background::kBNueCC, Selection::kCC, Detector::kFar, Detector::kNear, Background::kNueCC, Background::kNuMuCC, Background::kNuTauCC, Systematic::kOscProb, Background::kSelCC, NueRecord::mctrue, ANtpTruthInfoBeam::nonOscNuFlavor, ANtpTruthInfo::nuEnergy, NueExtrapHelper::PassCuts(), NueStandard::PassesCCSelection(), NueStandard::PassesPIDSelection(), NueStandard::PassesPreSelection(), ResetExtrapHistograms(), NueSystematic::SetSysValue(), Background::TranslateFromNueClass(), and NueSystematic::UpdateRecord().

Referenced by RunExtrapStudy(), and RunSystematicStudy().

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 }

void NueExtrapolationJB::Reset (  )  [inline]

Definition at line 17 of file NueExtrapolationJB.h.

00017 {};

void NueExtrapolationJB::ResetExtrapHistograms (  )  [private]

Definition at line 854 of file NueExtrapolationJB.cxx.

References bfld::AsString(), FNHistsM2::fDirectory, fExtrapMethod, FNHistsM2::fFD_numu_TrueEnergyBase, FNHistsM2::fFD_RecoEnergy, FNHistsM2::fFD_RecoVsTrue, FNHistsM2::fFD_RecoVsTrue_BNue, FNHistsM2::fFD_TrueEnergy, fHistMap, FNHistsM2::fND_RecoEnergy, FNHistsM2::fND_TrueEnergy, InitializeExtrapHistograms(), Background::kNuMuCC, and Background::kSelCC.

Referenced by PrepareExtrapHistograms().

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 }

void NueExtrapolationJB::RunExtrapStudy ( Selection::Selection_t  sel  ) 

Definition at line 636 of file NueExtrapolationJB.cxx.

References NueConvention::bnue, fCurrentSys, fData, fExtrapMethod, fname, fPredMap, fSystematics, NueSystematic::GetOscParams(), InitializeExtrapHistograms(), InitializePredictionHistograms(), Background::kBNueCC, Background::kNC, Background::kNueCC, Background::kNuMuCC, Background::kNuTauCC, LoadDataHistograms(), LoadFiles(), MakePrediction(), NueConvention::nue, NueConvention::numu, PrepareExtrapHistograms(), and total().

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 }

void NueExtrapolationJB::RunSystematicStudy ( vector< Selection::Selection_t > &  sel  ) 

Definition at line 601 of file NueExtrapolationJB.cxx.

References fCurrentSys, fname, fSystematics, Systematic::GetDefaultValue(), InitializeExtrapHistograms(), InitializeNeugen(), InitializePredictionHistograms(), Systematic::kKNO, Systematic::kMA_QE, Systematic::kMA_RES, LoadDataHistograms(), LoadFiles(), MakePrediction(), PrepareExtrapHistograms(), and WriteToFile().

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 }

void NueExtrapolationJB::RunSystematicStudy ( Selection::Selection_t  sel  ) 

Definition at line 592 of file NueExtrapolationJB.cxx.

00593 {
00594   vector<Selection::Selection_t> temp;
00595   temp.push_back(sel);
00596  
00597   return RunSystematicStudy(temp);
00598 }

void NueExtrapolationJB::SetCCRecoBins ( int  ,
double *   
)
void NueExtrapolationJB::SetCCRecoBins ( int  ,
double  ,
double   
)
void NueExtrapolationJB::SetDataFileName ( string  name  )  [inline]

Definition at line 126 of file NueExtrapolationJB.h.

References inDataFileName.

00126 {inDataFileName = name;}

void NueExtrapolationJB::SetExtrapMethod ( int  method  )  [inline]

Definition at line 48 of file NueExtrapolationJB.h.

References fExtrapMethod.

00048 { fExtrapMethod = method;};

void NueExtrapolationJB::SetNueRecoBins ( int  ,
double *   
)
void NueExtrapolationJB::SetNueRecoBins ( int  ,
double  ,
double   
)
void NueExtrapolationJB::SetOutputFile ( string  name  ) 

Definition at line 683 of file NueExtrapolationJB.cxx.

References outFileName.

00684 {
00685    outFileName = name;
00686 }

void NueExtrapolationJB::SetSeparation ( string  in  )  [inline]

Definition at line 125 of file NueExtrapolationJB.h.

References fSeparation.

00125 {fSeparation = in;}

void NueExtrapolationJB::SetTrueBins ( int  ,
double *   
)
void NueExtrapolationJB::SetTrueBins ( int  ,
double  ,
double   
)
void NueExtrapolationJB::UseMCAsCCData ( bool  in  )  [inline]

Definition at line 124 of file NueExtrapolationJB.h.

References fUseMCAsCCData.

00124 {fUseMCAsCCData = in; }

void NueExtrapolationJB::UseMCAsData ( bool  in  )  [inline]

Definition at line 123 of file NueExtrapolationJB.h.

References fUseMCAsData.

00123 {fUseMCAsData = in; }

void NueExtrapolationJB::WriteToFile ( Selection::Selection_t  sel  ) 

Definition at line 689 of file NueExtrapolationJB.cxx.

References Selection::AsString(), fCurrentSys, fDataMap, fFarCCPOT, NueExtrapHelper::fFarPOT, fHistMap, fNearCCPOT, NueExtrapHelper::fNearPOT, fPredMap, fSystematics, NueSystematic::GetName(), NueSystematic::MakeBranches(), and outFileName.

Referenced by RunSystematicStudy().

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 }


Member Data Documentation

std::map<TChain*, int> NueExtrapolationJB::fChainDataMap [private]

Definition at line 73 of file NueExtrapolationJB.h.

Referenced by AddChain(), and LoadFiles().

Definition at line 92 of file NueExtrapolationJB.h.

Referenced by AddChain(), GetPrediction(), NueExtrapolationJB(), and WriteToFile().

TH1D* NueExtrapolationJB::fMREEff [private]

Definition at line 116 of file NueExtrapolationJB.h.

Referenced by GetPrediction(), and LoadDataHistograms().

Definition at line 96 of file NueExtrapolationJB.h.

Referenced by GetCCRatio(), LoadFiles(), and NueExtrapolationJB().

Definition at line 93 of file NueExtrapolationJB.h.

Referenced by AddChain(), NueExtrapolationJB(), and WriteToFile().

Definition at line 95 of file NueExtrapolationJB.h.

Referenced by GetPrediction(), LoadFiles(), and NueExtrapolationJB().

Definition at line 111 of file NueExtrapolationJB.h.

Referenced by BuildAppTrueHistExact(), and LoadFiles().

Definition at line 103 of file NueExtrapolationJB.h.

TH1F* NueExtrapolationJB::fNueCCXSec[5] [private]

Definition at line 107 of file NueExtrapolationJB.h.

Referenced by BuildAppTrueHistExact(), BuildAppTrueHistFast(), and LoadFiles().

Definition at line 100 of file NueExtrapolationJB.h.

Definition at line 109 of file NueExtrapolationJB.h.

Referenced by BuildAppTrueHistExact(), and LoadFiles().

Definition at line 101 of file NueExtrapolationJB.h.

TH1F* NueExtrapolationJB::fNuMuCCXSec[5] [private]

Definition at line 105 of file NueExtrapolationJB.h.

Referenced by BuildAppTrueHistExact(), BuildAppTrueHistFast(), and LoadFiles().

Definition at line 98 of file NueExtrapolationJB.h.

Definition at line 110 of file NueExtrapolationJB.h.

Referenced by BuildAppTrueHistExact(), and LoadFiles().

Definition at line 102 of file NueExtrapolationJB.h.

Definition at line 106 of file NueExtrapolationJB.h.

Referenced by BuildAppTrueHistExact(), BuildAppTrueHistFast(), and LoadFiles().

Definition at line 99 of file NueExtrapolationJB.h.

Int_t NueExtrapolationJB::fNZBins [private]

Definition at line 89 of file NueExtrapolationJB.h.

std::string NueExtrapolationJB::fSeparation [private]

Definition at line 117 of file NueExtrapolationJB.h.

Referenced by LoadDataHistograms(), and SetSeparation().

Definition at line 85 of file NueExtrapolationJB.h.

Referenced by GetPrediction(), and NueExtrapolationJB().

TH1D* NueExtrapolationJB::fTauEff [private]

Definition at line 115 of file NueExtrapolationJB.h.

Referenced by GetPrediction(), and LoadDataHistograms().

Definition at line 87 of file NueExtrapolationJB.h.

Referenced by NueExtrapolationJB(), and UseMCAsCCData().

double NueExtrapolationJB::fXSecPos[2000] [private]

Definition at line 113 of file NueExtrapolationJB.h.

Referenced by BuildAppTrueHistExact(), BuildAppTrueHistFast(), and LoadFiles().

Double_t* NueExtrapolationJB::fZBins [private]
std::string NueExtrapolationJB::inDataFileName [private]

Definition at line 83 of file NueExtrapolationJB.h.

Referenced by LoadDataHistograms(), and SetDataFileName().

std::string NueExtrapolationJB::outFileName [private]

Definition at line 84 of file NueExtrapolationJB.h.

Referenced by NueExtrapolationJB(), SetOutputFile(), and WriteToFile().


The documentation for this class was generated from the following files:

Generated on 27 Apr 2017 for loon by  doxygen 1.6.1