ErrorCalc Class Reference

#include <ErrorCalc.h>

Inheritance diagram for ErrorCalc:
ErrorCalc_Joint

List of all members.

Public Member Functions

 ErrorCalc ()
virtual ~ErrorCalc ()
virtual void AddExtrap (Extrapolate2D *E)
void AddFNExtrapSyst (string systname, string farplustag, string nearplustag, string farminustag, string nearminustag, string nominaltag, string filetag, int flag=0, bool useAllRuns=true)
void SetSysFileDir (string dir)
void AddSpecialFNSyst (string systname, double err=0, string histname="NULL", string filetag="NULL", int flag=0)
void AddCovarianceMatrix (string systname, string file, string histname, int flag=0)
void SetNDDataFile (string s)
void CalculateFNExtrapError_SingleBin ()
void CalculateAppearanceExtrapError_SingleBin (Background::Background_t bg)
void CalculateSystErrorMatrix ()
void CalculateHOOError_SingleBin ()
void CalculateMRCCError_SingleBin ()
virtual void CalculateHOOError ()
void FNRatioError (string ratioerrfilename="FNRatioErrs.root")
void NDSystematics (string nderrfilename="NDSysts.root")
void FDSystematics ()
void MakeFracError_Lists ()
virtual void SetGridPred (int nbins, vector< vector< double > > nc, vector< vector< double > > cc, vector< vector< double > > bnue, vector< vector< double > > tau, vector< vector< double > > sig)
virtual void SetUseGrid (bool grid=false)
TH2D * GetShifts (Background::Background_t bg, bool plus1sigma, string systname, int extrap)

Public Attributes

TH2D * CovMatrix
TH2D * CovMatrix_Decomp
std::vector< TH1D * > SysBkgd_Plus1Sigma
std::vector< TH1D * > SysSig_Plus1Sigma

Protected Member Functions

virtual void Initialize ()
void InitializeSys ()
void InitializeDecompSys ()
void ReadSysFiles_FNExtrap (int n)
void ReadSysFiles_Appearance (int n)
virtual void ReadSysFiles_FNExtrap_AllRuns (int n)
virtual void ReadSysFiles_Appearance_AllRuns (int n)
virtual void ReadSpecialFiles (int n)
virtual void ReadCovarianceFiles (int n)
virtual void CalculateSystErrorMatrixGrid ()
virtual void CalculateSystErrorMatrixExtrap ()

Protected Attributes

vector< Extrapolate2D * > Extrap
MultiBinAnaHelper MBH
Bool_t Init
Bool_t InitSys
Bool_t InitDecompSys
string SysFileDir
string NDData_infile
TH2D * NDCovMatrix
vector< string > FNExtrap_FarPlusTag
vector< string > FNExtrap_NearPlusTag
vector< string > FNExtrap_FarMinusTag
vector< string > FNExtrap_NearMinusTag
vector< string > FNExtrap_StdTag
vector< string > FNExtrap_SystName
vector< string > FNExtrap_FileTag
vector< int > FNExtrap_Flag
vector< bool > FNExtrap_UseAllRuns
vector< string > SpecialSyst_SystName
vector< string > SpecialSyst_FileTag
vector< string > SpecialSyst_HistName
vector< double > SpecialSyst_Value
vector< int > SpecialSyst_Flag
vector< string > ExtraCovariance_SystName
vector< string > ExtraCovariance_File
vector< string > ExtraCovariance_HistName
vector< int > ExtraCovariance_Flag
std::map< string, vector< TH2D * > > FN_NC_Plus1Sigma
std::map< string, vector< TH2D * > > FN_NC_Minus1Sigma
std::map< string, vector< TH2D * > > FN_NuMuCC_Plus1Sigma
std::map< string, vector< TH2D * > > FN_NuMuCC_Minus1Sigma
std::map< string, vector< TH2D * > > FN_BNueCC_Plus1Sigma
std::map< string, vector< TH2D * > > FN_BNueCC_Minus1Sigma
std::map< string, vector< TH2D * > > FD_NC_Plus1Sigma
std::map< string, vector< TH2D * > > ND_NC_Plus1Sigma
std::map< string, vector< TH2D * > > FD_NuMuCC_Plus1Sigma
std::map< string, vector< TH2D * > > ND_NuMuCC_Plus1Sigma
std::map< string, vector< TH2D * > > FD_BNueCC_Plus1Sigma
std::map< string, vector< TH2D * > > ND_BNueCC_Plus1Sigma
std::map< string, vector< TH2D * > > FD_NC_Minus1Sigma
std::map< string, vector< TH2D * > > ND_NC_Minus1Sigma
std::map< string, vector< TH2D * > > FD_NuMuCC_Minus1Sigma
std::map< string, vector< TH2D * > > ND_NuMuCC_Minus1Sigma
std::map< string, vector< TH2D * > > FD_BNueCC_Minus1Sigma
std::map< string, vector< TH2D * > > ND_BNueCC_Minus1Sigma
std::map< string, vector< TH1D * > > Pred_CC_Fid_Plus1Sigma
std::map< string, vector< TH1D * > > Pred_CC_Fid_Minus1Sigma
std::map< string, vector< TH2D * > > NuTauCC_MC_Plus1Sigma
std::map< string, vector< TH2D * > > NuTauCC_MC_Minus1Sigma
std::map< string, vector< TH2D * > > NueCC_MC_Plus1Sigma
std::map< string, vector< TH2D * > > NueCC_MC_Minus1Sigma
std::map< string, TH2D * > ExtraCovariance
std::map
< Background::Background_t,
vector< TH1D * > > 
GridPred
bool UseGrid

Detailed Description

Definition at line 21 of file ErrorCalc.h.


Constructor & Destructor Documentation

ErrorCalc::ErrorCalc (  ) 

Definition at line 5 of file ErrorCalc.cxx.

References Init, InitDecompSys, InitSys, and UseGrid.

00006 {
00007   Init = false;
00008   InitSys = false;
00009   InitDecompSys = false;
00010   UseGrid = false;
00011   
00012   return;
00013 }

ErrorCalc::~ErrorCalc (  )  [virtual]

Definition at line 14 of file ErrorCalc.cxx.

00015 {
00016 }


Member Function Documentation

void ErrorCalc::AddCovarianceMatrix ( string  systname,
string  file,
string  histname,
int  flag = 0 
)

Definition at line 56 of file ErrorCalc.cxx.

References ExtraCovariance_File, ExtraCovariance_Flag, ExtraCovariance_HistName, and ExtraCovariance_SystName.

00057 {
00058   if(flag>7)
00059   {
00060     cout<<"Error in AddCovarianceMatrix(): 'flag' should be <=7.  Not adding this covariance matrix."<<endl;
00061     return;
00062   }
00063   
00064   ExtraCovariance_SystName.push_back(systname);
00065   ExtraCovariance_File.push_back(file);
00066   ExtraCovariance_HistName.push_back(histname);
00067   ExtraCovariance_Flag.push_back(flag);
00068   
00069   return;
00070 }

void ErrorCalc::AddExtrap ( Extrapolate2D E  )  [virtual]

Reimplemented in ErrorCalc_Joint.

Definition at line 17 of file ErrorCalc.cxx.

References Extrap.

00018 {
00019   Extrap.push_back(E);
00020   return;
00021 }

void ErrorCalc::AddFNExtrapSyst ( string  systname,
string  farplustag,
string  nearplustag,
string  farminustag,
string  nearminustag,
string  nominaltag,
string  filetag,
int  flag = 0,
bool  useAllRuns = true 
)

Definition at line 22 of file ErrorCalc.cxx.

References FNExtrap_FarMinusTag, FNExtrap_FarPlusTag, FNExtrap_FileTag, FNExtrap_Flag, FNExtrap_NearMinusTag, FNExtrap_NearPlusTag, FNExtrap_StdTag, FNExtrap_SystName, and FNExtrap_UseAllRuns.

00023 {
00024   if(flag>2)
00025   {
00026     cout<<"Error in AddFNExtrapSyst(): 'flag' should be <=2.  Not adding this systematic."<<endl;
00027     return;
00028   }
00029   FNExtrap_SystName.push_back(systname);
00030   FNExtrap_FarPlusTag.push_back(farplustag);
00031   FNExtrap_NearPlusTag.push_back(nearplustag);
00032   FNExtrap_FarMinusTag.push_back(farminustag);
00033   FNExtrap_NearMinusTag.push_back(nearminustag);
00034   FNExtrap_StdTag.push_back(stdtag);
00035   FNExtrap_FileTag.push_back(filetag);
00036   FNExtrap_Flag.push_back(flag);
00037   FNExtrap_UseAllRuns.push_back(useAllRuns);
00038   
00039   return;
00040 }

void ErrorCalc::AddSpecialFNSyst ( string  systname,
double  err = 0,
string  histname = "NULL",
string  filetag = "NULL",
int  flag = 0 
)

Definition at line 41 of file ErrorCalc.cxx.

References SpecialSyst_FileTag, SpecialSyst_Flag, SpecialSyst_HistName, SpecialSyst_SystName, and SpecialSyst_Value.

00042 {
00043   if(flag>4)
00044   {
00045     cout<<"Error in AddSpecialFNSyst(): 'flag' should be <=4.  Not adding this systematic."<<endl;
00046     return;
00047   }
00048   SpecialSyst_SystName.push_back(systname);
00049   SpecialSyst_FileTag.push_back(filetag);
00050   SpecialSyst_HistName.push_back(histname);
00051   SpecialSyst_Value.push_back(err);
00052   SpecialSyst_Flag.push_back(flag);
00053   
00054   return;
00055 }

void ErrorCalc::CalculateAppearanceExtrapError_SingleBin ( Background::Background_t  bg  ) 

Definition at line 2220 of file ErrorCalc.cxx.

References Background::AsString(), MuELoss::e, Extrap, Munits::g, Init, Initialize(), InitializeSys(), Background::kNueCC, Background::kNuTauCC, NueCC_MC_Minus1Sigma, NueCC_MC_Plus1Sigma, NuTauCC_MC_Minus1Sigma, and NuTauCC_MC_Plus1Sigma.

02221 {
02222   if(bg!=Background::kNueCC && bg!=Background::kNuTauCC)
02223   {
02224     cout<<"ErrorCalc::CalculateAppearanceExtrapError() works only for NueCC or NuTauCC.  Quitting..."<<endl;
02225     return;
02226   }
02227   
02228   Initialize();
02229   if(!Init) return;
02230   
02231   InitializeSys();
02232   
02233 //   TH3D *h;
02234   TH2D *g;
02235   double pl=0,mn=0,psum=0,nsum=0,sum=0;
02236   unsigned int ie;
02237   string systname;
02238   
02239   if(bg==Background::kNuTauCC) cout<<"Tau Table"<<endl;
02240   if(bg==Background::kNueCC) cout<<"Nue Table"<<endl;
02241   
02242   cout<<"\\begin{table}[htp]"<<endl;
02243   cout<<"\\caption{}"<<endl;
02244   cout<<"\\centering"<<endl;
02245   cout<<"\\begin{tabular}{c|cc}"<<endl;
02246   cout<<"Syst.&"<<string(Background::AsString(bg))<<"&\\\\"<<endl;
02247   cout<<"\\hline"<<endl;
02248   cout<<"\\hline"<<endl;
02249   
02250   psum=0;nsum=0;
02251   
02252   std::map<string, vector<TH2D*> >::iterator sysiter2;
02253   std::map<string, vector<TH2D*> >::iterator last2;
02254   if(bg==Background::kNueCC)
02255   {
02256     sysiter2 = NueCC_MC_Plus1Sigma.begin();
02257     last2  = NueCC_MC_Plus1Sigma.end();
02258   }
02259   else
02260   {
02261     sysiter2 = NuTauCC_MC_Plus1Sigma.begin();
02262     last2  = NuTauCC_MC_Plus1Sigma.end();
02263   }
02264   
02265   while(sysiter2!=last2)
02266   {
02267     systname = sysiter2->first;
02268   
02269     pl=0;
02270     mn=0;
02271     sum=0;
02272     for(ie=0;ie<Extrap.size();ie++)
02273     {
02274       if(bg==Background::kNuTauCC)
02275       {
02276         g = (TH2D*)Extrap[ie]->Pred[Background::kNuTauCC]->Clone();
02277         g->Multiply(NuTauCC_MC_Plus1Sigma[systname][ie]);
02278         pl += g->Integral();
02279       
02280         g = (TH2D*)Extrap[ie]->Pred[Background::kNuTauCC]->Clone();
02281         g->Multiply(NuTauCC_MC_Minus1Sigma[systname][ie]);
02282         mn += g->Integral();
02283       
02284         sum += Extrap[ie]->Pred[Background::kNuTauCC]->Integral();
02285       }
02286       if(bg==Background::kNueCC)
02287       {
02288         g = (TH2D*)Extrap[ie]->Pred[Background::kNueCC]->Clone();
02289         g->Multiply(NueCC_MC_Plus1Sigma[systname][ie]);
02290         pl += g->Integral();
02291       
02292         g = (TH2D*)Extrap[ie]->Pred[Background::kNueCC]->Clone();
02293         g->Multiply(NueCC_MC_Minus1Sigma[systname][ie]);
02294         mn += g->Integral();
02295       
02296         sum += Extrap[ie]->Pred[Background::kNueCC]->Integral();
02297       }
02298     }
02299     
02300     cout<<systname<<"&";
02301     cout<<100.*pl/sum<<"\\% &";
02302     if(TMath::Abs(pl+mn)<1e-10) cout<<"- \\\\"<<endl;
02303     else cout<<100.*mn/sum<<"\\% \\\\"<<endl;
02304     
02305     if((pl>0 && mn>0) || (pl<0 && mn<0))
02306     {
02307       if(TMath::Abs(pl)>TMath::Abs(mn))
02308       {
02309         mn = -1.*pl;
02310       }
02311       else
02312       {
02313         pl = -1.*mn;
02314       }
02315     }
02316     
02317     if(pl>0) psum+=(pl*pl);
02318     else nsum+=(pl*pl);
02319     if(mn>0) psum+=(mn*mn);
02320     else nsum+=(mn*mn);
02321     
02322     sysiter2++;
02323   }
02324   
02325   //no longer needed due to change in how Pred_CC_Fid_Plus/Minus1Sigma is used
02326 //   std::map<string, vector<TH1D*> >::iterator sysiter3 = Pred_CC_Fid_Plus1Sigma.begin();
02327 //   std::map<string, vector<TH1D*> >::iterator last3 = Pred_CC_Fid_Plus1Sigma.end();
02328 //   
02329 //   int nr,np,nt,it;
02330 //   
02331 //   while(sysiter3!=last3)
02332 //   {
02333 //     systname = sysiter3->first;
02334 //   
02335 //     pl=0;
02336 //     mn=0;
02337 //     sum=0;
02338 //     for(ie=0;ie<Extrap.size();ie++)
02339 //     {
02340 //       nr = Extrap[ie]->GetNReco();
02341 //       np = Extrap[ie]->GetNPID();
02342 //       nt = Extrap[ie]->GetNTrue();
02343 //       
02344 //       if(bg==Background::kNuTauCC)
02345 //       {
02346 //         h = (TH3D*)Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNuTau]->Clone();
02347 //         h->Add(Extrap[ie]->Pred_3D[Extrapolate2D::qBNueToNuTau]);
02348 //         sum += h->Integral();
02349 //       }
02350 //       else
02351 //       {
02352 //         h = (TH3D*)Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNue]->Clone();
02353 //         sum += h->Integral();
02354 //       }
02355 //       
02356 //       for(it=0;it<nt;it++)
02357 //       {
02358 //         pl+=h->Integral(1,np,1,nr,it+1,it+1)*Pred_CC_Fid_Plus1Sigma[systname][ie]->GetBinContent(it+1);
02359 //         
02360 //         mn+=h->Integral(1,np,1,nr,it+1,it+1)*Pred_CC_Fid_Minus1Sigma[systname][ie]->GetBinContent(it+1);
02361 //       }
02362 //     }
02363 //     
02364 //     cout<<systname<<"&";
02365 //     cout<<100.*pl/sum<<"\\% &";
02366 //     if(TMath::Abs(pl+mn)<1e-10) cout<<"- \\\\"<<endl;
02367 //     else cout<<100.*mn/sum<<"\\% \\\\"<<endl;
02368 //     
02369 //     if((pl>0 && mn>0) || (pl<0 && mn<0))
02370 //     {
02371 //       if(TMath::Abs(pl)>TMath::Abs(mn))
02372 //       {
02373 //         mn = -1.*pl;
02374 //       }
02375 //       else
02376 //       {
02377 //         pl = -1.*mn;
02378 //       }
02379 //     }
02380 //     
02381 //     if(pl>0) psum+=(pl*pl);
02382 //     else nsum+=(pl*pl);
02383 //     if(mn>0) psum+=(mn*mn);
02384 //     else nsum+=(mn*mn);
02385 //     
02386 //     sysiter3++;
02387 //   }
02388   
02389   psum = sqrt(psum);
02390   nsum = -1.*sqrt(nsum);
02391   
02392   cout<<"\\hline"<<endl;
02393   cout<<"Total Extrap&"<<100.*psum/sum<<"\\% &"<<100.*nsum/sum<<"\\% \\\\"<<endl;
02394   cout<<"\\end{tabular}"<<endl;
02395   cout<<"\\end{table}"<<endl;
02396   
02397   return;
02398 }

void ErrorCalc::CalculateFNExtrapError_SingleBin (  ) 

Definition at line 2073 of file ErrorCalc.cxx.

References NueConvention::bnue, MuELoss::e, Extrap, FN_BNueCC_Minus1Sigma, FN_BNueCC_Plus1Sigma, FN_NC_Minus1Sigma, FN_NC_Plus1Sigma, FN_NuMuCC_Minus1Sigma, FN_NuMuCC_Plus1Sigma, Init, Initialize(), InitializeSys(), Background::kBNueCC, Background::kNC, and Background::kNuMuCC.

02074 {
02075   Initialize();
02076   if(!Init) return;
02077   
02078   InitializeSys();
02079   
02080   string systname;
02081   
02082   unsigned int ie;
02083   TH2D *ch1,*ch2,*ch3;
02084   double bksum=0,pl=0,mn=0,bnue=0,bpl=0,bmn=0;
02085   double ncpl=0,ncmn=0,ccpl=0,ccmn=0;
02086   double nc=0,cc=0;
02087   double psum=0.,nsum=0.;
02088   
02089   std::map<string, vector<TH2D*> >::iterator systiter = FN_NC_Plus1Sigma.begin();
02090   std::map<string, vector<TH2D*> >::iterator last  = FN_NC_Plus1Sigma.end();
02091   
02092   cout<<"\\begin{table}[htp]"<<endl;
02093   cout<<"\\caption{}"<<endl;
02094   cout<<"\\centering"<<endl;
02095   cout<<"\\begin{tabular}{c|cc|cc|cc|cc}"<<endl;
02096   cout<<"Syst.&NC&&NuMuCC&&BNueCC&&Total&\\\\"<<endl;
02097   cout<<"\\hline"<<endl;
02098   cout<<"\\hline"<<endl;
02099   
02100   while(systiter!=last)
02101   {
02102     systname = systiter->first;
02103     
02104     pl=0;
02105     mn=0;
02106     bpl=0;
02107     bmn=0;
02108     ncpl=0;
02109     ncmn=0;
02110     ccpl=0;
02111     ccmn=0;
02112     bksum=0;
02113     bnue=0;
02114     cc=0;
02115     nc=0;
02116     for(ie=0;ie<Extrap.size();ie++)
02117     {
02118       ch1 = (TH2D*)Extrap[ie]->Pred[Background::kNC]->Clone();
02119       ch1->Multiply(FN_NC_Plus1Sigma[systname][ie]);
02120       ncpl+=ch1->Integral();
02121       ch2 = (TH2D*)Extrap[ie]->Pred[Background::kNuMuCC]->Clone();
02122       ch2->Multiply(FN_NuMuCC_Plus1Sigma[systname][ie]);
02123       ccpl+=ch2->Integral();
02124       ch1->Add(ch2);
02125       if(Extrap[ie]->GetFNforBeamNue())
02126       {
02127         ch3 = (TH2D*)Extrap[ie]->Pred[Background::kBNueCC]->Clone();
02128         ch3->Multiply(FN_BNueCC_Plus1Sigma[systname][ie]);
02129         bpl+=ch3->Integral();
02130         ch1->Add(ch3);
02131       }
02132       pl += ch1->Integral();
02133       
02134       ch1 = (TH2D*)Extrap[ie]->Pred[Background::kNC]->Clone();
02135       ch1->Multiply(FN_NC_Minus1Sigma[systname][ie]);
02136       ncmn+=ch1->Integral();
02137       ch2 = (TH2D*)Extrap[ie]->Pred[Background::kNuMuCC]->Clone();
02138       ch2->Multiply(FN_NuMuCC_Minus1Sigma[systname][ie]);
02139       ccmn+=ch2->Integral();
02140       ch1->Add(ch2);
02141       if(Extrap[ie]->GetFNforBeamNue())
02142       {
02143         ch3 = (TH2D*)Extrap[ie]->Pred[Background::kBNueCC]->Clone();
02144         ch3->Multiply(FN_BNueCC_Minus1Sigma[systname][ie]);
02145         bmn+=ch3->Integral();
02146         ch1->Add(ch3);
02147       }
02148       mn += ch1->Integral();
02149       
02150       bksum+=(Extrap[ie]->Pred[Background::kNC]->Integral()+Extrap[ie]->Pred[Background::kNuMuCC]->Integral());
02151       cc+=Extrap[ie]->Pred[Background::kNuMuCC]->Integral();
02152       nc+=Extrap[ie]->Pred[Background::kNC]->Integral();
02153       if(Extrap[ie]->GetFNforBeamNue())
02154       {
02155         bksum+=Extrap[ie]->Pred[Background::kBNueCC]->Integral();
02156         bnue+=Extrap[ie]->Pred[Background::kBNueCC]->Integral();
02157       }
02158     }
02159     
02160     if(Extrap[0]->GetFNforBeamNue())
02161     {
02162       cout<<systname<<"&";
02163       cout<<100.*ncpl/nc<<"\\% &";
02164       if(TMath::Abs(ncpl+ncmn)<1e-10) cout<<"-&";
02165       else cout<<100.*ncmn/nc<<"\\% &";
02166       cout<<100.*ccpl/cc<<"\\% &";
02167       if(TMath::Abs(ccpl+ccmn)<1e-10) cout<<"-&";
02168       else cout<<100.*ccmn/cc<<"\\% &";
02169       cout<<100.*bpl/bnue<<"\\% &";
02170       if(TMath::Abs(bpl+bmn)<1e-10) cout<<"-&";
02171       else cout<<100.*bmn/bnue<<"\\% &";
02172       cout<<100.*pl/bksum<<"\\% &";
02173       if(TMath::Abs(pl+mn)<1e-10) cout<<"- \\\\"<<endl;
02174       else cout<<100.*mn/bksum<<"\\% \\\\"<<endl;
02175     }
02176     else
02177     {
02178       cout<<systname<<"&";
02179       cout<<100.*ncpl/nc<<"\\% &";
02180       if(TMath::Abs(ncpl+ncmn)<1e-10) cout<<"-&";
02181       else cout<<100.*ncmn/nc<<"\\% &";
02182       cout<<100.*ccpl/cc<<"\\% &";
02183       if(TMath::Abs(ccpl+ccmn)<1e-10) cout<<"-&";
02184       else cout<<100.*ccmn/cc<<"\\% &";
02185       cout<<100.*pl/bksum<<"\\% &";
02186       if(TMath::Abs(pl+mn)<1e-10) cout<<"- \\\\"<<endl;
02187       else cout<<100.*mn/bksum<<"\\% \\\\"<<endl;
02188     }
02189     
02190     if((pl>0 && mn>0) || (pl<0 && mn<0))
02191     {
02192       if(TMath::Abs(pl)>TMath::Abs(mn))
02193       {
02194         mn = -1.*pl;
02195       }
02196       else
02197       {
02198         pl = -1.*mn;
02199       }
02200     }
02201     
02202     if(pl>0) psum+=(pl*pl);
02203     else nsum+=(pl*pl);
02204     if(mn>0) psum+=(mn*mn);
02205     else nsum+=(mn*mn);
02206     
02207     systiter++;
02208   }
02209   
02210   psum = sqrt(psum);
02211   nsum = -1.*sqrt(nsum);
02212   
02213   cout<<"\\hline"<<endl;
02214   cout<<"Total Extrap&&&&&&&"<<100.*psum/bksum<<"\\% &"<<100.*nsum/bksum<<"\\% \\\\"<<endl;
02215   cout<<"\\end{tabular}"<<endl;
02216   cout<<"\\end{table}"<<endl;
02217   
02218   return;
02219 }

void ErrorCalc::CalculateHOOError (  )  [virtual]

Reimplemented in ErrorCalc_Joint.

Definition at line 3770 of file ErrorCalc.cxx.

References CovMatrix_Decomp, Extrap, GridPred, Init, InitDecompSys, Initialize(), InitializeDecompSys(), Background::kBNueCC, Background::kNC, Background::kNuMuCC, Munits::m, n, NDCovMatrix, and UseGrid.

Referenced by NueFit2D::BinLikelihood(), NueFit2D::CalculateErrorMatrixInv(), NueFit2D::RunMultiBinPseudoExpts(), and NueFit2D::RunMultiBinPseudoExpts_MHDeltaFit().

03771 {
03772   //this function assumes the error matrix is normalized to the same POT that the denominator of the F/N ratio is (i.e.1e19)
03773   
03774   //note!!!! the binning convention used in the covariance matrix in the HOOHE files is:
03775   //bin = ienergy + ipid*Nenergy + irun*Nenergy*Npid + icomponent*Nenergy*Npid*Nrun
03776   //the FD prediction bin number convention is bin = ipid + ienergy*Npid
03777   
03778   Initialize();
03779   if(!Init) return;
03780   
03781   InitializeDecompSys();
03782   
03783   if(!InitDecompSys)
03784   {
03785     cout<<"Error in setting up CalculateHOOError().  Quitting..."<<endl;
03786     return;
03787   } 
03788   
03789   int n = NDCovMatrix->GetNbinsX();//dimension of error matrix = pid bins x energy bins x number of runs x number of components
03790   
03791   unsigned int npid = Extrap[0]->GetNPID();//number pid bins
03792   unsigned int nenergy = Extrap[0]->GetNReco();//number of reco energy bins
03793 //   int ncomp = 3;//number of background components
03794   unsigned int nrun = Extrap.size();//number of runs
03795   
03796   unsigned int i;
03797   int j,k,m;
03798   int ipid,ienergy,irun,icomp;
03799   unsigned int ipred;
03800   Background::Background_t bg = Background::kNC;
03801   double temp;
03802   
03803   vector< vector<double> > A;//[npid*nenergy] by [npid*nenergy*nrun*ncomp] matrix of FD predictions to go between the covariance matrix binning to the FD prediction binning
03804   for(i=0;i<npid*nenergy;i++)
03805   {
03806     A.push_back( vector<double>() );
03807     for(j=0;j<n;j++)
03808     {
03809       icomp = int(j/(npid*nenergy*nrun));
03810       irun = int((j - icomp*npid*nenergy*nrun)/(nenergy*npid));
03811       ipid = int((j - icomp*npid*nenergy*nrun - irun*nenergy*npid)/nenergy);
03812       ienergy = j%nenergy;
03813       
03814       ipred = ipid + ienergy*npid;//FD prediction bin corresponding to HOOHE covariance matrix bin j
03815       
03816       if(icomp==0) bg = Background::kNC;
03817       else if(icomp==1) bg = Background::kNuMuCC;
03818       else if(icomp==2) bg = Background::kBNueCC;
03819       
03820       A[i].push_back(0);//start out with zeros
03821       if(ipred==i) //if this FD prediction bin in the matrix A (index i) corresponds to the FD prediction bin from the jth bin of the Decomp covariance matrix (given by ipred)
03822       {
03823         if(UseGrid) A[i][j] = GridPred[bg][irun]->GetBinContent(ipred+1);
03824         else A[i][j] = Extrap[irun]->Pred[bg]->GetBinContent(ipid+1,ienergy+1);
03825       }
03826     }
03827   }
03828   
03829 //   cout<<"(NC+NuMuCC+BNueCC) background prediction:"<<endl;
03830 //   vector<double> F;
03831 //   for(i=0;i<npid*nenergy;i++)
03832 //   {
03833 //     F.push_back(0);
03834 //     for(j=0;j<n;j++)
03835 //     {
03836 //       F[i] += A[i][j];
03837 //     }
03838 //     cout<<F[i]<<endl;
03839 //   }
03840   
03841   CovMatrix_Decomp->Reset();
03842   for(unsigned int i=0;i<npid*nenergy;i++)
03843   {
03844     for(unsigned int j=0;j<npid*nenergy;j++)
03845     {
03846       temp=0;
03847       for(k=0;k<n;k++)
03848       {
03849         for(m=0;m<n;m++)
03850         {
03851           temp+=(NDCovMatrix->GetBinContent(k+1,m+1)*A[i][k]*A[j][m]);
03852         }
03853       }
03854       CovMatrix_Decomp->SetBinContent(i+1,j+1,temp);
03855     }
03856   }
03857   
03858   return;
03859 }

void ErrorCalc::CalculateHOOError_SingleBin (  ) 

Definition at line 2399 of file ErrorCalc.cxx.

References Extrap, gSystem(), Init, Initialize(), Background::kBNueCC, Background::kNC, Background::kNuMuCC, n, and NDData_infile.

02400 {
02401   //this function assumes the error matrix is normalized to the same POT that the denominator of the F/N ratio is (i.e.1e19)
02402   Initialize();
02403   if(!Init) return;
02404   
02405   if(gSystem->AccessPathName(gSystem->ExpandPathName(NDData_infile.c_str())))
02406   {
02407     cout<<"Failure to read ND Data file."<<endl;
02408     return;
02409   }
02410   
02411   TFile *f = new TFile(gSystem->ExpandPathName(NDData_infile.c_str()),"READ");
02412   string name = "CovMatrix_"+Extrap[0]->GetNDDataPID();
02413   TH2D *VN = (TH2D*)f->Get(name.c_str());
02414   
02415   int ip,ir;
02416   unsigned int ie;
02417   int nr,np;
02418   int n = VN->GetNbinsX();
02419   vector<double> A;
02420   
02421   for(ie=0;ie<Extrap.size();ie++)
02422   {
02423     nr = Extrap[ie]->GetNReco();
02424     np = Extrap[ie]->GetNPID();
02425     for(ir=0;ir<nr;ir++)
02426     {
02427       for(ip=0;ip<np;ip++)
02428       {
02429         A.push_back(Extrap[ie]->FNRatio[Background::kNC]->GetBinContent(ip+1,ir+1));
02430       }
02431     }
02432   }
02433   
02434   for(ie=0;ie<Extrap.size();ie++)
02435   {
02436     nr = Extrap[ie]->GetNReco();
02437     np = Extrap[ie]->GetNPID();
02438     for(ir=0;ir<nr;ir++)
02439     {
02440       for(ip=0;ip<np;ip++)
02441       {
02442         A.push_back(Extrap[ie]->FNRatio[Background::kNuMuCC]->GetBinContent(ip+1,ir+1));
02443       }
02444     }
02445   }
02446   
02447   for(ie=0;ie<Extrap.size();ie++)
02448   {
02449     nr = Extrap[ie]->GetNReco();
02450     np = Extrap[ie]->GetNPID();
02451     for(ir=0;ir<nr;ir++)
02452     {
02453       for(ip=0;ip<np;ip++)
02454       {
02455         A.push_back(Extrap[ie]->FNRatio[Background::kBNueCC]->GetBinContent(ip+1,ir+1));
02456       }
02457     }
02458   }
02459   
02460   double VF = 0;
02461   for(int i=0;i<n;i++)
02462   {
02463     for(int j=0;j<n;j++)
02464     {
02465       VF += (VN->GetBinContent(i+1,j+1)*A[i]*A[j]);
02466     }
02467   }
02468   
02469   double bksum=0;
02470   for(ie=0;ie<Extrap.size();ie++)
02471   {
02472     bksum+=(Extrap[ie]->Pred[Background::kNC]->Integral()+Extrap[ie]->Pred[Background::kNuMuCC]->Integral()+Extrap[ie]->Pred[Background::kBNueCC]->Integral());
02473   }
02474   
02475   cout<<"Total syst (NC+NuMuCC+BNueCC) due to HOO Input: "<<bksum<<" +- "<<100.*TMath::Sqrt(VF)/bksum<<"%"<<endl;
02476   
02477   return;
02478 }

void ErrorCalc::CalculateMRCCError_SingleBin (  ) 

Definition at line 2479 of file ErrorCalc.cxx.

References Extrap, Form(), gSystem(), Init, Initialize(), Background::kBNueCC, Background::kNC, Background::kNuMuCC, and NDData_infile.

02480 {
02481   //see DocDB- (Josh's syst note from the first analysis), section 6.2 for the explanation of this calculation.  Correlation between the CC and NC components is taken into account.
02482   
02483   Initialize();
02484   if(!Init) return;
02485   
02486   if(gSystem->AccessPathName(gSystem->ExpandPathName(NDData_infile.c_str())))
02487   {
02488     cout<<"Failure to read ND Data file."<<endl;
02489     return;
02490   }
02491   
02492   TFile *f = new TFile(gSystem->ExpandPathName(NDData_infile.c_str()),"READ");
02493   
02494   string p = Extrap[0]->GetNDDataPID();
02495   unsigned int ie;
02496   TH2D *temp = 0;
02497   
02498   TH2D *NDData_NuMuCC_Errors[3];
02499   TH2D *NDData_BNueCC_Errors[3];
02500   TH2D *NDData_Total[3];
02501   TH2D *NDData_NuMuCC[3];
02502   TH2D *NDData_BNueCC[3];
02503   for(ie=0;ie<3;ie++)
02504   {
02505     NDData_NuMuCC_Errors[ie] = (TH2D*)f->Get(Form("NDData_NuMuCC_Errors_%s_Run%i",p.c_str(),ie+1));
02506     NDData_BNueCC_Errors[ie] = (TH2D*)f->Get(Form("NDData_BNueCC_Errors_%s_Run%i",p.c_str(),ie+1));
02507     NDData_Total[ie] = (TH2D*)f->Get(Form("NDData_Total_%s_Run%i",p.c_str(),ie+1));
02508     NDData_NuMuCC[ie] = (TH2D*)f->Get(Form("NDData_NuMuCC_%s_Run%i",p.c_str(),ie+1));
02509     NDData_BNueCC[ie] = (TH2D*)f->Get(Form("NDData_BNueCC_%s_Run%i",p.c_str(),ie+1));
02510   }
02511   
02512   double numucc_syst=0;
02513   double bnuecc_syst=0;
02514   double total_stat=0;
02515   double numucc_stat=0;
02516   double bnuecc_stat=0;
02517   
02518   int ip,ir;
02519   int np,nr;
02520   double t;
02521   for(ie=0;ie<Extrap.size();ie++)
02522   {
02523     temp = (TH2D*)Extrap[ie]->FNRatio[Background::kNuMuCC]->Clone("temp");
02524     temp->Add(Extrap[ie]->FNRatio[Background::kNC],-1.0);
02525     temp->Multiply(NDData_NuMuCC_Errors[ie]);
02526     numucc_syst += temp->Integral();
02527     
02528     temp = (TH2D*)Extrap[ie]->FNRatio[Background::kBNueCC]->Clone("temp");
02529     temp->Add(Extrap[ie]->FNRatio[Background::kNC],-1.0);
02530     temp->Multiply(NDData_BNueCC_Errors[ie]);
02531     bnuecc_syst += temp->Integral();
02532     
02533     nr = Extrap[ie]->GetNReco();
02534     np = Extrap[ie]->GetNPID();
02535     for(ip=0;ip<np;ip++)
02536     {
02537       for(ir=0;ir<nr;ir++)
02538       {
02539         total_stat += (Extrap[ie]->FNRatio[Background::kNC]->GetBinContent(ip+1,ir+1)*Extrap[ie]->FNRatio[Background::kNC]->GetBinContent(ip+1,ir+1)*NDData_Total[ie]->GetBinError(ip+1,ir+1)*NDData_Total[ie]->GetBinError(ip+1,ir+1));
02540         
02541         t = Extrap[ie]->FNRatio[Background::kNuMuCC]->GetBinContent(ip+1,ir+1) - Extrap[ie]->FNRatio[Background::kNC]->GetBinContent(ip+1,ir+1);
02542         numucc_stat += (t*t*NDData_NuMuCC[ie]->GetBinError(ip+1,ir+1)*NDData_NuMuCC[ie]->GetBinError(ip+1,ir+1));
02543         
02544         t = Extrap[ie]->FNRatio[Background::kBNueCC]->GetBinContent(ip+1,ir+1) - Extrap[ie]->FNRatio[Background::kNC]->GetBinContent(ip+1,ir+1);
02545         bnuecc_stat += (t*t*NDData_BNueCC[ie]->GetBinError(ip+1,ir+1)*NDData_BNueCC[ie]->GetBinError(ip+1,ir+1));
02546       }
02547     }
02548   }
02549   
02550   double totalsyst = sqrt(numucc_syst*numucc_syst + bnuecc_syst*bnuecc_syst + total_stat + numucc_stat + bnuecc_stat);
02551   
02552   double bksum=0;
02553   for(ie=0;ie<Extrap.size();ie++)
02554   {
02555     bksum+=(Extrap[ie]->Pred[Background::kNC]->Integral()+Extrap[ie]->Pred[Background::kNuMuCC]->Integral()+Extrap[ie]->Pred[Background::kBNueCC]->Integral());
02556   }
02557   
02558   cout<<"Total syst (NC+NuMuCC+BNueCC) due to MRCC Input: "<<bksum<<" +- "<<100.*totalsyst/bksum<<"%"<<endl;
02559   
02560   return;
02561 }

void ErrorCalc::CalculateSystErrorMatrix (  ) 
void ErrorCalc::CalculateSystErrorMatrixExtrap (  )  [protected, virtual]

Reimplemented in ErrorCalc_Joint.

Definition at line 3299 of file ErrorCalc.cxx.

References CovMatrix, ExtraCovariance, ExtraCovariance_Flag, Extrap, FD_BNueCC_Plus1Sigma, FN_BNueCC_Plus1Sigma, FN_NC_Plus1Sigma, FN_NuMuCC_Plus1Sigma, gSystem(), Init, Initialize(), InitializeSys(), Background::kBNueCC, Background::kNC, Background::kNueCC, Background::kNuMuCC, Background::kNuTauCC, NDData_infile, NueCC_MC_Plus1Sigma, and NuTauCC_MC_Plus1Sigma.

Referenced by CalculateSystErrorMatrix().

03300 {
03301   //calculates the SYSTEMATIC part of the error matrix
03302   if(Extrap.size()==0)
03303   {
03304     cout<<"Error in ErrorCalc::CalculateSystErrorMatrix(): You must add an Extrapolate2D object.  Quitting..."<<endl;
03305     return;
03306   }
03307   
03308   if(gSystem->AccessPathName(gSystem->ExpandPathName(NDData_infile.c_str())))
03309   {
03310     cout<<"Failure to read ND Data file."<<endl;
03311     return;
03312   }
03313   
03314   Initialize();
03315   if(!Init) return;
03316   
03317   InitializeSys();
03318   
03319   string systname;
03320   
03321   int nPID = Extrap[0]->GetNPID();
03322   int nbins = Extrap[0]->GetNReco()*Extrap[0]->GetNPID();
03323 //   int nt = Extrap[0]->GetNTrue();
03324   
03325   std::map<string, vector<TH2D*> >::iterator systiter = FN_NC_Plus1Sigma.begin();
03326   std::map<string, vector<TH2D*> >::iterator last  = FN_NC_Plus1Sigma.end();
03327   
03328   unsigned int ie;
03329   int i,j;
03330   int ip,ir,jp,jr;
03331   double element;
03332   double di,dj;
03333   
03334   CovMatrix->Reset();
03335   
03336   for(i=0;i<nbins;i++)
03337   {
03338     ir = int(i/nPID);
03339     ip = i%nPID;
03340       
03341     for(j=0;j<nbins;j++)
03342     {
03343       jr = int(j/nPID);
03344       jp = j%nPID;
03345       
03346       element = 0;
03347       systiter = FN_NC_Plus1Sigma.begin();
03348       while(systiter!=last)
03349       {
03350         systname = systiter->first;
03351         
03352         di=0;
03353         dj=0;
03354         for(ie=0;ie<Extrap.size();ie++)
03355         {
03356           di += FN_NC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->Pred[Background::kNC]->GetBinContent(ip+1,ir+1) + FN_NuMuCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->Pred[Background::kNuMuCC]->GetBinContent(ip+1,ir+1) + NuTauCC_MC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->Pred[Background::kNuTauCC]->GetBinContent(ip+1,ir+1) + NueCC_MC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->Pred[Background::kNueCC]->GetBinContent(ip+1,ir+1);
03357           
03358           dj += FN_NC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*Extrap[ie]->Pred[Background::kNC]->GetBinContent(jp+1,jr+1) + FN_NuMuCC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*Extrap[ie]->Pred[Background::kNuMuCC]->GetBinContent(jp+1,jr+1) + NuTauCC_MC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*Extrap[ie]->Pred[Background::kNuTauCC]->GetBinContent(jp+1,jr+1) + NueCC_MC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*Extrap[ie]->Pred[Background::kNueCC]->GetBinContent(jp+1,jr+1);
03359           
03360           if(Extrap[0]->GetFNforBeamNue())
03361           {
03362           di += FN_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->Pred[Background::kBNueCC]->GetBinContent(ip+1,ir+1);
03363           
03364           dj += FN_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*Extrap[ie]->Pred[Background::kBNueCC]->GetBinContent(jp+1,jr+1);
03365             
03366           }
03367           else
03368           {
03369             di += FD_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->Pred[Background::kBNueCC]->GetBinContent(ip+1,ir+1);
03370             
03371             dj += FD_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*Extrap[ie]->Pred[Background::kBNueCC]->GetBinContent(jp+1,jr+1);
03372             
03373              }
03374         }
03375         
03376         element += (di*dj);
03377         
03378         systiter++;
03379       }
03380       
03381       CovMatrix->SetBinContent(i+1,j+1,element);
03382     }
03383   }
03384   
03385   std::map<string,TH2D*>::iterator coviter = ExtraCovariance.begin();
03386   std::map<string,TH2D*>::iterator covlast  = ExtraCovariance.end();
03387   double ni,nj;
03388   int k=0;
03389   int flag;
03390   
03391   k=0;
03392   while(coviter!=covlast)
03393   {
03394     systname = coviter->first;
03395     flag = ExtraCovariance_Flag.at(k);
03396     for(i=0;i<nbins;i++)
03397     {
03398       ir = int(i/nPID);
03399       ip = i%nPID;
03400       
03401       ni=0;
03402       for(ie=0;ie<Extrap.size();ie++)
03403       {
03404         if(flag==0 || flag==1 || flag==4 || flag==5)
03405         {
03406           ni+=Extrap[ie]->Pred[Background::kNC]->GetBinContent(ip+1,ir+1);
03407         }
03408         if(flag==0 || flag==1 || flag==4 || flag==6)
03409         {
03410           ni+=Extrap[ie]->Pred[Background::kNuMuCC]->GetBinContent(ip+1,ir+1);
03411         }
03412         if(flag==0 || flag==1 || flag==7)
03413         {
03414           ni+=Extrap[ie]->Pred[Background::kBNueCC]->GetBinContent(ip+1,ir+1);
03415         }
03416         if(flag==0 || flag==2)
03417         {
03418           ni+=Extrap[ie]->Pred[Background::kNuTauCC]->GetBinContent(ip+1,ir+1);
03419         }
03420         if(flag==0 || flag==3)
03421         {
03422           ni+=Extrap[ie]->Pred[Background::kNueCC]->GetBinContent(ip+1,ir+1);
03423         }
03424       }
03425       
03426       for(j=0;j<nbins;j++)
03427       {
03428         jr = int(j/nPID);
03429         jp = j%nPID;
03430         
03431         nj=0;
03432         for(ie=0;ie<Extrap.size();ie++)
03433         {
03434           if(flag==0 || flag==1 || flag==4 || flag==5)
03435           {
03436             nj+=Extrap[ie]->Pred[Background::kNC]->GetBinContent(jp+1,jr+1);
03437           }
03438           if(flag==0 || flag==1 || flag==4 || flag==6)
03439           {
03440             nj+=Extrap[ie]->Pred[Background::kNuMuCC]->GetBinContent(jp+1,jr+1);
03441           }
03442           if(flag==0 || flag==1 || flag==7)
03443           {
03444             nj+=Extrap[ie]->Pred[Background::kBNueCC]->GetBinContent(jp+1,jr+1);
03445           }
03446           if(flag==0 || flag==2)
03447           {
03448             nj+=Extrap[ie]->Pred[Background::kNuTauCC]->GetBinContent(jp+1,jr+1);
03449           }
03450           if(flag==0 || flag==3)
03451           {
03452             nj+=Extrap[ie]->Pred[Background::kNueCC]->GetBinContent(jp+1,jr+1);
03453           }
03454         }
03455         element = CovMatrix->GetBinContent(i+1,j+1);
03456         element += (ni*nj*ExtraCovariance[systname]->GetBinContent(i+1,j+1));
03457         CovMatrix->SetBinContent(i+1,j+1,element);
03458       }
03459     }
03460     
03461     coviter++;
03462     k++;
03463   }
03464   
03465   return;
03466 }

void ErrorCalc::CalculateSystErrorMatrixGrid (  )  [protected, virtual]

Reimplemented in ErrorCalc_Joint.

Definition at line 3467 of file ErrorCalc.cxx.

References CovMatrix, ExtraCovariance, ExtraCovariance_Flag, Extrap, FD_BNueCC_Plus1Sigma, FN_BNueCC_Plus1Sigma, FN_NC_Plus1Sigma, FN_NuMuCC_Plus1Sigma, GridPred, gSystem(), Init, Initialize(), InitializeSys(), Background::kBNueCC, Background::kNC, Background::kNueCC, Background::kNuMuCC, Background::kNuTauCC, NDData_infile, NueCC_MC_Plus1Sigma, and NuTauCC_MC_Plus1Sigma.

Referenced by CalculateSystErrorMatrix().

03468 {
03469   //calculates the SYSTEMATIC part of the error matrix
03470   if(Extrap.size()==0)
03471   {
03472     cout<<"Error in ErrorCalc::CalculateSystErrorMatrixGrid(): You must add an Extrapolate2D object.  Quitting..."<<endl;
03473     return;
03474   }
03475   
03476   if(gSystem->AccessPathName(gSystem->ExpandPathName(NDData_infile.c_str())))
03477   {
03478     cout<<"Failure to read ND Data file."<<endl;
03479     return;
03480   }
03481   
03482   Initialize();
03483   if(!Init) return;
03484   
03485   InitializeSys();
03486   
03487   string systname;
03488   
03489   int nPID = Extrap[0]->GetNPID();
03490   int nbins = Extrap[0]->GetNReco()*Extrap[0]->GetNPID();
03491   
03492   std::map<string, vector<TH2D*> >::iterator systiter = FN_NC_Plus1Sigma.begin();
03493   std::map<string, vector<TH2D*> >::iterator last  = FN_NC_Plus1Sigma.end();
03494   
03495   unsigned int ie;
03496   int i,j;
03497   int ip,ir,jp,jr;
03498   double element;
03499   double di,dj;
03500   
03501   CovMatrix->Reset();
03502   
03503   for(i=0;i<nbins;i++)
03504   {
03505     ir = int(i/nPID);
03506     ip = i%nPID;
03507     
03508     for(j=0;j<nbins;j++)
03509     {
03510       jr = int(j/nPID);
03511       jp = j%nPID;
03512       
03513       element = 0;
03514       systiter = FN_NC_Plus1Sigma.begin();
03515       while(systiter!=last)
03516       {
03517         systname = systiter->first;
03518         
03519         di=0;
03520         dj=0;
03521         for(ie=0;ie<Extrap.size();ie++)
03522         {
03523           di += FN_NC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*GridPred[Background::kNC][ie]->GetBinContent(i+1) + FN_NuMuCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*GridPred[Background::kNuMuCC][ie]->GetBinContent(i+1) + NuTauCC_MC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*GridPred[Background::kNuTauCC][ie]->GetBinContent(i+1) + NueCC_MC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*GridPred[Background::kNueCC][ie]->GetBinContent(i+1);
03524           
03525           dj += FN_NC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*GridPred[Background::kNC][ie]->GetBinContent(j+1) + FN_NuMuCC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*GridPred[Background::kNuMuCC][ie]->GetBinContent(j+1) + NuTauCC_MC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*GridPred[Background::kNuTauCC][ie]->GetBinContent(j+1) + NueCC_MC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*GridPred[Background::kNueCC][ie]->GetBinContent(j+1);
03526           
03527           if(Extrap[0]->GetFNforBeamNue())
03528           {
03529             di += FN_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*GridPred[Background::kBNueCC][ie]->GetBinContent(i+1);
03530             
03531             dj += FN_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*GridPred[Background::kBNueCC][ie]->GetBinContent(j+1);
03532             
03533           }
03534           else
03535           {
03536             di += FD_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*GridPred[Background::kBNueCC][ie]->GetBinContent(i+1);
03537             
03538             dj += FD_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*GridPred[Background::kBNueCC][ie]->GetBinContent(j+1);
03539             
03540           }
03541         }
03542         
03543         element += (di*dj);
03544         
03545         systiter++;
03546       }
03547       
03548       CovMatrix->SetBinContent(i+1,j+1,element);
03549     }
03550   }
03551   
03552   std::map<string,TH2D*>::iterator coviter = ExtraCovariance.begin();
03553   std::map<string,TH2D*>::iterator covlast  = ExtraCovariance.end();
03554   double ni,nj;
03555   int k=0;
03556   int flag;
03557   
03558   k=0;
03559   while(coviter!=covlast)
03560   {
03561     systname = coviter->first;
03562     flag = ExtraCovariance_Flag.at(k);
03563     for(i=0;i<nbins;i++)
03564     {
03565       
03566       ni=0;
03567       for(ie=0;ie<Extrap.size();ie++)
03568       {
03569         if(flag==0 || flag==1 || flag==4 || flag==5)
03570         {
03571           ni+=GridPred[Background::kNC][ie]->GetBinContent(i+1);
03572         }
03573         if(flag==0 || flag==1 || flag==4 || flag==6)
03574         {
03575           ni+=GridPred[Background::kNuMuCC][ie]->GetBinContent(i+1);
03576         }
03577         if(flag==0 || flag==1 || flag==7)
03578         {
03579           ni+=GridPred[Background::kBNueCC][ie]->GetBinContent(i+1);
03580         }
03581         if(flag==0 || flag==2)
03582         {
03583           ni+=GridPred[Background::kNuTauCC][ie]->GetBinContent(i+1);
03584         }
03585         if(flag==0 || flag==3)
03586         {
03587           ni+=GridPred[Background::kNueCC][ie]->GetBinContent(i+1);
03588         }
03589       }
03590       
03591       for(j=0;j<nbins;j++)
03592       {
03593         nj=0;
03594         for(ie=0;ie<Extrap.size();ie++)
03595         {
03596           if(flag==0 || flag==1 || flag==4 || flag==5)
03597           {
03598             nj+=GridPred[Background::kNC][ie]->GetBinContent(j+1);
03599           }
03600           if(flag==0 || flag==1 || flag==4 || flag==6)
03601           {
03602             nj+=GridPred[Background::kNuMuCC][ie]->GetBinContent(j+1);
03603           }
03604           if(flag==0 || flag==1 || flag==7)
03605           {
03606             nj+=GridPred[Background::kBNueCC][ie]->GetBinContent(j+1);
03607           }
03608           if(flag==0 || flag==2)
03609           {
03610             nj+=GridPred[Background::kNuTauCC][ie]->GetBinContent(j+1);
03611           }
03612           if(flag==0 || flag==3)
03613           {
03614             nj+=GridPred[Background::kNueCC][ie]->GetBinContent(j+1);
03615           }
03616         }
03617         element = CovMatrix->GetBinContent(i+1,j+1);
03618         element += (ni*nj*ExtraCovariance[systname]->GetBinContent(i+1,j+1));
03619         CovMatrix->SetBinContent(i+1,j+1,element);
03620       }
03621     }
03622     
03623     coviter++;
03624     k++;
03625   }
03626   
03627   return;
03628 }

void ErrorCalc::FDSystematics (  ) 

Definition at line 3075 of file ErrorCalc.cxx.

References NueConvention::bnue, MuELoss::e, Extrap, FD_BNueCC_Minus1Sigma, FD_BNueCC_Plus1Sigma, FD_NC_Minus1Sigma, FD_NC_Plus1Sigma, FD_NuMuCC_Minus1Sigma, FD_NuMuCC_Plus1Sigma, Init, Initialize(), InitializeSys(), Background::kBNueCC, Background::kNC, and Background::kNuMuCC.

03076 {
03077   Initialize();
03078   if(!Init) return;
03079   
03080   InitializeSys();
03081   
03082   int i;
03083   unsigned int ie;
03084   
03085   int nr = Extrap[0]->GetNReco();
03086   int np = Extrap[0]->GetNPID();
03087   double *r = new double[nr+1];
03088   double *p = new double[np+1];
03089   for(i=0;i<nr+1;i++)
03090   {
03091     r[i] = Extrap[0]->Pred_TotalBkgd->GetYaxis()->GetBinLowEdge(i+1);
03092   }
03093   for(i=0;i<np+1;i++)
03094   {
03095     p[i] = Extrap[0]->Pred_TotalBkgd->GetXaxis()->GetBinLowEdge(i+1);
03096   }
03097   
03098   string systname;
03099   double pl=0,mn=0;
03100   double bksum=0,cc=0,nc=0,bnue=0;
03101   double ccpl=0,ncpl=0,bpl=0;
03102   double ccmn=0,ncmn=0,bmn=0;
03103   double psum=0,nsum=0;
03104   double ccpsum=0,ccnsum=0;
03105   double ncpsum=0,ncnsum=0;
03106   double bpsum=0,bnsum=0;
03107   TH2D *ch1,*ch2,*ch3;
03108   
03109   std::map<string, vector<TH2D*> >::iterator systiter = FD_NC_Plus1Sigma.begin();
03110   std::map<string, vector<TH2D*> >::iterator last  = FD_NC_Plus1Sigma.end();
03111   
03112   cout<<"FD Syst Table"<<endl;
03113   
03114   ofstream myfile("Tables/fdtable.tex");
03115   myfile.precision(2);
03116   myfile<<fixed;
03117   
03118   myfile<<"\\begin{table}[htp]"<<endl;
03119   myfile<<"\\caption{"<<Extrap[0]->GetPID()<<", Far Detector MC}"<<endl;
03120   myfile<<"\\centering"<<endl;
03121   myfile<<"\\begin{tabular}{c|cc|cc|cc|cc}"<<endl;
03122   myfile<<"Syst.&NC&&NuMuCC&&BNueCC&&Total&\\\\"<<endl;
03123   myfile<<"\\hline"<<endl;
03124   myfile<<"\\hline"<<endl;
03125   
03126   systiter = FD_NC_Plus1Sigma.begin();
03127   psum=0;
03128   nsum=0;
03129   ccpsum=0;
03130   ccnsum=0;
03131   ncpsum=0;
03132   ncnsum=0;
03133   bpsum=0;
03134   bnsum=0;
03135   while(systiter!=last)
03136   {
03137     systname = systiter->first;
03138     
03139     pl=0;
03140     mn=0;
03141     bpl=0;
03142     bmn=0;
03143     ncpl=0;
03144     ncmn=0;
03145     ccpl=0;
03146     ccmn=0;
03147     bksum=0;
03148     bnue=0;
03149     cc=0;
03150     nc=0;
03151     for(ie=0;ie<Extrap.size();ie++)
03152     {
03153       ch1 = (TH2D*)Extrap[ie]->FDMC[Background::kNC]->Clone();
03154       ch1->Multiply(FD_NC_Plus1Sigma[systname][ie]);
03155       ncpl+=ch1->Integral();
03156       ch2 = (TH2D*)Extrap[ie]->FDMC[Background::kNuMuCC]->Clone();
03157       ch2->Multiply(FD_NuMuCC_Plus1Sigma[systname][ie]);
03158       ccpl+=ch2->Integral();
03159       ch1->Add(ch2);
03160       ch3 = (TH2D*)Extrap[ie]->FDMC[Background::kBNueCC]->Clone();
03161       ch3->Multiply(FD_BNueCC_Plus1Sigma[systname][ie]);
03162       bpl+=ch3->Integral();
03163       ch1->Add(ch3);
03164       pl += ch1->Integral();
03165       
03166       ch1 = (TH2D*)Extrap[ie]->FDMC[Background::kNC]->Clone();
03167       ch1->Multiply(FD_NC_Minus1Sigma[systname][ie]);
03168       ncmn+=ch1->Integral();
03169       ch2 = (TH2D*)Extrap[ie]->FDMC[Background::kNuMuCC]->Clone();
03170       ch2->Multiply(FD_NuMuCC_Minus1Sigma[systname][ie]);
03171       ccmn+=ch2->Integral();
03172       ch1->Add(ch2);
03173       ch3 = (TH2D*)Extrap[ie]->FDMC[Background::kBNueCC]->Clone();
03174       ch3->Multiply(FD_BNueCC_Minus1Sigma[systname][ie]);
03175       bmn+=ch3->Integral();
03176       ch1->Add(ch3);
03177       mn += ch1->Integral();
03178       
03179       bksum+=(Extrap[ie]->FDMC[Background::kNuMuCC]->Integral() + Extrap[ie]->FDMC[Background::kNC]->Integral() + Extrap[ie]->FDMC[Background::kBNueCC]->Integral());
03180       cc+=Extrap[ie]->FDMC[Background::kNuMuCC]->Integral();
03181       nc+=Extrap[ie]->FDMC[Background::kNC]->Integral();
03182       bnue+=Extrap[ie]->FDMC[Background::kBNueCC]->Integral();
03183     }
03184     
03185     myfile<<systname<<"&";
03186     myfile<<100.*ncpl/nc<<"\\% &";
03187     if(TMath::Abs(ncpl+ncmn)<1e-20) myfile<<"-&";
03188     else myfile<<100.*ncmn/nc<<"\\% &";
03189     myfile<<100.*ccpl/cc<<"\\% &";
03190     if(TMath::Abs(ccpl+ccmn)<1e-20) myfile<<"-&";
03191     else myfile<<100.*ccmn/cc<<"\\% &";
03192     myfile<<100.*bpl/bnue<<"\\% &";
03193     if(TMath::Abs(bpl+bmn)<1e-20) myfile<<"-&";
03194     else myfile<<100.*bmn/bnue<<"\\% &";
03195     myfile<<100.*pl/bksum<<"\\% &";
03196     if(TMath::Abs(pl+mn)<1e-20) myfile<<"- \\\\"<<endl;
03197     else myfile<<100.*mn/bksum<<"\\% \\\\"<<endl;
03198     
03199     if((pl>0 && mn>0) || (pl<0 && mn<0))
03200     {
03201       if(TMath::Abs(pl)>TMath::Abs(mn))
03202       {
03203         mn = -1.*pl;
03204       }
03205       else
03206       {
03207         pl = -1.*mn;
03208       }
03209     }
03210     if((ccpl>0 && ccmn>0) || (ccpl<0 && ccmn<0))
03211     {
03212       if(TMath::Abs(ccpl)>TMath::Abs(ccmn))
03213       {
03214         ccmn = -1.*ccpl;
03215       }
03216       else
03217       {
03218         ccpl = -1.*ccmn;
03219       }
03220     }
03221     if((ncpl>0 && ncmn>0) || (ncpl<0 && ncmn<0))
03222     {
03223       if(TMath::Abs(ncpl)>TMath::Abs(ncmn))
03224       {
03225         ncmn = -1.*ncpl;
03226       }
03227       else
03228       {
03229         ncpl = -1.*ncmn;
03230       }
03231     }
03232     if((bpl>0 && bmn>0) || (bpl<0 && bmn<0))
03233     {
03234       if(TMath::Abs(bpl)>TMath::Abs(bmn))
03235       {
03236         bmn = -1.*bpl;
03237       }
03238       else
03239       {
03240         bpl = -1.*bmn;
03241       }
03242     }
03243     
03244     if(pl>0) psum+=(pl*pl);
03245     else nsum+=(pl*pl);
03246     if(mn>0) psum+=(mn*mn);
03247     else nsum+=(mn*mn);
03248     
03249     if(ccpl>0) ccpsum+=(ccpl*ccpl);
03250     else ccnsum+=(ccpl*ccpl);
03251     if(ccmn>0) ccpsum+=(ccmn*ccmn);
03252     else ccnsum+=(ccmn*ccmn);
03253     
03254     if(ncpl>0) ncpsum+=(ncpl*ncpl);
03255     else ncnsum+=(ncpl*ncpl);
03256     if(ncmn>0) ncpsum+=(ncmn*ncmn);
03257     else ncnsum+=(ncmn*ncmn);
03258     
03259     if(bpl>0) bpsum+=(bpl*bpl);
03260     else bnsum+=(bpl*bpl);
03261     if(bmn>0) bpsum+=(bmn*bmn);
03262     else bnsum+=(bmn*bmn);
03263     
03264     systiter++;
03265   }
03266   
03267   psum = sqrt(psum);
03268   nsum = -1.*sqrt(nsum);
03269   
03270   ccpsum = sqrt(ccpsum);
03271   ccnsum = -1.*sqrt(ccnsum);
03272   
03273   ncpsum = sqrt(ncpsum);
03274   ncnsum = -1.*sqrt(ncnsum);
03275   
03276   bpsum = sqrt(bpsum);
03277   bnsum = -1.*sqrt(bnsum);
03278   
03279   myfile<<"\\hline"<<endl;
03280   myfile<<"Total Extrap&"<<100.*ncpsum/nc<<"&"<<100.*ncnsum/nc<<"&"<<100.*ccpsum/cc<<"&"<<100.*ccnsum/cc<<"&"<<100.*bpsum/bnue<<"&"<<100.*bnsum/bnue<<"&"<<100.*psum/bksum<<"\\% &"<<100.*nsum/bksum<<"\\% \\\\"<<endl;
03281   myfile<<"\\end{tabular}"<<endl;
03282   myfile<<"\\end{table}"<<endl;
03283   myfile.close();
03284   
03285   return;
03286 }

void ErrorCalc::FNRatioError ( string  ratioerrfilename = "FNRatioErrs.root"  ) 

Definition at line 2562 of file ErrorCalc.cxx.

References Extrap, FN_BNueCC_Minus1Sigma, FN_BNueCC_Plus1Sigma, FN_NC_Minus1Sigma, FN_NC_Plus1Sigma, FN_NuMuCC_Minus1Sigma, FN_NuMuCC_Plus1Sigma, Form(), gSystem(), Init, Initialize(), and InitializeSys().

02563 {
02564   Initialize();
02565   if(!Init) return;
02566   
02567   InitializeSys();
02568   
02569   vector<TH2D*> FN_NC_Plus_Total;
02570   vector<TH2D*> FN_NC_Minus_Total;
02571   vector<TH2D*> FN_NuMuCC_Plus_Total;
02572   vector<TH2D*> FN_NuMuCC_Minus_Total;
02573   vector<TH2D*> FN_BNueCC_Plus_Total;
02574   vector<TH2D*> FN_BNueCC_Minus_Total;
02575   
02576   int i,ip,ir;
02577   unsigned int ie;
02578   
02579   int nr = Extrap[0]->GetNReco();
02580   int np = Extrap[0]->GetNPID();
02581   double *r = new double[nr+1];
02582   double *p = new double[np+1];
02583   for(i=0;i<nr+1;i++)
02584   {
02585     r[i] = Extrap[0]->Pred_TotalBkgd->GetYaxis()->GetBinLowEdge(i+1);
02586   }
02587   for(i=0;i<np+1;i++)
02588   {
02589     p[i] = Extrap[0]->Pred_TotalBkgd->GetXaxis()->GetBinLowEdge(i+1);
02590   }
02591   
02592   for(ie=0;ie<Extrap.size();ie++)
02593   {
02594     FN_NC_Plus_Total.push_back(new TH2D(Form("FN_NC_Plus_Total_Run%i",ie+1),"",np,p,nr,r));
02595     FN_NC_Minus_Total.push_back(new TH2D(Form("FN_NC_Minus_Total_Run%i",ie+1),"",np,p,nr,r));
02596     FN_NuMuCC_Plus_Total.push_back(new TH2D(Form("FN_NuMuCC_Plus_Total_Run%i",ie+1),"",np,p,nr,r));
02597     FN_NuMuCC_Minus_Total.push_back(new TH2D(Form("FN_NuMuCC_Minus_Total_Run%i",ie+1),"",np,p,nr,r));
02598     FN_BNueCC_Plus_Total.push_back(new TH2D(Form("FN_BNueCC_Plus_Total_Run%i",ie+1),"",np,p,nr,r));
02599     FN_BNueCC_Minus_Total.push_back(new TH2D(Form("FN_BNueCC_Minus_Total_Run%i",ie+1),"",np,p,nr,r));
02600   }
02601   
02602   string systname;
02603   double pl,mn,temp;
02604   
02605   std::map<string, vector<TH2D*> >::iterator systiter = FN_NC_Plus1Sigma.begin();
02606   std::map<string, vector<TH2D*> >::iterator last  = FN_NC_Plus1Sigma.end();
02607   
02608   for(ie=0;ie<Extrap.size();ie++)
02609   {
02610     for(ip=0;ip<np;ip++)
02611     {
02612       for(ir=0;ir<nr;ir++)
02613       {
02614         pl=0;
02615         mn=0;
02616         systiter = FN_NC_Plus1Sigma.begin();
02617         while(systiter!=last)
02618         {
02619           systname = systiter->first;
02620           
02621           temp = FN_NC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1);
02622           if(temp>0) pl+=temp*temp;
02623           else mn+=temp*temp;
02624           
02625           temp = FN_NC_Minus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1);
02626           if(temp>0) pl+=temp*temp;
02627           else mn+=temp*temp;
02628           
02629           systiter++;
02630         }
02631         FN_NC_Plus_Total[ie]->SetBinContent(ip+1,ir+1,sqrt(pl));
02632         FN_NC_Minus_Total[ie]->SetBinContent(ip+1,ir+1,-1.*sqrt(mn));
02633         
02634         pl=0;
02635         mn=0;
02636         systiter = FN_NC_Plus1Sigma.begin();
02637         while(systiter!=last)
02638         {
02639           systname = systiter->first;
02640           
02641           temp = FN_NuMuCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1);
02642           if(temp>0) pl+=temp*temp;
02643           else mn+=temp*temp;
02644           
02645           temp = FN_NuMuCC_Minus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1);
02646           if(temp>0) pl+=temp*temp;
02647           else mn+=temp*temp;
02648           
02649           systiter++;
02650         }
02651         FN_NuMuCC_Plus_Total[ie]->SetBinContent(ip+1,ir+1,sqrt(pl));
02652         FN_NuMuCC_Minus_Total[ie]->SetBinContent(ip+1,ir+1,-1.*sqrt(mn));
02653         
02654         pl=0;
02655         mn=0;
02656         systiter = FN_NC_Plus1Sigma.begin();
02657         while(systiter!=last)
02658         {
02659           systname = systiter->first;
02660           
02661           temp = FN_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1);
02662           if(temp>0) pl+=temp*temp;
02663           else mn+=temp*temp;
02664           
02665           temp = FN_BNueCC_Minus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1);
02666           if(temp>0) pl+=temp*temp;
02667           else mn+=temp*temp;
02668           
02669           systiter++;
02670         }
02671         FN_BNueCC_Plus_Total[ie]->SetBinContent(ip+1,ir+1,sqrt(pl));
02672         FN_BNueCC_Minus_Total[ie]->SetBinContent(ip+1,ir+1,-1.*sqrt(mn));
02673       }
02674     }
02675   }
02676   
02677   TFile *f = new TFile(gSystem->ExpandPathName(ratioerrfilename.c_str()),"RECREATE");
02678   for(ie=0;ie<Extrap.size();ie++)
02679   {
02680     FN_NC_Plus_Total[ie]->Write();
02681     FN_NC_Minus_Total[ie]->Write();
02682     FN_NuMuCC_Plus_Total[ie]->Write();
02683     FN_NuMuCC_Minus_Total[ie]->Write();
02684     FN_BNueCC_Plus_Total[ie]->Write();
02685     FN_BNueCC_Minus_Total[ie]->Write();
02686   }
02687   f->Close();
02688   
02689   return;
02690 }

TH2D * ErrorCalc::GetShifts ( Background::Background_t  bg,
bool  plus1sigma,
string  systname,
int  extrap 
)

Definition at line 227 of file ErrorCalc.cxx.

References Extrap, FN_BNueCC_Minus1Sigma, FN_BNueCC_Plus1Sigma, FN_NC_Minus1Sigma, FN_NC_Plus1Sigma, FN_NuMuCC_Minus1Sigma, FN_NuMuCC_Plus1Sigma, Init, Initialize(), InitializeSys(), Background::kBNueCC, Background::kNC, Background::kNueCC, Background::kNuMuCC, Background::kNuTauCC, NueCC_MC_Minus1Sigma, NueCC_MC_Plus1Sigma, NuTauCC_MC_Minus1Sigma, and NuTauCC_MC_Plus1Sigma.

00228 {
00229   TH2D *hist = new TH2D("hist","hist",1,0,1,1,0,1);
00230   Initialize();
00231   if(!Init) return hist;
00232   
00233   InitializeSys();
00234   
00235   if(FN_NC_Plus1Sigma.count(systname)==0)
00236   {
00237     cout<<"Unrecognized systname."<<endl;
00238     return hist;
00239   }
00240   if(unsigned(extrap)>=Extrap.size())
00241   {
00242     cout<<"Unrecognized extrap number"<<endl;
00243     return hist;
00244   }
00245   
00246   if(bg==Background::kNC)
00247   {
00248     if(plus1sigma==true)
00249     {
00250       return FN_NC_Plus1Sigma[systname][extrap];
00251     }
00252     else
00253     {
00254       return FN_NC_Minus1Sigma[systname][extrap];
00255     }
00256   }
00257   else if(bg==Background::kNuMuCC)
00258   {
00259     if(plus1sigma==true)
00260     {
00261       return FN_NuMuCC_Plus1Sigma[systname][extrap];
00262     }
00263     else
00264     {
00265       return FN_NuMuCC_Minus1Sigma[systname][extrap];
00266     }
00267   }
00268   else if(bg==Background::kBNueCC)
00269   {
00270     if(plus1sigma==true)
00271     {
00272       return FN_BNueCC_Plus1Sigma[systname][extrap];
00273     }
00274     else
00275     {
00276       return FN_BNueCC_Minus1Sigma[systname][extrap];
00277     }
00278   }
00279   else if(bg==Background::kNuTauCC)
00280   {
00281     if(plus1sigma==true)
00282     {
00283       return NuTauCC_MC_Plus1Sigma[systname][extrap];
00284     }
00285     else
00286     {
00287       return NuTauCC_MC_Minus1Sigma[systname][extrap];
00288     }
00289   }
00290   else if(bg==Background::kNueCC)
00291   {
00292     if(plus1sigma==true)
00293     {
00294       return NueCC_MC_Plus1Sigma[systname][extrap];
00295     }
00296     else
00297     {
00298       return NueCC_MC_Minus1Sigma[systname][extrap];
00299     }
00300   }
00301   else
00302   {
00303     cout<<"Unknown background type."<<endl;
00304     return hist;
00305   }
00306 }

void ErrorCalc::Initialize (  )  [protected, virtual]

Reimplemented in ErrorCalc_Joint.

Definition at line 71 of file ErrorCalc.cxx.

References Background::AsString(), CovMatrix, CovMatrix_Decomp, Extrap, Form(), GridPred, Init, Background::kBNueCC, Background::kNC, Background::kNueCC, Background::kNuMuCC, and Background::kNuTauCC.

Referenced by CalculateAppearanceExtrapError_SingleBin(), CalculateFNExtrapError_SingleBin(), CalculateHOOError(), CalculateHOOError_SingleBin(), CalculateMRCCError_SingleBin(), CalculateSystErrorMatrixExtrap(), CalculateSystErrorMatrixGrid(), FDSystematics(), FNRatioError(), GetShifts(), MakeFracError_Lists(), NDSystematics(), and SetGridPred().

00072 {
00073   if(Init) return;
00074   
00075   cout.precision(5);
00076   cout<<fixed;
00077   
00078   unsigned int ie;
00079   for(ie=0;ie<Extrap.size();ie++)
00080   {
00081     Extrap[ie]->GetPrediction();
00082     if(Extrap[ie]->GetOscFlag())
00083     {
00084       cout<<"Error in ErrorCalc::Initialize(): You have called Extrapolate2D::OscillatePrediction() before initializing the systematics!  Quitting..."<<endl;
00085       return;
00086     }
00087   }
00088   
00089   gROOT->cd();
00090   int nbins = Extrap[0]->GetNReco()*Extrap[0]->GetNPID();
00091   CovMatrix = new TH2D("CovMatrix","",nbins,-0.5,nbins-0.5,nbins,-0.5,nbins-0.5);
00092   CovMatrix_Decomp = new TH2D("CovMatrix_Decomp","",nbins,-0.5,nbins-0.5,nbins,-0.5,nbins-0.5);
00093   
00094   string st;
00095   for(ie=0;ie<Extrap.size();ie++)
00096   {
00097     st = string(Background::AsString(Background::kNC));  
00098     GridPred[Background::kNC].push_back(new TH1D(Form("GridPred_%s_%i",st.c_str(),ie),"",nbins,-0.5,nbins-0.5));
00099     st = string(Background::AsString(Background::kNuMuCC));  
00100     GridPred[Background::kNuMuCC].push_back(new TH1D(Form("GridPred_%s_%i",st.c_str(),ie),"",nbins,-0.5,nbins-0.5));
00101     st = string(Background::AsString(Background::kBNueCC));  
00102     GridPred[Background::kBNueCC].push_back(new TH1D(Form("GridPred_%s_%i",st.c_str(),ie),"",nbins,-0.5,nbins-0.5));
00103     st = string(Background::AsString(Background::kNuTauCC));  
00104     GridPred[Background::kNuTauCC].push_back(new TH1D(Form("GridPred_%s_%i",st.c_str(),ie),"",nbins,-0.5,nbins-0.5));
00105     st = string(Background::AsString(Background::kNueCC));  
00106     GridPred[Background::kNueCC].push_back(new TH1D(Form("GridPred_%s_%i",st.c_str(),ie),"",nbins,-0.5,nbins-0.5));
00107   }
00108   
00109   Init = true;
00110   
00111   return;
00112 }

void ErrorCalc::InitializeDecompSys (  )  [protected]

Definition at line 143 of file ErrorCalc.cxx.

References Extrap, gSystem(), InitDecompSys, Background::kBNueCC, Background::kNC, Background::kNuMuCC, n, NDCovMatrix, and NDData_infile.

Referenced by CalculateHOOError().

00144 {
00145   if(InitDecompSys) return;
00146   
00147   NDCovMatrix=0;
00148   
00149   if(gSystem->AccessPathName(gSystem->ExpandPathName(NDData_infile.c_str())))
00150   {
00151     cout<<"Error in InitializeDecompSys(): Failure to read ND Data file."<<endl;
00152     return;
00153   }
00154   
00155   if(Extrap.size()==0) 
00156   {
00157     cout<<"Error in InitializeDecompSys(): No Extrapolate2D objects have been added."<<endl;
00158     return;
00159   }
00160   
00161   unsigned int npid = Extrap[0]->GetNPID();//number pid bins
00162   unsigned int nenergy = Extrap[0]->GetNReco();//number of reco energy bins
00163   int ncomp = 3;//number of background components
00164   unsigned int nrun = Extrap.size();//number of runs
00165   int ntot=npid*nenergy*ncomp*nrun;
00166   
00167   TFile *f = new TFile(gSystem->ExpandPathName(NDData_infile.c_str()),"READ");
00168   string name = "CovMatrix_"+Extrap[0]->GetNDDataPID();
00169   
00170   if(f->Read(name.c_str())==0)
00171   {
00172     cout<<"Warning in InitializeDecompSys(): Can't read ND covariance matrix from file.  Setting to zero..."<<endl;
00173     NDCovMatrix = new TH2D("NDCovMatrix","",ntot,-0.5,ntot-0.5,ntot,-0.5,ntot-0.5);
00174     InitDecompSys=true;
00175     gROOT->cd();
00176     return;
00177   }
00178   
00179   NDCovMatrix = (TH2D*)f->Get(name.c_str());
00180   
00181   int n = NDCovMatrix->GetNbinsX();//dimension of error matrix = pid bins x energy bins x number of runs x number of components
00182   if(n!=ntot)
00183   {
00184     cout<<"Error in InitializeDecompSys(): The number of elements in the ND covariance matrix doesn't seem to match your setup, i.e. # PID bins x # energy bins x # runs x 3(NC,CC,Bnue)."<<endl;
00185     return;
00186   }
00187   
00188   int i,j;
00189   int ipid,ienergy,irun,icomp;
00190   
00191   Background::Background_t bg = Background::kNC;
00192   vector<double> N;//vector of near detector data with same binning as covariance matrix
00193   for(i=0;i<n;i++)
00194   {
00195     icomp = int(i/(npid*nenergy*nrun));
00196     irun = int((i - icomp*npid*nenergy*nrun)/(nenergy*npid));
00197     ipid = int((i - icomp*npid*nenergy*nrun - irun*nenergy*npid)/nenergy);
00198     ienergy = i%nenergy;
00199     
00200     if(icomp==0) bg = Background::kNC;
00201     else if(icomp==1) bg = Background::kNuMuCC;
00202     else if(icomp==2) bg = Background::kBNueCC;
00203     
00204     N.push_back(Extrap[irun]->NDData[bg]->GetBinContent(ipid+1,ienergy+1));
00205   }
00206   
00207   double temp;
00208   
00209   //make the ND covariance matrix a fractional covariance - this way we don't need to do (F/N)(Delta NData) to get the uncertainty, we do the equivalent thing: FPred*(Delta NData/NData) = (F/N)(NData)*(Delta NData/NData) = (F/N)(Delta NData)
00210   //this way it can be used with grid predictions for FC (where we have only the predictions, not F/N ratio)
00211   for(i=0;i<n;i++)
00212   {
00213     for(j=0;j<n;j++)
00214     {
00215       if (N[i]!=0&&N[j]!=0) {temp = NDCovMatrix->GetBinContent(i+1,j+1)/(N[i]*N[j]);}
00216       else temp = 0;
00217       NDCovMatrix->SetBinContent(i+1,j+1,temp);
00218     }
00219   }
00220   
00221   InitDecompSys=true;
00222   
00223   gROOT->cd();
00224   
00225   return;
00226 }

void ErrorCalc::InitializeSys (  )  [protected]
void ErrorCalc::MakeFracError_Lists (  ) 

Definition at line 3629 of file ErrorCalc.cxx.

References Extrap, FD_BNueCC_Plus1Sigma, FN_BNueCC_Plus1Sigma, FN_NC_Plus1Sigma, FN_NuMuCC_Plus1Sigma, Form(), gSystem(), Init, Initialize(), InitializeSys(), Background::kBNueCC, Background::kNC, Background::kNueCC, Background::kNuMuCC, Background::kNuTauCC, NDData_infile, NueCC_MC_Plus1Sigma, NuTauCC_MC_Plus1Sigma, SysBkgd_Plus1Sigma, and SysSig_Plus1Sigma.

03630 {
03631   //calculates Fractional Errors for the Standard Likelihood
03632   if(Extrap.size()==0)
03633   {
03634     cout<<"Error in ErrorCalc::CalculateSystErrorMatrix(): You must add an Extrapolate2D object.  Quitting..."<<endl;
03635     return;
03636   }
03637   
03638   if(gSystem->AccessPathName(gSystem->ExpandPathName(NDData_infile.c_str())))
03639   {
03640     cout<<"Failure to read ND Data file."<<endl;
03641     return;
03642   }
03643 
03644   Initialize();
03645   if(!Init) return;
03646   
03647   InitializeSys();
03648 
03649   string systname;
03650   
03651   int nPID = Extrap[0]->GetNPID();
03652   int nbins = Extrap[0]->GetNReco()*Extrap[0]->GetNPID();
03653 //   int nt = Extrap[0]->GetNTrue();
03654   
03655   std::map<string, vector<TH2D*> >::iterator systiter = FN_NC_Plus1Sigma.begin();
03656   std::map<string, vector<TH2D*> >::iterator last  = FN_NC_Plus1Sigma.end();
03657   
03658   unsigned int ie;
03659   int i;
03660   unsigned int j;
03661   int ip,ir;
03662   double element_bkg;
03663   double element_sig;
03664   double tot_bkg, tot_sig;
03665   double dnc, dnumu, dbnue;
03666   double dnue, dtau;
03667   double *bins = new double[nbins+1];;
03668 
03669   //binning for new histograms:
03670 
03671   for (int k=0; k<nbins+1; k++){
03672     bins[k] = Extrap[0]->Pred_TotalBkgd_VsBinNumber->GetBinLowEdge(k+1);
03673   }
03674 
03675  //May or may not use this:
03676   //SysBkgd_Plus1Sigma.clear();
03677   //SysSig_Plus1Sigma.clear();
03678 
03679   i=0;
03680   j=0;
03681 
03682   systiter = FN_NC_Plus1Sigma.begin();
03683   while(systiter!=last)
03684     {
03685       systname = systiter->first;
03686 
03687       if (j>=SysBkgd_Plus1Sigma.size()){
03688       SysBkgd_Plus1Sigma.push_back(new TH1D(Form("bkg_sys_%s",systname.c_str()),"",nbins,bins));
03689       } else SysBkgd_Plus1Sigma[j]->Reset();
03690 
03691       if (j>=SysSig_Plus1Sigma.size()){
03692       SysSig_Plus1Sigma.push_back(new TH1D(Form("sig_sys_%s",systname.c_str()),"",nbins,bins));
03693       } else SysSig_Plus1Sigma[j]->Reset();
03694 
03695       if (SysBkgd_Plus1Sigma.size() != SysSig_Plus1Sigma.size()){
03696         cout << "ERROR: Got out of phase somehow in FracErr Calc." << endl;
03697         return;
03698       }
03699 
03700   
03701       for(i=0;i<nbins;i++)
03702         {
03703           ir = int(i/nPID);
03704           ip = i%nPID;
03705       
03706       
03707           element_sig = 0;
03708           element_bkg = 0;
03709 
03710           dnc=0;        
03711           dnumu=0;        
03712           dbnue=0;        
03713           dnue=0;
03714           dtau=0;
03715 
03716           tot_bkg=0;
03717           tot_sig=0;
03718 
03719           element_sig=0;
03720           element_bkg=0;
03721 
03722 
03723           for(ie=0;ie<Extrap.size();ie++)
03724             {
03725               dnc += FN_NC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->Pred[Background::kNC]->GetBinContent(ip+1,ir+1);
03726               dnumu += FN_NuMuCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->Pred[Background::kNuMuCC]->GetBinContent(ip+1,ir+1);
03727               dtau += NuTauCC_MC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->Pred[Background::kNuTauCC]->GetBinContent(ip+1,ir+1);
03728               dnue += NueCC_MC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->Pred[Background::kNueCC]->GetBinContent(ip+1,ir+1);
03729               
03730           
03731               if(Extrap[0]->GetFNforBeamNue())
03732                 {
03733                 dbnue += FN_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->Pred[Background::kBNueCC]->GetBinContent(ip+1,ir+1);
03734                 
03735                }
03736               else
03737               {
03738                  dbnue += FD_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->Pred[Background::kBNueCC]->GetBinContent(ip+1,ir+1);
03739                   
03740               }
03741 
03742               //Calculate total prediction in that in bin, too:       
03743               tot_bkg += Extrap[ie]->Pred_TotalBkgd_VsBinNumber->GetBinContent(i+1);
03744               tot_sig += Extrap[ie]->Pred_Signal_VsBinNumber->GetBinContent(i+1);
03745             }       
03746         
03747           if (tot_sig>0){
03748             element_sig = dnue/tot_sig;
03749           } else element_sig = 0;
03750 
03751           if (tot_bkg>0){
03752             element_bkg = (dnc + dnumu + dtau + dbnue)/tot_bkg;
03753           } else element_bkg = 0;
03754 
03755           SysBkgd_Plus1Sigma[j]->SetBinContent(i+1,element_bkg);
03756           SysSig_Plus1Sigma[j]->SetBinContent(i+1,element_sig);
03757 
03758         }
03759 
03760       j++;
03761       systiter++;
03762       
03763 
03764     }
03765 
03766   
03767   return;
03768 }

void ErrorCalc::NDSystematics ( string  nderrfilename = "NDSysts.root"  ) 

Definition at line 2691 of file ErrorCalc.cxx.

References NueConvention::bnue, MuELoss::e, Extrap, Form(), gSystem(), Init, Initialize(), InitializeSys(), Background::kBNueCC, Background::kNC, Background::kNuMuCC, ND_BNueCC_Minus1Sigma, ND_BNueCC_Plus1Sigma, ND_NC_Minus1Sigma, ND_NC_Plus1Sigma, ND_NuMuCC_Minus1Sigma, and ND_NuMuCC_Plus1Sigma.

02692 {
02693   Initialize();
02694   if(!Init) return;
02695   
02696   InitializeSys();
02697   
02698   vector<TH2D*> ND_NC_Plus_Total;
02699   vector<TH2D*> ND_NC_Minus_Total;
02700   vector<TH2D*> ND_NuMuCC_Plus_Total;
02701   vector<TH2D*> ND_NuMuCC_Minus_Total;
02702   vector<TH2D*> ND_BNueCC_Plus_Total;
02703   vector<TH2D*> ND_BNueCC_Minus_Total;
02704   vector<TH2D*> ND_Total_Plus_Total;
02705   vector<TH2D*> ND_Total_Minus_Total;
02706   TH2D *ND_Total_Plus_Total_AllRuns;
02707   TH2D *ND_Total_Minus_Total_AllRuns;
02708   
02709   int i,ip,ir;
02710   unsigned int ie;
02711   
02712   int nr = Extrap[0]->GetNReco();
02713   int np = Extrap[0]->GetNPID();
02714   double *r = new double[nr+1];
02715   double *p = new double[np+1];
02716   for(i=0;i<nr+1;i++)
02717   {
02718     r[i] = Extrap[0]->Pred_TotalBkgd->GetYaxis()->GetBinLowEdge(i+1);
02719   }
02720   for(i=0;i<np+1;i++)
02721   {
02722     p[i] = Extrap[0]->Pred_TotalBkgd->GetXaxis()->GetBinLowEdge(i+1);
02723   }
02724   
02725   for(ie=0;ie<Extrap.size();ie++)
02726   {
02727     ND_NC_Plus_Total.push_back(new TH2D(Form("ND_NC_Plus_Total_Run%i",ie+1),"",np,p,nr,r));
02728     ND_NC_Minus_Total.push_back(new TH2D(Form("ND_NC_Minus_Total_Run%i",ie+1),"",np,p,nr,r));
02729     ND_NuMuCC_Plus_Total.push_back(new TH2D(Form("ND_NuMuCC_Plus_Total_Run%i",ie+1),"",np,p,nr,r));
02730     ND_NuMuCC_Minus_Total.push_back(new TH2D(Form("ND_NuMuCC_Minus_Total_Run%i",ie+1),"",np,p,nr,r));
02731     ND_BNueCC_Plus_Total.push_back(new TH2D(Form("ND_BNueCC_Plus_Total_Run%i",ie+1),"",np,p,nr,r));
02732     ND_BNueCC_Minus_Total.push_back(new TH2D(Form("ND_BNueCC_Minus_Total_Run%i",ie+1),"",np,p,nr,r));
02733     ND_Total_Plus_Total.push_back(new TH2D(Form("ND_Total_Plus_Total_Run%i",ie+1),"",np,p,nr,r));
02734     ND_Total_Minus_Total.push_back(new TH2D(Form("ND_Total_Minus_Total_Run%i",ie+1),"",np,p,nr,r));
02735   }
02736   ND_Total_Plus_Total_AllRuns = new TH2D("ND_Total_Plus_Total_AllRuns","",np,p,nr,r);
02737   ND_Total_Minus_Total_AllRuns = new TH2D("ND_Total_Minus_Total_AllRuns","",np,p,nr,r);
02738   
02739   string systname;
02740   double pl=0,mn=0,temp=0;
02741   double bksum=0,cc=0,nc=0,bnue=0;
02742   double ccpl=0,ncpl=0,bpl=0;
02743   double ccmn=0,ncmn=0,bmn=0;
02744   double psum=0,nsum=0;
02745   double ccpsum=0,ccnsum=0;
02746   double ncpsum=0,ncnsum=0;
02747   double bpsum=0,bnsum=0;
02748   TH2D *ch1,*ch2,*ch3;
02749   
02750   std::map<string, vector<TH2D*> >::iterator systiter = ND_NC_Plus1Sigma.begin();
02751   std::map<string, vector<TH2D*> >::iterator last  = ND_NC_Plus1Sigma.end();
02752   
02753   for(ie=0;ie<Extrap.size();ie++)
02754   {
02755     for(ip=0;ip<np;ip++)
02756     {
02757       for(ir=0;ir<nr;ir++)
02758       {
02759         pl=0;
02760         mn=0;
02761         systiter = ND_NC_Plus1Sigma.begin();
02762         while(systiter!=last)
02763         {
02764           systname = systiter->first;
02765           
02766           temp = ND_NC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1);
02767           if(temp>0) pl+=temp*temp;
02768           else mn+=temp*temp;
02769           
02770           temp = ND_NC_Minus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1);
02771           if(temp>0) pl+=temp*temp;
02772           else mn+=temp*temp;
02773           
02774           systiter++;
02775         }
02776         ND_NC_Plus_Total[ie]->SetBinContent(ip+1,ir+1,sqrt(pl));
02777         ND_NC_Minus_Total[ie]->SetBinContent(ip+1,ir+1,-1.*sqrt(mn));
02778         
02779         pl=0;
02780         mn=0;
02781         systiter = ND_NC_Plus1Sigma.begin();
02782         while(systiter!=last)
02783         {
02784           systname = systiter->first;
02785           
02786           temp = ND_NuMuCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1);
02787           if(temp>0) pl+=temp*temp;
02788           else mn+=temp*temp;
02789           
02790           temp = ND_NuMuCC_Minus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1);
02791           if(temp>0) pl+=temp*temp;
02792           else mn+=temp*temp;
02793           
02794           systiter++;
02795         }
02796         ND_NuMuCC_Plus_Total[ie]->SetBinContent(ip+1,ir+1,sqrt(pl));
02797         ND_NuMuCC_Minus_Total[ie]->SetBinContent(ip+1,ir+1,-1.*sqrt(mn));
02798         
02799         pl=0;
02800         mn=0;
02801         systiter = ND_NC_Plus1Sigma.begin();
02802         while(systiter!=last)
02803         {
02804           systname = systiter->first;
02805           
02806           temp = ND_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1);
02807           if(temp>0) pl+=temp*temp;
02808           else mn+=temp*temp;
02809           
02810           temp = ND_BNueCC_Minus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1);
02811           if(temp>0) pl+=temp*temp;
02812           else mn+=temp*temp;
02813           
02814           systiter++;
02815         }
02816         ND_BNueCC_Plus_Total[ie]->SetBinContent(ip+1,ir+1,sqrt(pl));
02817         ND_BNueCC_Minus_Total[ie]->SetBinContent(ip+1,ir+1,-1.*sqrt(mn));
02818         
02819         pl=0;
02820         mn=0;
02821         systiter = ND_NC_Plus1Sigma.begin();
02822         while(systiter!=last)
02823         {
02824           systname = systiter->first;
02825           
02826           temp = ND_NC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->NDMC[Background::kNC]->GetBinContent(ip+1,ir+1) + ND_NuMuCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->NDMC[Background::kNuMuCC]->GetBinContent(ip+1,ir+1) + ND_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->NDMC[Background::kBNueCC]->GetBinContent(ip+1,ir+1);
02827           if(temp>0) pl+=temp*temp;
02828           else mn+=temp*temp;
02829           
02830           temp = ND_NC_Minus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->NDMC[Background::kNC]->GetBinContent(ip+1,ir+1) + ND_NuMuCC_Minus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->NDMC[Background::kNuMuCC]->GetBinContent(ip+1,ir+1) + ND_BNueCC_Minus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->NDMC[Background::kBNueCC]->GetBinContent(ip+1,ir+1);
02831           if(temp>0) pl+=temp*temp;
02832           else mn+=temp*temp;
02833           
02834           systiter++;
02835         }
02836         temp = Extrap[ie]->NDMC[Background::kNC]->GetBinContent(ip+1,ir+1) + Extrap[ie]->NDMC[Background::kNuMuCC]->GetBinContent(ip+1,ir+1) + Extrap[ie]->NDMC[Background::kBNueCC]->GetBinContent(ip+1,ir+1);
02837         
02838         ND_Total_Plus_Total[ie]->SetBinContent(ip+1,ir+1,sqrt(pl)/temp);
02839         ND_Total_Minus_Total[ie]->SetBinContent(ip+1,ir+1,-1.*sqrt(mn)/temp);
02840       }
02841     }
02842   }
02843   
02844   for(ip=0;ip<np;ip++)
02845   {
02846     for(ir=0;ir<nr;ir++)
02847     {
02848       pl=0;
02849       mn=0;
02850       systiter = ND_NC_Plus1Sigma.begin();
02851       while(systiter!=last)
02852       {
02853         systname = systiter->first;
02854         temp=0;
02855         for(ie=0;ie<Extrap.size();ie++)
02856         {
02857           temp += (ND_NC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->NDMC[Background::kNC]->GetBinContent(ip+1,ir+1) + ND_NuMuCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->NDMC[Background::kNuMuCC]->GetBinContent(ip+1,ir+1) + ND_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->NDMC[Background::kBNueCC]->GetBinContent(ip+1,ir+1));
02858         }
02859         if(temp>0) pl+=temp*temp;
02860         else mn+=temp*temp;
02861         
02862         temp=0;
02863         for(ie=0;ie<Extrap.size();ie++)
02864         { 
02865           temp += (ND_NC_Minus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->NDMC[Background::kNC]->GetBinContent(ip+1,ir+1) + ND_NuMuCC_Minus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->NDMC[Background::kNuMuCC]->GetBinContent(ip+1,ir+1) + ND_BNueCC_Minus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->NDMC[Background::kBNueCC]->GetBinContent(ip+1,ir+1));
02866         }
02867         if(temp>0) pl+=temp*temp;
02868         else mn+=temp*temp;
02869         
02870         systiter++;
02871       }
02872       
02873       temp=0;
02874       for(ie=0;ie<Extrap.size();ie++)
02875       {
02876         temp += (Extrap[ie]->NDMC[Background::kNC]->GetBinContent(ip+1,ir+1) + Extrap[ie]->NDMC[Background::kNuMuCC]->GetBinContent(ip+1,ir+1) + Extrap[ie]->NDMC[Background::kBNueCC]->GetBinContent(ip+1,ir+1));
02877       }
02878       ND_Total_Plus_Total_AllRuns->SetBinContent(ip+1,ir+1,sqrt(pl)/temp);
02879       ND_Total_Minus_Total_AllRuns->SetBinContent(ip+1,ir+1,-1.*sqrt(mn)/temp);
02880     }
02881   }
02882   
02883   TFile *f = new TFile(gSystem->ExpandPathName(nderrfilename.c_str()),"RECREATE");
02884   for(ie=0;ie<Extrap.size();ie++)
02885   {
02886     ND_NC_Plus_Total[ie]->Write();
02887     ND_NC_Minus_Total[ie]->Write();
02888     ND_NuMuCC_Plus_Total[ie]->Write();
02889     ND_NuMuCC_Minus_Total[ie]->Write();
02890     ND_BNueCC_Plus_Total[ie]->Write();
02891     ND_BNueCC_Minus_Total[ie]->Write();
02892     ND_Total_Plus_Total[ie]->Write();
02893     ND_Total_Minus_Total[ie]->Write();
02894   }
02895   ND_Total_Plus_Total_AllRuns->Write();
02896   ND_Total_Minus_Total_AllRuns->Write();
02897   
02898   f->Close();
02899   
02900   cout<<"ND Syst Table"<<endl;
02901   
02902   ofstream myfile("Tables/ndtable.tex");
02903   myfile.precision(2);
02904   myfile<<fixed;
02905   
02906   myfile<<"\\begin{table}[htp]"<<endl;
02907   myfile<<"\\caption{}"<<endl;
02908   myfile<<"\\centering"<<endl;
02909   myfile<<"\\begin{tabular}{c|cc|cc|cc|cc}"<<endl;
02910   myfile<<"Syst.&NC&&NuMuCC&&BNueCC&&Total&\\\\"<<endl;
02911   myfile<<"\\hline"<<endl;
02912   myfile<<"\\hline"<<endl;
02913   
02914   systiter = ND_NC_Plus1Sigma.begin();
02915   psum=0;
02916   nsum=0;
02917   ccpsum=0;
02918   ccnsum=0;
02919   ncpsum=0;
02920   ncnsum=0;
02921   bpsum=0;
02922   bnsum=0;
02923   while(systiter!=last)
02924   {
02925     systname = systiter->first;
02926     
02927     pl=0;
02928     mn=0;
02929     bpl=0;
02930     bmn=0;
02931     ncpl=0;
02932     ncmn=0;
02933     ccpl=0;
02934     ccmn=0;
02935     bksum=0;
02936     bnue=0;
02937     cc=0;
02938     nc=0;
02939     for(ie=0;ie<Extrap.size();ie++)
02940     {
02941       ch1 = (TH2D*)Extrap[ie]->NDMC[Background::kNC]->Clone();
02942       ch1->Multiply(ND_NC_Plus1Sigma[systname][ie]);
02943       ncpl+=ch1->Integral();
02944       ch2 = (TH2D*)Extrap[ie]->NDMC[Background::kNuMuCC]->Clone();
02945       ch2->Multiply(ND_NuMuCC_Plus1Sigma[systname][ie]);
02946       ccpl+=ch2->Integral();
02947       ch1->Add(ch2);
02948       ch3 = (TH2D*)Extrap[ie]->NDMC[Background::kBNueCC]->Clone();
02949       ch3->Multiply(ND_BNueCC_Plus1Sigma[systname][ie]);
02950       bpl+=ch3->Integral();
02951       ch1->Add(ch3);
02952       pl += ch1->Integral();
02953       
02954       ch1 = (TH2D*)Extrap[ie]->NDMC[Background::kNC]->Clone();
02955       ch1->Multiply(ND_NC_Minus1Sigma[systname][ie]);
02956       ncmn+=ch1->Integral();
02957       ch2 = (TH2D*)Extrap[ie]->NDMC[Background::kNuMuCC]->Clone();
02958       ch2->Multiply(ND_NuMuCC_Minus1Sigma[systname][ie]);
02959       ccmn+=ch2->Integral();
02960       ch1->Add(ch2);
02961       ch3 = (TH2D*)Extrap[ie]->NDMC[Background::kBNueCC]->Clone();
02962       ch3->Multiply(ND_BNueCC_Minus1Sigma[systname][ie]);
02963       bmn+=ch3->Integral();
02964       ch1->Add(ch3);
02965       mn += ch1->Integral();
02966       
02967       bksum+=(Extrap[ie]->NDMC[Background::kNuMuCC]->Integral() + Extrap[ie]->NDMC[Background::kNC]->Integral() + Extrap[ie]->NDMC[Background::kBNueCC]->Integral());
02968       cc+=Extrap[ie]->NDMC[Background::kNuMuCC]->Integral();
02969       nc+=Extrap[ie]->NDMC[Background::kNC]->Integral();
02970       bnue+=Extrap[ie]->NDMC[Background::kBNueCC]->Integral();
02971     }
02972     
02973     myfile<<systname<<"&";
02974     myfile<<100.*ncpl/nc<<"\\% &";
02975     if(TMath::Abs(ncpl+ncmn)<1e-20) myfile<<"-&";
02976     else myfile<<100.*ncmn/nc<<"\\% &";
02977     myfile<<100.*ccpl/cc<<"\\% &";
02978     if(TMath::Abs(ccpl+ccmn)<1e-20) myfile<<"-&";
02979     else myfile<<100.*ccmn/cc<<"\\% &";
02980     myfile<<100.*bpl/bnue<<"\\% &";
02981     if(TMath::Abs(bpl+bmn)<1e-20) myfile<<"-&";
02982     else myfile<<100.*bmn/bnue<<"\\% &";
02983     myfile<<100.*pl/bksum<<"\\% &";
02984     if(TMath::Abs(pl+mn)<1e-20) myfile<<"- \\\\"<<endl;
02985     else myfile<<100.*mn/bksum<<"\\% \\\\"<<endl;
02986     
02987     if((pl>0 && mn>0) || (pl<0 && mn<0))
02988     {
02989       if(TMath::Abs(pl)>TMath::Abs(mn))
02990       {
02991         mn = -1.*pl;
02992       }
02993       else
02994       {
02995         pl = -1.*mn;
02996       }
02997     }
02998     if((ccpl>0 && ccmn>0) || (ccpl<0 && ccmn<0))
02999     {
03000       if(TMath::Abs(ccpl)>TMath::Abs(ccmn))
03001       {
03002         ccmn = -1.*ccpl;
03003       }
03004       else
03005       {
03006         ccpl = -1.*ccmn;
03007       }
03008     }
03009     if((ncpl>0 && ncmn>0) || (ncpl<0 && ncmn<0))
03010     {
03011       if(TMath::Abs(ncpl)>TMath::Abs(ncmn))
03012       {
03013         ncmn = -1.*ncpl;
03014       }
03015       else
03016       {
03017         ncpl = -1.*ncmn;
03018       }
03019     }
03020     if((bpl>0 && bmn>0) || (bpl<0 && bmn<0))
03021     {
03022       if(TMath::Abs(bpl)>TMath::Abs(bmn))
03023       {
03024         bmn = -1.*bpl;
03025       }
03026       else
03027       {
03028         bpl = -1.*bmn;
03029       }
03030     }
03031     
03032     if(pl>0) psum+=(pl*pl);
03033     else nsum+=(pl*pl);
03034     if(mn>0) psum+=(mn*mn);
03035     else nsum+=(mn*mn);
03036     
03037     if(ccpl>0) ccpsum+=(ccpl*ccpl);
03038     else ccnsum+=(ccpl*ccpl);
03039     if(ccmn>0) ccpsum+=(ccmn*ccmn);
03040     else ccnsum+=(ccmn*ccmn);
03041     
03042     if(ncpl>0) ncpsum+=(ncpl*ncpl);
03043     else ncnsum+=(ncpl*ncpl);
03044     if(ncmn>0) ncpsum+=(ncmn*ncmn);
03045     else ncnsum+=(ncmn*ncmn);
03046     
03047     if(bpl>0) bpsum+=(bpl*bpl);
03048     else bnsum+=(bpl*bpl);
03049     if(bmn>0) bpsum+=(bmn*bmn);
03050     else bnsum+=(bmn*bmn);
03051     
03052     systiter++;
03053   }
03054   
03055   psum = sqrt(psum);
03056   nsum = -1.*sqrt(nsum);
03057   
03058   ccpsum = sqrt(ccpsum);
03059   ccnsum = -1.*sqrt(ccnsum);
03060   
03061   ncpsum = sqrt(ncpsum);
03062   ncnsum = -1.*sqrt(ncnsum);
03063   
03064   bpsum = sqrt(bpsum);
03065   bnsum = -1.*sqrt(bnsum);
03066   
03067   myfile<<"\\hline"<<endl;
03068   myfile<<"Total Extrap&"<<100.*ncpsum/nc<<"&"<<100.*ncnsum/nc<<"&"<<100.*ccpsum/cc<<"&"<<100.*ccnsum/cc<<"&"<<100.*bpsum/bnue<<"&"<<100.*bnsum/bnue<<"&"<<100.*psum/bksum<<"\\% &"<<100.*nsum/bksum<<"\\% \\\\"<<endl;
03069   myfile<<"\\end{tabular}"<<endl;
03070   myfile<<"\\end{table}"<<endl;
03071   myfile.close();
03072   
03073   return;
03074 }

void ErrorCalc::ReadCovarianceFiles ( int  n  )  [protected, virtual]

Reimplemented in ErrorCalc_Joint.

Definition at line 2045 of file ErrorCalc.cxx.

References ExtraCovariance, ExtraCovariance_File, ExtraCovariance_HistName, ExtraCovariance_SystName, Extrap, and gSystem().

Referenced by InitializeSys().

02046 {
02047   string systname = ExtraCovariance_SystName.at(n);
02048   string file = ExtraCovariance_File.at(n);
02049   string hist = ExtraCovariance_HistName.at(n);
02050 //   int flag = ExtraCovariance_Flag.at(n);
02051   
02052   if(gSystem->AccessPathName(gSystem->ExpandPathName(file.c_str())))
02053   {
02054     cout<<"Can't read "<<file<<endl;
02055     return;
02056   }
02057   
02058   TFile *f = new TFile(file.c_str(),"READ");
02059   TH2D *h = (TH2D*)f->Get(hist.c_str());
02060   
02061   int nbins = Extrap[0]->GetNPID()*Extrap[0]->GetNReco();
02062   
02063   if(h->GetNbinsX() != nbins)
02064   {
02065     cout<<"Error in ReadCovarianceFiles(): # bins from matrix in file doesn't match extrapolation bins!"<<endl;
02066     return;
02067   }
02068   
02069   ExtraCovariance[systname] = h;
02070   
02071   return;
02072 }

void ErrorCalc::ReadSpecialFiles ( int  n  )  [protected, virtual]

Reimplemented in ErrorCalc_Joint.

Definition at line 1868 of file ErrorCalc.cxx.

References Extrap, FD_BNueCC_Minus1Sigma, FD_BNueCC_Plus1Sigma, FD_NC_Minus1Sigma, FD_NC_Plus1Sigma, FD_NuMuCC_Minus1Sigma, FD_NuMuCC_Plus1Sigma, FN_BNueCC_Minus1Sigma, FN_BNueCC_Plus1Sigma, FN_NC_Minus1Sigma, FN_NC_Plus1Sigma, FN_NuMuCC_Minus1Sigma, FN_NuMuCC_Plus1Sigma, Form(), gSystem(), ND_BNueCC_Minus1Sigma, ND_BNueCC_Plus1Sigma, ND_NC_Minus1Sigma, ND_NC_Plus1Sigma, ND_NuMuCC_Minus1Sigma, ND_NuMuCC_Plus1Sigma, NueCC_MC_Minus1Sigma, NueCC_MC_Plus1Sigma, NuTauCC_MC_Minus1Sigma, NuTauCC_MC_Plus1Sigma, Pred_CC_Fid_Minus1Sigma, Pred_CC_Fid_Plus1Sigma, Extrapolate2D::qNuMuToNuTau, SpecialSyst_FileTag, SpecialSyst_Flag, SpecialSyst_HistName, SpecialSyst_SystName, SpecialSyst_Value, and SysFileDir.

Referenced by InitializeSys().

01869 {
01870   string systname = SpecialSyst_SystName.at(n);
01871   TH2D *h;
01872   unsigned int ie;
01873   int i,nr,np,nt;
01874   double *r,*p,*t;
01875   nr = Extrap[0]->GetNReco();
01876   np = Extrap[0]->GetNPID();
01877   nt = Extrap[0]->GetNTrue();
01878   r = new double[nr+1];
01879   p = new double[np+1];
01880   t = new double[nt+1];
01881   for(i=0;i<nr+1;i++)
01882   {
01883     r[i] = Extrap[0]->Pred_TotalBkgd->GetYaxis()->GetBinLowEdge(i+1);
01884   }
01885   for(i=0;i<np+1;i++)
01886   {
01887     p[i] = Extrap[0]->Pred_TotalBkgd->GetXaxis()->GetBinLowEdge(i+1);
01888   }
01889   for(i=0;i<nt+1;i++)
01890   {
01891     t[i] = Extrap[0]->Pred_3D[Extrapolate2D::qNuMuToNuTau]->GetZaxis()->GetBinLowEdge(i+1);
01892   }
01893   
01894   //create empty histograms for all shifts
01895   for(ie=0;ie<Extrap.size();ie++)
01896   {
01897     FN_NC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01898     FN_NC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01899     
01900     FN_NuMuCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NuMuCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01901     FN_NuMuCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NuMuCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01902     
01903     FN_BNueCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_BNueCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01904     FN_BNueCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_BNueCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01905     
01906     ND_NC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01907     ND_NC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01908     
01909     ND_NuMuCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NuMuCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01910     ND_NuMuCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NuMuCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01911     
01912     ND_BNueCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_BNueCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01913     ND_BNueCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_BNueCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01914     
01915     FD_NC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01916     FD_NC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01917     
01918     FD_NuMuCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NuMuCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01919     FD_NuMuCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NuMuCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01920     
01921     FD_BNueCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_BNueCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01922     FD_BNueCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_BNueCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01923     
01924     NueCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01925     NueCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01926     NuTauCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01927     NuTauCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01928     
01929     Pred_CC_Fid_Plus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Plus1",systname.c_str(),ie+1),"",nt,t));
01930     Pred_CC_Fid_Minus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Minus1",systname.c_str(),ie+1),"",nt,t));
01931   }
01932   
01933   int flag = SpecialSyst_Flag.at(n);
01934   
01935   if(SpecialSyst_Value.at(n)>0)//if value is set, use that number in all bins
01936   {
01937     h = (TH2D*)Extrap[0]->Pred_TotalBkgd->Clone("htemp");
01938     h->Reset();
01939     for(int ip=0;ip<Extrap[0]->GetNPID();ip++)
01940     {
01941       for(int ir=0;ir<Extrap[0]->GetNReco();ir++)
01942       {
01943         h->SetBinContent(ip+1,ir+1,SpecialSyst_Value.at(n));
01944       }
01945     }
01946     
01947     for(unsigned int ie=0;ie<Extrap.size();ie++)
01948     {
01949       if(flag==0 || flag==1)//for NC,NuMuCC,BNueCC
01950       {
01951         FN_NC_Plus1Sigma[systname][ie]->Add(h);
01952         FN_NC_Minus1Sigma[systname][ie]->Add(h,-1);
01953         
01954         FN_NuMuCC_Plus1Sigma[systname][ie]->Add(h);
01955         FN_NuMuCC_Minus1Sigma[systname][ie]->Add(h,-1);
01956         
01957         FN_BNueCC_Plus1Sigma[systname][ie]->Add(h);
01958         FN_BNueCC_Minus1Sigma[systname][ie]->Add(h,-1);
01959         
01960         FD_BNueCC_Plus1Sigma[systname][ie]->Add(h);
01961         FD_BNueCC_Minus1Sigma[systname][ie]->Add(h,-1);
01962       }
01963       if(flag==0 || flag==2)
01964       {
01965         NuTauCC_MC_Plus1Sigma[systname][ie]->Add(h);
01966         NuTauCC_MC_Minus1Sigma[systname][ie]->Add(h,-1);
01967       }
01968       if(flag==0 || flag==3)
01969       {
01970         NueCC_MC_Plus1Sigma[systname][ie]->Add(h);
01971         NueCC_MC_Minus1Sigma[systname][ie]->Add(h,-1);
01972       }
01973       if(flag==4)//for NC and NuMuCC (fake HOO)
01974       {
01975         FN_NC_Plus1Sigma[systname][ie]->Add(h);
01976         FN_NC_Minus1Sigma[systname][ie]->Add(h,-1);
01977         
01978         FN_NuMuCC_Plus1Sigma[systname][ie]->Add(h);
01979         FN_NuMuCC_Minus1Sigma[systname][ie]->Add(h,-1);
01980 
01981       }
01982 
01983 
01984     }
01985     
01986     return;
01987   }
01988   
01989   //if value is 0, read from a file
01990   string filetag,histname;
01991   filetag=SpecialSyst_FileTag.at(n);
01992   histname = SpecialSyst_HistName.at(n);
01993   
01994   TFile *f;
01995   
01996   for(unsigned int ie=0;ie<Extrap.size();ie++)
01997   {
01998     if(gSystem->AccessPathName(gSystem->ExpandPathName(Form("%s/%s_%i.root",SysFileDir.c_str(),filetag.c_str(),Extrap[ie]->GetRunPeriod()))))
01999     {
02000       cout<<"Warning: "<<Form("%s/%s_%i.root",SysFileDir.c_str(),filetag.c_str(),Extrap[ie]->GetRunPeriod())<<" doesn't exist."<<endl;
02001       continue;
02002     }
02003     
02004     f = new TFile(gSystem->ExpandPathName(Form("%s/%s_%i.root",SysFileDir.c_str(),filetag.c_str(),Extrap[ie]->GetRunPeriod())),"READ");
02005   
02006     h = (TH2D*)f->Get(histname.c_str());
02007     
02008     if(flag==0 || flag==1)//for NC,NuMuCC,BNueCC
02009     {
02010       FN_NC_Plus1Sigma[systname][ie]->Add(h);
02011       FN_NC_Minus1Sigma[systname][ie]->Add(h,-1);
02012         
02013       FN_NuMuCC_Plus1Sigma[systname][ie]->Add(h);
02014       FN_NuMuCC_Minus1Sigma[systname][ie]->Add(h,-1);
02015         
02016       FN_BNueCC_Plus1Sigma[systname][ie]->Add(h);
02017       FN_BNueCC_Minus1Sigma[systname][ie]->Add(h,-1);
02018         
02019       FD_BNueCC_Plus1Sigma[systname][ie]->Add(h);
02020       FD_BNueCC_Minus1Sigma[systname][ie]->Add(h,-1);
02021     }
02022     if(flag==0 || flag==2)
02023     {
02024       NuTauCC_MC_Plus1Sigma[systname][ie]->Add(h);
02025       NuTauCC_MC_Minus1Sigma[systname][ie]->Add(h,-1);
02026     }
02027     if(flag==0 || flag==3)
02028     {
02029       NueCC_MC_Plus1Sigma[systname][ie]->Add(h);
02030       NueCC_MC_Minus1Sigma[systname][ie]->Add(h,-1);
02031     }
02032     if(flag==4)//for NC and NuMuCC (fake HOO)
02033       {
02034         FN_NC_Plus1Sigma[systname][ie]->Add(h);
02035         FN_NC_Minus1Sigma[systname][ie]->Add(h,-1);
02036 
02037         FN_NuMuCC_Plus1Sigma[systname][ie]->Add(h);
02038         FN_NuMuCC_Minus1Sigma[systname][ie]->Add(h,-1);
02039       }
02040 
02041   }
02042   
02043   return;
02044 }

void ErrorCalc::ReadSysFiles_Appearance ( int  n  )  [protected]

Definition at line 895 of file ErrorCalc.cxx.

References Background::AsString(), Extrap, FNExtrap_FarMinusTag, FNExtrap_FarPlusTag, FNExtrap_FileTag, FNExtrap_Flag, FNExtrap_NearMinusTag, FNExtrap_NearPlusTag, FNExtrap_StdTag, FNExtrap_SystName, Form(), OscFit::GetRunPeriod(), gSystem(), it, Background::kNueCC, Background::kNuTauCC, MBH, NueCC_MC_Minus1Sigma, NueCC_MC_Plus1Sigma, NuTauCC_MC_Minus1Sigma, NuTauCC_MC_Plus1Sigma, Pred_CC_Fid_Minus1Sigma, Pred_CC_Fid_Plus1Sigma, Extrapolate2D::qNuMuToNue, Extrapolate2D::qNuMuToNuTau, MultiBinAnaHelper::Rebin2DHist(), and SysFileDir.

Referenced by InitializeSys().

00896 {
00897   int i;
00898   
00899   string plf,pln,mnf,mnn,std,systname,filetag;
00900   systname = FNExtrap_SystName.at(n);
00901   plf = FNExtrap_FarPlusTag.at(n);
00902   pln = FNExtrap_NearPlusTag.at(n);
00903   mnf = FNExtrap_FarMinusTag.at(n);
00904   mnn = FNExtrap_NearMinusTag.at(n);
00905   std = FNExtrap_StdTag.at(n);
00906   filetag=FNExtrap_FileTag.at(n);
00907   int flag = FNExtrap_Flag.at(n);
00908   
00909   vector<string> FidFiles;
00910   vector<string> CClikeFiles;
00911   vector<string> PreFiles;
00912   
00913   unsigned int ie;
00914   for(ie=0;ie<Extrap.size();ie++)
00915   {
00916     if(flag==1)//systs related to NuMuCC-like prediction
00917     {
00918       FidFiles.push_back(SysFileDir+"/SysFile2D_"+filetag+"_Fid_"+Form("%i",Extrap[ie]->GetRunPeriod())+".root");
00919       CClikeFiles.push_back(SysFileDir+"/SysFile2D_"+filetag+"_CCLike_"+Form("%i",Extrap[ie]->GetRunPeriod())+".root");
00920       if(gSystem->AccessPathName(gSystem->ExpandPathName(FidFiles[ie].c_str())))
00921       {
00922         cout<<"Warning: "<<FidFiles[ie]<<" doesn't exist."<<endl;
00923       }
00924       if(gSystem->AccessPathName(gSystem->ExpandPathName(CClikeFiles[ie].c_str())))
00925       {
00926         cout<<"Warning: "<<CClikeFiles[ie]<<" doesn't exist."<<endl;
00927       }
00928     }
00929     if(flag==2)//uncertainty in selected tau or signal FD MC
00930     {
00931       PreFiles.push_back(SysFileDir+"/SysFile2D_"+filetag+"_Pre_"+Form("%i",Extrap[ie]->GetRunPeriod())+".root");
00932       if(gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[ie].c_str())))
00933       {
00934         cout<<"Warning: "<<PreFiles[ie]<<" doesn't exist."<<endl;
00935       }
00936     }
00937   }
00938   
00939   unsigned int j;
00940   vector<Background::Background_t> app;
00941   app.push_back(Background::kNueCC);
00942   app.push_back(Background::kNuTauCC);
00943   
00944   string name;
00945   int np,nr,nt;
00946   double *r,*p,*t;
00947   TFile *fpre,*ffid,*fcc;
00948   double npot_pre,fpot_pre;
00949   double npot_fid,fpot_fid;
00950   double npot_cc,fpot_cc;
00951   TTree *treepre,*treefid,*treecc;
00952   double nPOTNear,nPOTFar;
00953   TH2D *prestd,*prepl,*premn;
00954   double sum,temp;
00955   int ir,it,ip;
00956   TH1D *ccstd,*ccpl,*ccmn;
00957   TH1D *data;
00958   TH2D *farfidccstd,*farcclikestd,*nearcclikestd;
00959   TH2D *farfidccpl,*farcclikepl,*nearcclikepl;
00960   TH2D *farfidccmn,*farcclikemn,*nearcclikemn;
00961   TH2D *h2;
00962   double sum_plus_tau,sum_minus_tau,sum_plus_nue,sum_minus_nue;
00963   string ccname = "CC";
00964   
00965   for(ie=0;ie<Extrap.size();ie++)
00966   {
00967 
00968     nr = Extrap[ie]->GetNReco();
00969     np = Extrap[ie]->GetNPID();
00970     nt = Extrap[ie]->GetNTrue();
00971     r = new double[nr+1];
00972     p = new double[np+1];
00973     t = new double[nt+1];
00974     for(i=0;i<nr+1;i++)
00975     {
00976       r[i] = Extrap[ie]->Pred_TotalBkgd->GetYaxis()->GetBinLowEdge(i+1);
00977     }
00978     for(i=0;i<np+1;i++)
00979     {
00980       p[i] = Extrap[ie]->Pred_TotalBkgd->GetXaxis()->GetBinLowEdge(i+1);
00981     }
00982     for(i=0;i<nt+1;i++)
00983     {
00984       t[i] = Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNuTau]->GetZaxis()->GetBinLowEdge(i+1);
00985     }
00986     nPOTFar=Extrap[ie]->GetFarPOT();
00987     nPOTNear=Extrap[ie]->GetNearPOT();
00988     if(flag==0)
00989     {
00990       //create empty histograms for overall signal/tau shifts
00991       NueCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00992       NueCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00993       NuTauCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00994       NuTauCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00995       
00996       //creating empty numuCC shifts
00997       Pred_CC_Fid_Plus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Plus1",systname.c_str(),ie+1),"",nt,t));
00998       Pred_CC_Fid_Minus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Minus1",systname.c_str(),ie+1),"",nt,t));
00999     }
01000     if(flag==1)
01001     {
01002       //create empty histograms for overall signal/tau shifts
01003       NueCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01004       NueCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01005       NuTauCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01006       NuTauCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01007       
01008       if(gSystem->AccessPathName(gSystem->ExpandPathName(CClikeFiles[ie].c_str())) || gSystem->AccessPathName(gSystem->ExpandPathName(FidFiles[ie].c_str())))//if CClike or Fid file does NOT exist
01009       {
01010         Pred_CC_Fid_Plus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Plus1",systname.c_str(),ie+1),"",nt,t));
01011         Pred_CC_Fid_Minus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Minus1",systname.c_str(),ie+1),"",nt,t));
01012         
01013         continue;
01014       }
01015       
01016       fcc = new TFile(gSystem->ExpandPathName(CClikeFiles[ie].c_str()),"READ");
01017       treecc = (TTree*)fcc->Get("paramtree");
01018       treecc->SetBranchAddress("nearPOT",&npot_cc);
01019       treecc->SetBranchAddress("farPOT",&fpot_cc);
01020       treecc->GetEntry(0);
01021       
01022       ffid = new TFile(gSystem->ExpandPathName(FidFiles[ie].c_str()),"READ");
01023       treefid = (TTree*)ffid->Get("paramtree");
01024       treefid->SetBranchAddress("nearPOT",&npot_fid);
01025       treefid->SetBranchAddress("farPOT",&fpot_fid);
01026       treefid->GetEntry(0);
01027       
01028       data = (TH1D*)Extrap[ie]->NDData_Reco_CClike->Clone(Form("data_%i",ie));
01029       
01030       name = "NuMuCC_" + std + "_Fid/FD_TrueVsReco";
01031       farfidccstd = (TH2D*)ffid->Get(name.c_str());
01032       farfidccstd->Scale(nPOTFar/fpot_fid);
01033       
01034       name = "NuMuCC_" + std + "_" + ccname + "/FD_TrueVsReco";
01035       farcclikestd = (TH2D*)fcc->Get(name.c_str());
01036       name = "NC_" + std + "_" + ccname + "/FD_TrueVsReco";
01037       h2 = (TH2D*)fcc->Get(name.c_str());
01038       farcclikestd->Add(h2);
01039       name = "BNueCC_" + std + "_" + ccname + "/FD_TrueVsReco";
01040       h2 = (TH2D*)fcc->Get(name.c_str());
01041       farcclikestd->Add(h2);
01042       farcclikestd->Scale(nPOTFar/fpot_cc);
01043       
01044       name = "NuMuCC_" + std + "_" + ccname + "/ND_TrueVsReco";
01045       nearcclikestd = (TH2D*)fcc->Get(name.c_str());
01046       name = "NC_" + std + "_" + ccname + "/ND_TrueVsReco";
01047       h2 = (TH2D*)fcc->Get(name.c_str());
01048       nearcclikestd->Add(h2);
01049       name = "BNueCC_" + std + "_" + ccname + "/ND_TrueVsReco";
01050       h2 = (TH2D*)fcc->Get(name.c_str());
01051       nearcclikestd->Add(h2);
01052       nearcclikestd->Scale(nPOTNear/npot_cc);
01053       
01054       ccstd = new TH1D("ccstd","",nt,t);
01055       for(it=0;it<nt;it++)
01056       {
01057         sum=0;
01058         for(ir=0;ir<farcclikestd->GetNbinsX();ir++)
01059         {
01060           if(nearcclikestd->Integral(ir+1,ir+1,1,nt)>0)
01061           {
01062             sum += (data->GetBinContent(ir+1)/nearcclikestd->Integral(ir+1,ir+1,1,nt))*farcclikestd->GetBinContent(ir+1,it+1);
01063           }
01064         }
01065         temp=0;
01066         if(farcclikestd->Integral(1,farcclikestd->GetNbinsX(),it+1,it+1)>0)
01067         {
01068           temp = farfidccstd->Integral(1,farfidccstd->GetNbinsX(),it+1,it+1)*(sum/farcclikestd->Integral(1,farcclikestd->GetNbinsX(),it+1,it+1));
01069         }
01070         ccstd->SetBinContent(it+1,temp);
01071       }
01072       
01073       name = "NuMuCC_" + plf + "_Fid/FD_TrueVsReco";
01074       if(plf==std)
01075       {
01076         farfidccpl = (TH2D*)farfidccstd->Clone(name.c_str());
01077       }
01078       else
01079       {
01080         farfidccpl = (TH2D*)ffid->Get(name.c_str());
01081         farfidccpl->Scale(nPOTFar/fpot_fid);
01082       }
01083       
01084       name = "NuMuCC_" + plf + "_" + ccname + "/FD_TrueVsReco";
01085       if(plf==std)
01086       {
01087         farcclikepl = (TH2D*)farcclikestd->Clone(name.c_str());
01088       }
01089       else
01090       {
01091         farcclikepl = (TH2D*)fcc->Get(name.c_str());
01092         name = "NC_" + plf + "_" + ccname + "/FD_TrueVsReco";
01093         h2 = (TH2D*)fcc->Get(name.c_str());
01094         farcclikepl->Add(h2);
01095         name = "BNueCC_" + plf + "_" + ccname + "/FD_TrueVsReco";
01096         h2 = (TH2D*)fcc->Get(name.c_str());
01097         farcclikepl->Add(h2);
01098         farcclikepl->Scale(nPOTFar/fpot_cc);
01099       }
01100       
01101       name = "NuMuCC_" + pln + "_" + ccname + "/ND_TrueVsReco";
01102       if(pln==std)
01103       {
01104         nearcclikepl = (TH2D*)nearcclikestd->Clone(name.c_str());
01105       }
01106       else
01107       {
01108         nearcclikepl = (TH2D*)fcc->Get(name.c_str());
01109         name = "NC_" + pln + "_" + ccname + "/ND_TrueVsReco";
01110         h2 = (TH2D*)fcc->Get(name.c_str());
01111         nearcclikepl->Add(h2);
01112         name = "BNueCC_" + pln + "_" + ccname + "/ND_TrueVsReco";
01113         h2 = (TH2D*)fcc->Get(name.c_str());
01114         nearcclikepl->Add(h2);
01115         nearcclikepl->Scale(nPOTNear/npot_cc);
01116       }
01117       
01118       ccpl = new TH1D("ccpl","",nt,t);
01119       for(it=0;it<nt;it++)
01120       {
01121         sum=0;
01122         for(ir=0;ir<farcclikepl->GetNbinsX();ir++)
01123         {
01124           if(nearcclikepl->Integral(ir+1,ir+1,1,nt)>0)
01125           {
01126             sum += (data->GetBinContent(ir+1)/nearcclikepl->Integral(ir+1,ir+1,1,nt))*farcclikepl->GetBinContent(ir+1,it+1);
01127           }
01128         }
01129         temp=0;
01130         if(farcclikepl->Integral(1,farcclikepl->GetNbinsX(),it+1,it+1)>0)
01131         {
01132           temp = farfidccpl->Integral(1,farfidccpl->GetNbinsX(),it+1,it+1)*(sum/farcclikepl->Integral(1,farcclikepl->GetNbinsX(),it+1,it+1));
01133         }
01134         ccpl->SetBinContent(it+1,temp);
01135       }
01136       
01137       name = "NuMuCC_" + mnf + "_Fid/FD_TrueVsReco";
01138       if(mnf==std)
01139       {
01140         farfidccmn = (TH2D*)farfidccstd->Clone(name.c_str());
01141       }
01142       else if(mnf==plf)
01143       {
01144         farfidccmn = (TH2D*)farfidccpl->Clone(name.c_str());
01145       }
01146       else
01147       {
01148         farfidccmn = (TH2D*)ffid->Get(name.c_str());
01149         farfidccmn->Scale(nPOTFar/fpot_fid);
01150       }
01151       
01152       name = "NuMuCC_" + mnf + "_" + ccname + "/FD_TrueVsReco";
01153       if(mnf==std)
01154       {
01155         farcclikemn = (TH2D*)farcclikestd->Clone(name.c_str());
01156       }
01157       else if(mnf==plf)
01158       {
01159         farcclikemn = (TH2D*)farcclikepl->Clone(name.c_str());
01160       }
01161       else
01162       {
01163         farcclikemn = (TH2D*)fcc->Get(name.c_str());
01164         name = "NC_" + mnf + "_" + ccname + "/FD_TrueVsReco";
01165         h2 = (TH2D*)fcc->Get(name.c_str());
01166         farcclikemn->Add(h2);
01167         name = "BNueCC_" + mnf + "_" + ccname + "/FD_TrueVsReco";
01168         h2 = (TH2D*)fcc->Get(name.c_str());
01169         farcclikemn->Add(h2);
01170         farcclikemn->Scale(nPOTFar/fpot_cc);
01171       }
01172       
01173       name = "NuMuCC_" + mnn + "_" + ccname + "/ND_TrueVsReco";
01174       if(mnn==std)
01175       {
01176         nearcclikemn = (TH2D*)nearcclikestd->Clone(name.c_str());
01177       }
01178       else if(mnn==pln)
01179       {
01180         nearcclikemn = (TH2D*)nearcclikepl->Clone(name.c_str());
01181       }
01182       else
01183       {
01184         nearcclikemn = (TH2D*)fcc->Get(name.c_str());
01185         name = "NC_" + mnn + "_" + ccname + "/ND_TrueVsReco";
01186         h2 = (TH2D*)fcc->Get(name.c_str());
01187         nearcclikemn->Add(h2);
01188         name = "BNueCC_" + mnn + "_" + ccname + "/ND_TrueVsReco";
01189         h2 = (TH2D*)fcc->Get(name.c_str());
01190         nearcclikemn->Add(h2);
01191         nearcclikemn->Scale(nPOTNear/npot_cc);
01192       }
01193       
01194       ccmn = new TH1D("ccmn","",nt,t);
01195       for(it=0;it<nt;it++)
01196       {
01197         sum=0;
01198         for(ir=0;ir<farcclikemn->GetNbinsX();ir++)
01199         {
01200           if(nearcclikemn->Integral(ir+1,ir+1,1,nt)>0)
01201           {
01202             sum += (data->GetBinContent(ir+1)/nearcclikemn->Integral(ir+1,ir+1,1,nt))*farcclikemn->GetBinContent(ir+1,it+1);
01203           }
01204         }
01205         temp=0;
01206         if(farcclikemn->Integral(1,farcclikemn->GetNbinsX(),it+1,it+1)>0)
01207         {
01208           temp = farfidccmn->Integral(1,farfidccmn->GetNbinsX(),it+1,it+1)*(sum/farcclikemn->Integral(1,farcclikemn->GetNbinsX(),it+1,it+1));
01209         }
01210         ccmn->SetBinContent(it+1,temp);
01211       }
01212       
01213       ccpl->Add(ccstd,-1.);
01214       ccpl->Divide(ccstd);
01215       ccmn->Add(ccstd,-1.);
01216       ccmn->Divide(ccstd);
01217       
01218       if(plf==mnf && pln==mnn)
01219       {
01220         ccmn->Scale(-1.);
01221       }
01222       
01223       Pred_CC_Fid_Plus1Sigma[systname].push_back((TH1D*)ccpl->Clone(Form("%s_%i_Pred_CC_Fid_Plus1",systname.c_str(),ie+1)));
01224       Pred_CC_Fid_Minus1Sigma[systname].push_back((TH1D*)ccmn->Clone(Form("%s_%i_Pred_CC_Fid_Minus1",systname.c_str(),ie+1)));
01225       
01226       for(ip=0;ip<np;ip++)
01227       {
01228         for(ir=0;ir<nr;ir++)
01229         {
01230           sum_plus_tau=0;
01231           sum_minus_tau=0;
01232           sum_plus_nue=0;
01233           sum_minus_nue=0;
01234           for(it=0;it<nt;it++)
01235           {
01236             sum_plus_tau += (Pred_CC_Fid_Plus1Sigma[systname][ie]->GetBinContent(it+1)*Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNuTau]->GetBinContent(ip+1,ir+1,it+1));
01237             
01238             sum_minus_tau += (Pred_CC_Fid_Minus1Sigma[systname][ie]->GetBinContent(it+1)*Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNuTau]->GetBinContent(ip+1,ir+1,it+1));
01239             
01240             sum_plus_nue += (Pred_CC_Fid_Plus1Sigma[systname][ie]->GetBinContent(it+1)*Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNue]->GetBinContent(ip+1,ir+1,it+1));
01241             
01242             sum_minus_nue += (Pred_CC_Fid_Minus1Sigma[systname][ie]->GetBinContent(it+1)*Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNue]->GetBinContent(ip+1,ir+1,it+1));
01243           }
01244           //make fractional errors again
01245           sum_plus_tau = sum_plus_tau/Extrap[ie]->Pred[Background::kNuTauCC]->GetBinContent(ip+1,ir+1);
01246           sum_minus_tau = sum_minus_tau/Extrap[ie]->Pred[Background::kNuTauCC]->GetBinContent(ip+1,ir+1);
01247           sum_plus_nue = sum_plus_nue/Extrap[ie]->Pred[Background::kNueCC]->GetBinContent(ip+1,ir+1);
01248           sum_minus_nue = sum_minus_nue/Extrap[ie]->Pred[Background::kNueCC]->GetBinContent(ip+1,ir+1);
01249           
01250           NuTauCC_MC_Plus1Sigma[systname][ie]->SetBinContent(ip+1,ir+1,sum_plus_tau);
01251           NuTauCC_MC_Minus1Sigma[systname][ie]->SetBinContent(ip+1,ir+1,sum_minus_tau);
01252           NueCC_MC_Plus1Sigma[systname][ie]->SetBinContent(ip+1,ir+1,sum_plus_nue);
01253           NueCC_MC_Minus1Sigma[systname][ie]->SetBinContent(ip+1,ir+1,sum_minus_nue);
01254         }
01255       }
01256     }
01257     
01258     if(flag==2)
01259     {
01260       //creating empty numuCC shifts
01261       Pred_CC_Fid_Plus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Plus1",systname.c_str(),ie+1),"",nt,t));
01262       Pred_CC_Fid_Minus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Minus1",systname.c_str(),ie+1),"",nt,t));
01263       
01264       if(gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[ie].c_str())))//if Pre file does NOT exist
01265       {
01266         NueCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01267         NueCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01268         NuTauCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01269         NuTauCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01270         
01271         continue;
01272       }
01273       
01274       fpre = new TFile(gSystem->ExpandPathName(PreFiles[ie].c_str()),"READ");
01275       treepre = (TTree*)fpre->Get("paramtree");
01276       treepre->SetBranchAddress("nearPOT",&npot_pre);
01277       treepre->SetBranchAddress("farPOT",&fpot_pre);
01278       treepre->GetEntry(0);
01279       
01280       for(j=0;j<app.size();j++)
01281       {
01282         name = string(Background::AsString(app[j])) + "_" + std + "_Presel/FD_RecoVs" + Extrap[0]->GetPID();
01283         prestd = (TH2D*)fpre->Get(name.c_str());
01284         prestd->Scale(nPOTFar/fpot_pre);
01285         MBH.Rebin2DHist(prestd,np,p,nr,r);
01286       
01287         name = string(Background::AsString(app[j])) + "_" + plf + "_Presel/FD_RecoVs" + Extrap[0]->GetPID();
01288         if(plf==std)
01289         {
01290           prepl = (TH2D*)prestd->Clone(name.c_str());
01291         }
01292         else
01293         {
01294           prepl = (TH2D*)fpre->Get(name.c_str());
01295           prepl->Scale(nPOTFar/fpot_pre);
01296           MBH.Rebin2DHist(prepl,np,p,nr,r);
01297         }
01298       
01299         name = string(Background::AsString(app[j])) + "_" + mnf + "_Presel/FD_RecoVs" + Extrap[0]->GetPID();
01300         if(mnf==std)
01301         {
01302           premn = (TH2D*)prestd->Clone(name.c_str());
01303         }
01304         else if(mnf==plf)
01305         {
01306           premn = (TH2D*)prepl->Clone(name.c_str());
01307         }
01308         else
01309         {
01310           premn = (TH2D*)fpre->Get(name.c_str());
01311           premn->Scale(nPOTFar/fpot_pre);
01312           MBH.Rebin2DHist(premn,np,p,nr,r);
01313         }
01314         
01315         prepl->Add(prestd,-1.);
01316         prepl->Divide(prestd);
01317         premn->Add(prestd,-1.);
01318         premn->Divide(prestd);
01319         
01320         if(plf==mnf && pln==mnn)
01321         {
01322           premn->Scale(-1.);
01323         }
01324         
01325         if(app[j]==Background::kNueCC)
01326         {
01327           NueCC_MC_Plus1Sigma[systname].push_back((TH2D*)prepl->Clone(Form("%s_%i_NueCC_MC_Plus1",systname.c_str(),ie+1)));
01328           NueCC_MC_Minus1Sigma[systname].push_back((TH2D*)premn->Clone(Form("%s_%i_NueCC_MC_Minus1",systname.c_str(),ie+1)));
01329         }
01330         else if(app[j]==Background::kNuTauCC)
01331         {
01332           NuTauCC_MC_Plus1Sigma[systname].push_back((TH2D*)prepl->Clone(Form("%s_%i_NuTauCC_MC_Plus1",systname.c_str(),ie+1)));
01333           NuTauCC_MC_Minus1Sigma[systname].push_back((TH2D*)premn->Clone(Form("%s_%i_NuTauCC_MC_Minus1",systname.c_str(),ie+1)));
01334         }
01335       }
01336     }
01337   }
01338   
01339   return;
01340 }

void ErrorCalc::ReadSysFiles_Appearance_AllRuns ( int  n  )  [protected, virtual]

Reimplemented in ErrorCalc_Joint.

Definition at line 1341 of file ErrorCalc.cxx.

References Background::AsString(), Extrap, FNExtrap_FarMinusTag, FNExtrap_FarPlusTag, FNExtrap_FileTag, FNExtrap_Flag, FNExtrap_NearMinusTag, FNExtrap_NearPlusTag, FNExtrap_StdTag, FNExtrap_SystName, Form(), OscFit::GetRunPeriod(), gSystem(), it, Background::kNueCC, Background::kNuTauCC, MBH, NueCC_MC_Minus1Sigma, NueCC_MC_Plus1Sigma, NuTauCC_MC_Minus1Sigma, NuTauCC_MC_Plus1Sigma, Pred_CC_Fid_Minus1Sigma, Pred_CC_Fid_Plus1Sigma, Extrapolate2D::qNuMuToNue, Extrapolate2D::qNuMuToNuTau, MultiBinAnaHelper::Rebin2DHist(), and SysFileDir.

Referenced by InitializeSys().

01342 {
01343   int i;
01344   
01345   string plf,pln,mnf,mnn,std,systname,filetag;
01346   systname = FNExtrap_SystName.at(n);
01347   plf = FNExtrap_FarPlusTag.at(n);
01348   pln = FNExtrap_NearPlusTag.at(n);
01349   mnf = FNExtrap_FarMinusTag.at(n);
01350   mnn = FNExtrap_NearMinusTag.at(n);
01351   std = FNExtrap_StdTag.at(n);
01352   filetag=FNExtrap_FileTag.at(n);
01353   int flag = FNExtrap_Flag.at(n);
01354   
01355   vector<string> FidFiles;
01356   vector<string> CClikeFiles;
01357   vector<string> PreFiles;
01358   
01359   unsigned int ie;
01360   for(ie=0;ie<Extrap.size();ie++)
01361   {
01362     if(flag==1)//systs related to NuMuCC-like prediction
01363     {
01364       FidFiles.push_back(SysFileDir+"/SysFile2D_"+filetag+"_Fid_"+Form("%i",Extrap[ie]->GetRunPeriod())+".root");
01365       CClikeFiles.push_back(SysFileDir+"/SysFile2D_"+filetag+"_CCLike_"+Form("%i",Extrap[ie]->GetRunPeriod())+".root");
01366       if(gSystem->AccessPathName(gSystem->ExpandPathName(FidFiles[ie].c_str())))
01367       {
01368         cout<<"Warning: "<<FidFiles[ie]<<" doesn't exist."<<endl;
01369       }
01370       if(gSystem->AccessPathName(gSystem->ExpandPathName(CClikeFiles[ie].c_str())))
01371       {
01372         cout<<"Warning: "<<CClikeFiles[ie]<<" doesn't exist."<<endl;
01373       }
01374     }
01375     else if(flag==2)//uncertainty in selected tau or signal FD MC
01376     {
01377       PreFiles.push_back(SysFileDir+"/SysFile2D_"+filetag+"_Pre_"+Form("%i",Extrap[ie]->GetRunPeriod())+".root");
01378       if(gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[ie].c_str())))
01379       {
01380         cout<<"Warning: "<<PreFiles[ie]<<" doesn't exist."<<endl;
01381       }
01382     }
01383   }
01384   
01385   unsigned int j;
01386   vector<Background::Background_t> app;
01387   app.push_back(Background::kNueCC);
01388   app.push_back(Background::kNuTauCC);
01389   
01390   string name;
01391   int np,nr,nt;
01392   double *r,*p,*t;
01393   TFile *fpre,*ffid,*fcc;
01394   double npot_pre,fpot_pre;
01395   double npot_fid,fpot_fid;
01396   double npot_cc,fpot_cc;
01397   TTree *treepre,*treefid,*treecc;
01398   double nPOTNear,nPOTFar;
01399   TH2D *prestd=0,*prepl=0,*premn=0;
01400   TH2D *totalprestd=0,*totalprepl=0,*totalpremn=0;
01401   double sum,temp;
01402   int ir,it,ip;
01403   TH1D *ccstd=0,*ccpl=0,*ccmn=0;
01404   TH1D *data=0;
01405   TH2D *farfidccstd=0,*farcclikestd=0,*nearcclikestd=0;
01406   TH2D *farfidccpl=0,*farcclikepl=0,*nearcclikepl=0;
01407   TH2D *farfidccmn=0,*farcclikemn=0,*nearcclikemn=0;
01408   TH1D *totaldata=0;
01409   TH2D *totalfarfidccstd=0,*totalfarcclikestd=0,*totalnearcclikestd=0;
01410   TH2D *totalfarfidccpl=0,*totalfarcclikepl=0,*totalnearcclikepl=0;
01411   TH2D *totalfarfidccmn=0,*totalfarcclikemn=0,*totalnearcclikemn=0;
01412   TH2D *h2;
01413   int nfiles;
01414   double sum_plus_tau,sum_minus_tau,sum_plus_nue,sum_minus_nue;
01415   string ccname = "CC";
01416   
01417   nr = Extrap[0]->GetNReco();
01418   np = Extrap[0]->GetNPID();
01419   nt = Extrap[0]->GetNTrue();
01420   r = new double[nr+1];
01421   p = new double[np+1];
01422   t = new double[nt+1];
01423   for(i=0;i<nr+1;i++)
01424   {
01425     r[i] = Extrap[0]->Pred_TotalBkgd->GetYaxis()->GetBinLowEdge(i+1);
01426   }
01427   for(i=0;i<np+1;i++)
01428   {
01429     p[i] = Extrap[0]->Pred_TotalBkgd->GetXaxis()->GetBinLowEdge(i+1);
01430   }
01431   for(i=0;i<nt+1;i++)
01432   {
01433     t[i] = Extrap[0]->Pred_3D[Extrapolate2D::qNuMuToNuTau]->GetZaxis()->GetBinLowEdge(i+1);
01434   }
01435   
01436   if(flag==0)
01437   {
01438     for(ie=0;ie<Extrap.size();ie++)
01439     {
01440        //create empty histograms for overall signal/tau shifts
01441       NueCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01442       NueCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01443       NuTauCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01444       NuTauCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01445       
01446        //creating empty numuCC shifts
01447       Pred_CC_Fid_Plus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Plus1",systname.c_str(),ie+1),"",nt,t));
01448       Pred_CC_Fid_Minus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Minus1",systname.c_str(),ie+1),"",nt,t));
01449     }
01450   }
01451   
01452   if(flag==1)
01453   {
01454     //create empty histograms for overall signal/tau shifts
01455     for(ie=0;ie<Extrap.size();ie++)
01456     {
01457       NueCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01458       NueCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01459       NuTauCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01460       NuTauCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01461     }
01462     
01463     if(gSystem->AccessPathName(gSystem->ExpandPathName(CClikeFiles[0].c_str())) && gSystem->AccessPathName(gSystem->ExpandPathName(CClikeFiles[1].c_str())) && gSystem->AccessPathName(gSystem->ExpandPathName(CClikeFiles[2].c_str())) && gSystem->AccessPathName(gSystem->ExpandPathName(FidFiles[0].c_str())) && gSystem->AccessPathName(gSystem->ExpandPathName(FidFiles[1].c_str())) && gSystem->AccessPathName(gSystem->ExpandPathName(FidFiles[2].c_str())))//if none of these files exist
01464     {
01465       for(ie=0;ie<Extrap.size();ie++)
01466       {
01467         Pred_CC_Fid_Plus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Plus1",systname.c_str(),ie+1),"",nt,t));
01468         Pred_CC_Fid_Minus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Minus1",systname.c_str(),ie+1),"",nt,t));
01469       }
01470       
01471       return;
01472     }
01473     
01474     nfiles=0;
01475     for(ie=0;ie<3;ie++)
01476     {
01477       if(gSystem->AccessPathName(gSystem->ExpandPathName(CClikeFiles[ie].c_str())) && gSystem->AccessPathName(gSystem->ExpandPathName(FidFiles[ie].c_str())))//if files don't exist
01478       { 
01479         continue;
01480       }
01481       
01482       fcc = new TFile(gSystem->ExpandPathName(CClikeFiles[ie].c_str()),"READ");
01483       treecc = (TTree*)fcc->Get("paramtree");
01484       treecc->SetBranchAddress("nearPOT",&npot_cc);
01485       treecc->SetBranchAddress("farPOT",&fpot_cc);
01486       treecc->GetEntry(0);
01487       
01488       ffid = new TFile(gSystem->ExpandPathName(FidFiles[ie].c_str()),"READ");
01489       treefid = (TTree*)ffid->Get("paramtree");
01490       treefid->SetBranchAddress("nearPOT",&npot_fid);
01491       treefid->SetBranchAddress("farPOT",&fpot_fid);
01492       treefid->GetEntry(0);
01493       
01494       data = (TH1D*)Extrap[ie]->NDData_Reco_CClike->Clone(Form("data_%i",ie));
01495       //should scale data!
01496       if(nfiles==0) totaldata = (TH1D*)data->Clone();
01497       else totaldata->Add(data);
01498       
01499       nPOTFar=Extrap[ie]->GetFarPOT();
01500       nPOTNear=Extrap[ie]->GetNearPOT();
01501       
01502       name = "NuMuCC_" + std + "_Fid/FD_TrueVsReco";
01503       farfidccstd = (TH2D*)ffid->Get(name.c_str());
01504       farfidccstd->Scale(nPOTFar/fpot_fid);
01505       if(nfiles==0) totalfarfidccstd = (TH2D*)farfidccstd->Clone();
01506       else totalfarfidccstd->Add(farfidccstd);
01507       
01508       name = "NuMuCC_" + std + "_" + ccname + "/FD_TrueVsReco";
01509       farcclikestd = (TH2D*)fcc->Get(name.c_str());
01510       name = "NC_" + std + "_" + ccname + "/FD_TrueVsReco";
01511       h2 = (TH2D*)fcc->Get(name.c_str());
01512       farcclikestd->Add(h2);
01513       name = "BNueCC_" + std + "_" + ccname + "/FD_TrueVsReco";
01514       h2 = (TH2D*)fcc->Get(name.c_str());
01515       farcclikestd->Add(h2);
01516       farcclikestd->Scale(nPOTFar/fpot_cc);
01517       if(nfiles==0) totalfarcclikestd = (TH2D*)farcclikestd->Clone();
01518       else totalfarcclikestd->Add(farcclikestd);
01519       
01520       name = "NuMuCC_" + std + "_" + ccname + "/ND_TrueVsReco";
01521       nearcclikestd = (TH2D*)fcc->Get(name.c_str());
01522       name = "NC_" + std + "_" + ccname + "/ND_TrueVsReco";
01523       h2 = (TH2D*)fcc->Get(name.c_str());
01524       nearcclikestd->Add(h2);
01525       name = "BNueCC_" + std + "_" + ccname + "/ND_TrueVsReco";
01526       h2 = (TH2D*)fcc->Get(name.c_str());
01527       nearcclikestd->Add(h2);
01528       nearcclikestd->Scale(nPOTNear/npot_cc);
01529       if(nfiles==0) totalnearcclikestd = (TH2D*)nearcclikestd->Clone();
01530       else totalnearcclikestd->Add(nearcclikestd);
01531       
01532       name = "NuMuCC_" + plf + "_Fid/FD_TrueVsReco";
01533       if(plf==std)
01534       {
01535         farfidccpl = (TH2D*)farfidccstd->Clone(name.c_str());
01536       }
01537       else
01538       {
01539         farfidccpl = (TH2D*)ffid->Get(name.c_str());
01540         farfidccpl->Scale(nPOTFar/fpot_fid);
01541       }
01542       if(nfiles==0) totalfarfidccpl = (TH2D*)farfidccpl->Clone();
01543       else totalfarfidccpl->Add(farfidccpl);
01544       
01545       name = "NuMuCC_" + plf + "_" + ccname + "/FD_TrueVsReco";
01546       if(plf==std)
01547       {
01548         farcclikepl = (TH2D*)farcclikestd->Clone(name.c_str());
01549       }
01550       else
01551       {
01552         farcclikepl = (TH2D*)fcc->Get(name.c_str());
01553         name = "NC_" + plf + "_" + ccname + "/FD_TrueVsReco";
01554         h2 = (TH2D*)fcc->Get(name.c_str());
01555         farcclikepl->Add(h2);
01556         name = "BNueCC_" + plf + "_" + ccname + "/FD_TrueVsReco";
01557         h2 = (TH2D*)fcc->Get(name.c_str());
01558         farcclikepl->Add(h2);
01559         farcclikepl->Scale(nPOTFar/fpot_cc);
01560       }
01561       if(nfiles==0) totalfarcclikepl = (TH2D*)farcclikepl->Clone();
01562       else totalfarcclikepl->Add(farcclikepl);
01563       
01564       name = "NuMuCC_" + pln + "_" + ccname + "/ND_TrueVsReco";
01565       if(pln==std)
01566       {
01567         nearcclikepl = (TH2D*)nearcclikestd->Clone(name.c_str());
01568       }
01569       else
01570       {
01571         nearcclikepl = (TH2D*)fcc->Get(name.c_str());
01572         name = "NC_" + pln + "_" + ccname + "/ND_TrueVsReco";
01573         h2 = (TH2D*)fcc->Get(name.c_str());
01574         nearcclikepl->Add(h2);
01575         name = "BNueCC_" + pln + "_" + ccname + "/ND_TrueVsReco";
01576         h2 = (TH2D*)fcc->Get(name.c_str());
01577         nearcclikepl->Add(h2);
01578         nearcclikepl->Scale(nPOTNear/npot_cc);
01579       }
01580       if(nfiles==0) totalnearcclikepl = (TH2D*)nearcclikepl->Clone();
01581       else totalnearcclikepl->Add(nearcclikepl);
01582       
01583       name = "NuMuCC_" + mnf + "_Fid/FD_TrueVsReco";
01584       if(mnf==std)
01585       {
01586         farfidccmn = (TH2D*)farfidccstd->Clone(name.c_str());
01587       }
01588       else if(mnf==plf)
01589       {
01590         farfidccmn = (TH2D*)farfidccpl->Clone(name.c_str());
01591       }
01592       else
01593       {
01594         farfidccmn = (TH2D*)ffid->Get(name.c_str());
01595         farfidccmn->Scale(nPOTFar/fpot_fid);
01596       }
01597       if(nfiles==0) totalfarfidccmn = (TH2D*)farfidccmn->Clone();
01598       else totalfarfidccmn->Add(farfidccmn);
01599       
01600       name = "NuMuCC_" + mnf + "_" + ccname + "/FD_TrueVsReco";
01601       if(mnf==std)
01602       {
01603         farcclikemn = (TH2D*)farcclikestd->Clone(name.c_str());
01604       }
01605       else if(mnf==plf)
01606       {
01607         farcclikemn = (TH2D*)farcclikepl->Clone(name.c_str());
01608       }
01609       else
01610       {
01611         farcclikemn = (TH2D*)fcc->Get(name.c_str());
01612         name = "NC_" + mnf + "_" + ccname + "/FD_TrueVsReco";
01613         h2 = (TH2D*)fcc->Get(name.c_str());
01614         farcclikemn->Add(h2);
01615         name = "BNueCC_" + mnf + "_" + ccname + "/FD_TrueVsReco";
01616         h2 = (TH2D*)fcc->Get(name.c_str());
01617         farcclikemn->Add(h2);
01618         farcclikemn->Scale(nPOTFar/fpot_cc);
01619       }
01620       if(nfiles==0) totalfarcclikemn = (TH2D*)farcclikemn->Clone();
01621       else totalfarcclikemn->Add(farcclikemn);
01622       
01623       name = "NuMuCC_" + mnn + "_" + ccname + "/ND_TrueVsReco";
01624       if(mnn==std)
01625       {
01626         nearcclikemn = (TH2D*)nearcclikestd->Clone(name.c_str());
01627       }
01628       else if(mnn==pln)
01629       {
01630         nearcclikemn = (TH2D*)nearcclikepl->Clone(name.c_str());
01631       }
01632       else
01633       {
01634         nearcclikemn = (TH2D*)fcc->Get(name.c_str());
01635         name = "NC_" + mnn + "_" + ccname + "/ND_TrueVsReco";
01636         h2 = (TH2D*)fcc->Get(name.c_str());
01637         nearcclikemn->Add(h2);
01638         name = "BNueCC_" + mnn + "_" + ccname + "/ND_TrueVsReco";
01639         h2 = (TH2D*)fcc->Get(name.c_str());
01640         nearcclikemn->Add(h2);
01641         nearcclikemn->Scale(nPOTNear/npot_cc);
01642       }
01643       if(nfiles==0) totalnearcclikemn = (TH2D*)nearcclikemn->Clone();
01644       else totalnearcclikemn->Add(nearcclikemn);
01645       
01646       nfiles++;
01647     }
01648     
01649     ccstd = new TH1D("ccstd","",nt,t);
01650     for(it=0;it<nt;it++)
01651     {
01652       sum=0;
01653       for(ir=0;ir<totalfarcclikestd->GetNbinsX();ir++)
01654       {
01655         if(totalnearcclikestd->Integral(ir+1,ir+1,1,nt)>0)
01656         {
01657           sum += (totaldata->GetBinContent(ir+1)/totalnearcclikestd->Integral(ir+1,ir+1,1,nt))*totalfarcclikestd->GetBinContent(ir+1,it+1);
01658         }
01659       }
01660       temp=0;
01661       if(totalfarcclikestd->Integral(1,totalfarcclikestd->GetNbinsX(),it+1,it+1)>0)
01662       {
01663         temp = totalfarfidccstd->Integral(1,totalfarfidccstd->GetNbinsX(),it+1,it+1)*(sum/totalfarcclikestd->Integral(1,totalfarcclikestd->GetNbinsX(),it+1,it+1));
01664       }
01665       ccstd->SetBinContent(it+1,temp);
01666     }
01667     
01668     ccpl = new TH1D("ccpl","",nt,t);
01669     for(it=0;it<nt;it++)
01670     {
01671       sum=0;
01672       for(ir=0;ir<totalfarcclikepl->GetNbinsX();ir++)
01673       {
01674         if(totalnearcclikepl->Integral(ir+1,ir+1,1,nt)>0)
01675         {
01676           sum += (totaldata->GetBinContent(ir+1)/totalnearcclikepl->Integral(ir+1,ir+1,1,nt))*totalfarcclikepl->GetBinContent(ir+1,it+1);
01677         }
01678       }
01679       temp=0;
01680       if(totalfarcclikepl->Integral(1,totalfarcclikepl->GetNbinsX(),it+1,it+1)>0)
01681       {
01682         temp = totalfarfidccpl->Integral(1,totalfarfidccpl->GetNbinsX(),it+1,it+1)*(sum/totalfarcclikepl->Integral(1,totalfarcclikepl->GetNbinsX(),it+1,it+1));
01683       }
01684       ccpl->SetBinContent(it+1,temp);
01685     }
01686     
01687     ccmn = new TH1D("ccmn","",nt,t);
01688     for(it=0;it<nt;it++)
01689     {
01690       sum=0;
01691       for(ir=0;ir<totalfarcclikemn->GetNbinsX();ir++)
01692       {
01693         if(totalnearcclikemn->Integral(ir+1,ir+1,1,nt)>0)
01694         {
01695           sum += (totaldata->GetBinContent(ir+1)/totalnearcclikemn->Integral(ir+1,ir+1,1,nt))*totalfarcclikemn->GetBinContent(ir+1,it+1);
01696         }
01697       }
01698       temp=0;
01699       if(totalfarcclikemn->Integral(1,totalfarcclikemn->GetNbinsX(),it+1,it+1)>0)
01700       {
01701         temp = totalfarfidccmn->Integral(1,totalfarfidccmn->GetNbinsX(),it+1,it+1)*(sum/totalfarcclikemn->Integral(1,totalfarcclikemn->GetNbinsX(),it+1,it+1));
01702       }
01703       ccmn->SetBinContent(it+1,temp);
01704     }
01705     
01706     ccpl->Add(ccstd,-1.);
01707     ccpl->Divide(ccstd);
01708     ccmn->Add(ccstd,-1.);
01709     ccmn->Divide(ccstd);
01710     
01711     if(plf==mnf && pln==mnn)
01712     {
01713       ccmn->Scale(-1.);
01714     }
01715     
01716     for(ie=0;ie<Extrap.size();ie++)
01717     {
01718       Pred_CC_Fid_Plus1Sigma[systname].push_back((TH1D*)ccpl->Clone(Form("%s_%i_Pred_CC_Fid_Plus1",systname.c_str(),ie+1)));
01719       Pred_CC_Fid_Minus1Sigma[systname].push_back((TH1D*)ccmn->Clone(Form("%s_%i_Pred_CC_Fid_Minus1",systname.c_str(),ie+1)));
01720     }
01721     
01722     for(ie=0;ie<Extrap.size();ie++)
01723     {
01724       for(ip=0;ip<np;ip++)
01725       {
01726         for(ir=0;ir<nr;ir++)
01727         {
01728           sum_plus_tau=0;
01729           sum_minus_tau=0;
01730           sum_plus_nue=0;
01731           sum_minus_nue=0;
01732           for(it=0;it<nt;it++)
01733           {
01734             sum_plus_tau += (Pred_CC_Fid_Plus1Sigma[systname][ie]->GetBinContent(it+1)*Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNuTau]->GetBinContent(ip+1,ir+1,it+1));
01735             
01736             sum_minus_tau += (Pred_CC_Fid_Minus1Sigma[systname][ie]->GetBinContent(it+1)*Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNuTau]->GetBinContent(ip+1,ir+1,it+1));
01737             
01738             sum_plus_nue += (Pred_CC_Fid_Plus1Sigma[systname][ie]->GetBinContent(it+1)*Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNue]->GetBinContent(ip+1,ir+1,it+1));
01739             
01740             sum_minus_nue += (Pred_CC_Fid_Minus1Sigma[systname][ie]->GetBinContent(it+1)*Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNue]->GetBinContent(ip+1,ir+1,it+1));
01741           }
01742           //make fractional errors again
01743           sum_plus_tau = sum_plus_tau/Extrap[ie]->Pred[Background::kNuTauCC]->GetBinContent(ip+1,ir+1);
01744           sum_minus_tau = sum_minus_tau/Extrap[ie]->Pred[Background::kNuTauCC]->GetBinContent(ip+1,ir+1);
01745           sum_plus_nue = sum_plus_nue/Extrap[ie]->Pred[Background::kNueCC]->GetBinContent(ip+1,ir+1);
01746           sum_minus_nue = sum_minus_nue/Extrap[ie]->Pred[Background::kNueCC]->GetBinContent(ip+1,ir+1);
01747           
01748           NuTauCC_MC_Plus1Sigma[systname][ie]->SetBinContent(ip+1,ir+1,sum_plus_tau);
01749           NuTauCC_MC_Minus1Sigma[systname][ie]->SetBinContent(ip+1,ir+1,sum_minus_tau);
01750           NueCC_MC_Plus1Sigma[systname][ie]->SetBinContent(ip+1,ir+1,sum_plus_nue);
01751           NueCC_MC_Minus1Sigma[systname][ie]->SetBinContent(ip+1,ir+1,sum_minus_nue);
01752         }
01753       }
01754     }
01755   }
01756   
01757   if(flag==2)
01758   {
01759     //create empty histograms for numuCC shift
01760     for(ie=0;ie<Extrap.size();ie++)
01761     {
01762       Pred_CC_Fid_Plus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Plus1",systname.c_str(),ie+1),"",nt,t));
01763       Pred_CC_Fid_Minus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Minus1",systname.c_str(),ie+1),"",nt,t));
01764     }
01765     
01766     if(gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[0].c_str())) && gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[1].c_str())) && gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[2].c_str())))//if none of these files exist
01767     {
01768       for(ie=0;ie<Extrap.size();ie++)
01769       {
01770         NueCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01771         NueCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01772         NuTauCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01773         NuTauCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01774       }
01775       
01776       return;
01777     }
01778     
01779     for(j=0;j<app.size();j++)
01780     {
01781       nfiles=0;
01782       for(ie=0;ie<3;ie++)
01783       {
01784         if(gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[ie].c_str())))//if file doesn't exist
01785         {
01786           continue;
01787         }
01788         
01789         fpre = new TFile(gSystem->ExpandPathName(PreFiles[ie].c_str()),"READ");
01790         treepre = (TTree*)fpre->Get("paramtree");
01791         treepre->SetBranchAddress("nearPOT",&npot_pre);
01792         treepre->SetBranchAddress("farPOT",&fpot_pre);
01793         treepre->GetEntry(0);
01794 
01795         nPOTFar=Extrap[ie]->GetFarPOT();
01796         nPOTNear=Extrap[ie]->GetNearPOT();
01797         
01798         name = string(Background::AsString(app[j])) + "_" + std + "_Presel/FD_RecoVs" + Extrap[0]->GetPID();
01799         prestd = (TH2D*)fpre->Get(name.c_str());
01800         prestd->Scale(nPOTFar/fpot_pre);
01801         MBH.Rebin2DHist(prestd,np,p,nr,r);
01802         if(nfiles==0) totalprestd = (TH2D*)prestd->Clone();
01803         else totalprestd->Add(prestd);
01804         
01805         name = string(Background::AsString(app[j])) + "_" + plf + "_Presel/FD_RecoVs" + Extrap[0]->GetPID();
01806         if(plf==std)
01807         {
01808           prepl = (TH2D*)prestd->Clone(name.c_str());
01809         }
01810         else
01811         {
01812           prepl = (TH2D*)fpre->Get(name.c_str());
01813           prepl->Scale(nPOTFar/fpot_pre);
01814           MBH.Rebin2DHist(prepl,np,p,nr,r);
01815         }
01816         if(nfiles==0) totalprepl = (TH2D*)prepl->Clone();
01817         else totalprepl->Add(prepl);
01818         
01819         name = string(Background::AsString(app[j])) + "_" + mnf + "_Presel/FD_RecoVs" + Extrap[0]->GetPID();
01820         if(mnf==std)
01821         {
01822           premn = (TH2D*)prestd->Clone(name.c_str());
01823         }
01824         else if(mnf==plf)
01825         {
01826           premn = (TH2D*)prepl->Clone(name.c_str());
01827         }
01828         else
01829         {
01830           premn = (TH2D*)fpre->Get(name.c_str());
01831           premn->Scale(nPOTFar/fpot_pre);
01832           MBH.Rebin2DHist(premn,np,p,nr,r);
01833         }
01834         if(nfiles==0) totalpremn = (TH2D*)premn->Clone();
01835         else totalpremn->Add(premn);
01836         
01837         nfiles++;
01838       }
01839       
01840       totalprepl->Add(totalprestd,-1.);
01841       totalprepl->Divide(totalprestd);
01842       totalpremn->Add(totalprestd,-1.);
01843       totalpremn->Divide(totalprestd);
01844       
01845       if(plf==mnf && pln==mnn)
01846       {
01847         totalpremn->Scale(-1.);
01848       }
01849       
01850       for(ie=0;ie<Extrap.size();ie++)
01851       {
01852         if(app[j]==Background::kNueCC)
01853         {
01854           NueCC_MC_Plus1Sigma[systname].push_back((TH2D*)totalprepl->Clone(Form("%s_%i_NueCC_MC_Plus1",systname.c_str(),ie+1)));
01855           NueCC_MC_Minus1Sigma[systname].push_back((TH2D*)totalpremn->Clone(Form("%s_%i_NueCC_MC_Minus1",systname.c_str(),ie+1)));
01856         }
01857         else if(app[j]==Background::kNuTauCC)
01858         {
01859           NuTauCC_MC_Plus1Sigma[systname].push_back((TH2D*)totalprepl->Clone(Form("%s_%i_NuTauCC_MC_Plus1",systname.c_str(),ie+1)));
01860           NuTauCC_MC_Minus1Sigma[systname].push_back((TH2D*)totalpremn->Clone(Form("%s_%i_NuTauCC_MC_Minus1",systname.c_str(),ie+1)));
01861         }
01862       }
01863     }
01864   }
01865   
01866   return;
01867 }

void ErrorCalc::ReadSysFiles_FNExtrap ( int  n  )  [protected]

Definition at line 307 of file ErrorCalc.cxx.

References Background::AsString(), Extrap, FD_BNueCC_Minus1Sigma, FD_BNueCC_Plus1Sigma, FD_NC_Minus1Sigma, FD_NC_Plus1Sigma, FD_NuMuCC_Minus1Sigma, FD_NuMuCC_Plus1Sigma, FN_BNueCC_Minus1Sigma, FN_BNueCC_Plus1Sigma, FN_NC_Minus1Sigma, FN_NC_Plus1Sigma, FN_NuMuCC_Minus1Sigma, FN_NuMuCC_Plus1Sigma, FNExtrap_FarMinusTag, FNExtrap_FarPlusTag, FNExtrap_FileTag, FNExtrap_Flag, FNExtrap_NearMinusTag, FNExtrap_NearPlusTag, FNExtrap_StdTag, FNExtrap_SystName, Form(), OscFit::GetRunPeriod(), gSystem(), Background::kBNueCC, Background::kNC, Background::kNuMuCC, MBH, ND_BNueCC_Minus1Sigma, ND_BNueCC_Plus1Sigma, ND_NC_Minus1Sigma, ND_NC_Plus1Sigma, ND_NuMuCC_Minus1Sigma, ND_NuMuCC_Plus1Sigma, MultiBinAnaHelper::Rebin2DHist(), and SysFileDir.

Referenced by InitializeSys().

00308 {
00309   int i;
00310   
00311   string plf,pln,mnf,mnn,std,systname,filetag;
00312   systname = FNExtrap_SystName.at(n);
00313   plf = FNExtrap_FarPlusTag.at(n);
00314   pln = FNExtrap_NearPlusTag.at(n);
00315   mnf = FNExtrap_FarMinusTag.at(n);
00316   mnn = FNExtrap_NearMinusTag.at(n);
00317   std = FNExtrap_StdTag.at(n);
00318   filetag=FNExtrap_FileTag.at(n);
00319   int flag = FNExtrap_Flag.at(n);
00320   
00321   vector<string> PreFiles;
00322   
00323   unsigned int ie;
00324   for(ie=0;ie<Extrap.size();ie++)
00325   {
00326     PreFiles.push_back(SysFileDir+"/SysFile2D_"+filetag+"_Pre_"+Form("%i",Extrap[ie]->GetRunPeriod())+".root");
00327     if(gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[ie].c_str())))
00328     {
00329       if(flag!=1) cout<<"Warning: "<<PreFiles[ie]<<" doesn't exist."<<endl;
00330     }
00331   }
00332   
00333   unsigned int j;
00334   vector<Background::Background_t> bgs;
00335   bgs.push_back(Background::kNC);
00336   bgs.push_back(Background::kNuMuCC);
00337   bgs.push_back(Background::kBNueCC);
00338   
00339   string name;
00340   int np,nr;
00341   double *r,*p;
00342   TFile *fpre;
00343   TH2D *farstd,*farpl,*farmn;
00344   TH2D *nearstd,*nearpl,*nearmn;
00345   TH2D *fnpl, *fnmn;
00346   double npot_pre,fpot_pre;
00347   TTree *treepre;
00348   double nPOTNear,nPOTFar;
00349   
00350   for(ie=0;ie<Extrap.size();ie++)
00351   {
00352     nr = Extrap[ie]->GetNReco();
00353     np = Extrap[ie]->GetNPID();
00354     r = new double[nr+1];
00355     p = new double[np+1];
00356     for(i=0;i<nr+1;i++)
00357     {
00358       r[i] = Extrap[ie]->Pred_TotalBkgd->GetYaxis()->GetBinLowEdge(i+1);
00359     }
00360     for(i=0;i<np+1;i++)
00361     {
00362       p[i] = Extrap[ie]->Pred_TotalBkgd->GetXaxis()->GetBinLowEdge(i+1);
00363     }
00364     nPOTFar=Extrap[ie]->GetFarPOT();
00365     nPOTNear=Extrap[ie]->GetNearPOT();
00366     
00367     //flag==1 means systematics related to numuCC extrapolation for signal/tau prediction - by definition 0 for these components
00368     if(flag==1 || gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[ie].c_str())))//if Pre file does NOT exist
00369     {
00370       FN_NC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00371       FN_NC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00372       
00373       FN_NuMuCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NuMuCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00374       FN_NuMuCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NuMuCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00375       
00376       FN_BNueCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_BNueCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00377       FN_BNueCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_BNueCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00378       
00379       ND_NC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00380       ND_NC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00381       
00382       ND_NuMuCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NuMuCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00383       ND_NuMuCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NuMuCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00384       
00385       ND_BNueCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_BNueCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00386       ND_BNueCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_BNueCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00387       
00388       FD_NC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00389       FD_NC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00390       
00391       FD_NuMuCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NuMuCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00392       FD_NuMuCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NuMuCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00393       
00394       FD_BNueCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_BNueCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00395       FD_BNueCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_BNueCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00396       
00397       continue;
00398     }
00399     
00400     fpre = new TFile(gSystem->ExpandPathName(PreFiles[ie].c_str()),"READ");
00401     
00402     treepre = (TTree*)fpre->Get("paramtree");
00403     treepre->SetBranchAddress("nearPOT",&npot_pre);
00404     treepre->SetBranchAddress("farPOT",&fpot_pre);
00405     treepre->GetEntry(0);
00406     
00407     for(j=0;j<bgs.size();j++)
00408     {
00409       name = string(Background::AsString(bgs[j])) + "_" + std + "_Presel/FD_RecoVs" + Extrap[0]->GetPID();
00410       farstd = (TH2D*)fpre->Get(name.c_str());
00411       farstd->Scale(nPOTFar/fpot_pre);
00412       MBH.Rebin2DHist(farstd,np,p,nr,r);
00413       
00414       name = string(Background::AsString(bgs[j])) + "_" + plf + "_Presel/FD_RecoVs" + Extrap[0]->GetPID();
00415       if(plf==std)
00416       {
00417         farpl = (TH2D*)farstd->Clone(name.c_str());
00418       }
00419       else
00420       {
00421         farpl = (TH2D*)fpre->Get(name.c_str());
00422         farpl->Scale(nPOTFar/fpot_pre);
00423         MBH.Rebin2DHist(farpl,np,p,nr,r);
00424       }
00425       
00426       name = string(Background::AsString(bgs[j])) + "_" + mnf + "_Presel/FD_RecoVs" + Extrap[0]->GetPID();
00427       if(mnf==std)
00428       {
00429         farmn = (TH2D*)farstd->Clone(name.c_str());
00430       }
00431       else if(mnf==plf)
00432       {
00433         farmn = (TH2D*)farpl->Clone(name.c_str());
00434       }
00435       else
00436       {
00437         farmn = (TH2D*)fpre->Get(name.c_str());
00438         farmn->Scale(nPOTFar/fpot_pre);
00439         MBH.Rebin2DHist(farmn,np,p,nr,r);
00440       }
00441       
00442       name = string(Background::AsString(bgs[j])) + "_" + std + "_Presel/ND_RecoVs" + Extrap[0]->GetPID();
00443       nearstd = (TH2D*)fpre->Get(name.c_str());
00444       nearstd->Scale(nPOTNear/npot_pre);
00445       MBH.Rebin2DHist(nearstd,np,p,nr,r);
00446       
00447       name = string(Background::AsString(bgs[j])) + "_" + pln + "_Presel/ND_RecoVs" + Extrap[0]->GetPID();
00448       if(pln==std)
00449       {
00450         nearpl = (TH2D*)nearstd->Clone(name.c_str());
00451       }
00452       else
00453       {
00454         nearpl = (TH2D*)fpre->Get(name.c_str());
00455         nearpl->Scale(nPOTNear/npot_pre);
00456         MBH.Rebin2DHist(nearpl,np,p,nr,r);
00457       }
00458       
00459       name = string(Background::AsString(bgs[j])) + "_" + mnn + "_Presel/ND_RecoVs" + Extrap[0]->GetPID();
00460       if(mnn==std)
00461       {
00462         nearmn=(TH2D*)nearstd->Clone(name.c_str());
00463       }
00464       else if(mnn==pln)
00465       {
00466         nearmn=(TH2D*)nearpl->Clone(name.c_str());
00467       }
00468       else
00469       {
00470         nearmn = (TH2D*)fpre->Get(name.c_str());
00471         nearmn->Scale(nPOTNear/npot_pre);
00472         MBH.Rebin2DHist(nearmn,np,p,nr,r);
00473       }
00474       
00475       if(bgs[j]==Background::kNC)
00476       {
00477         ND_NC_Plus1Sigma[systname].push_back((TH2D*)nearpl->Clone(Form("%s_%i_ND_NC_Plus1",systname.c_str(),ie+1)));
00478         ND_NC_Plus1Sigma[systname][ie]->Add(nearstd,-1.0);
00479         ND_NC_Plus1Sigma[systname][ie]->Divide(nearstd);
00480         ND_NC_Minus1Sigma[systname].push_back((TH2D*)nearmn->Clone(Form("%s_%i_ND_NC_Minus1",systname.c_str(),ie+1)));
00481         ND_NC_Minus1Sigma[systname][ie]->Add(nearstd,-1.0);
00482         ND_NC_Minus1Sigma[systname][ie]->Divide(nearstd);
00483         if(plf==mnf && pln==mnn)
00484         {
00485           ND_NC_Minus1Sigma[systname][ie]->Scale(-1.);
00486         }
00487         
00488         FD_NC_Plus1Sigma[systname].push_back((TH2D*)farpl->Clone(Form("%s_%i_FD_NC_Plus1",systname.c_str(),ie+1)));
00489         FD_NC_Plus1Sigma[systname][ie]->Add(farstd,-1.0);
00490         FD_NC_Plus1Sigma[systname][ie]->Divide(farstd);
00491         FD_NC_Minus1Sigma[systname].push_back((TH2D*)farmn->Clone(Form("%s_%i_FD_NC_Minus1",systname.c_str(),ie+1)));
00492         FD_NC_Minus1Sigma[systname][ie]->Add(farstd,-1.0);
00493         FD_NC_Minus1Sigma[systname][ie]->Divide(farstd);
00494         if(plf==mnf && pln==mnn)
00495         {
00496           FD_NC_Minus1Sigma[systname][ie]->Scale(-1.);
00497         }
00498       }
00499       else if(bgs[j]==Background::kNuMuCC)
00500       {
00501         ND_NuMuCC_Plus1Sigma[systname].push_back((TH2D*)nearpl->Clone(Form("%s_%i_ND_NuMuCC_Plus1",systname.c_str(),ie+1)));
00502         ND_NuMuCC_Plus1Sigma[systname][ie]->Add(nearstd,-1.0);
00503         ND_NuMuCC_Plus1Sigma[systname][ie]->Divide(nearstd);
00504         ND_NuMuCC_Minus1Sigma[systname].push_back((TH2D*)nearmn->Clone(Form("%s_%i_ND_NuMuCC_Minus1",systname.c_str(),ie+1)));
00505         ND_NuMuCC_Minus1Sigma[systname][ie]->Add(nearstd,-1.0);
00506         ND_NuMuCC_Minus1Sigma[systname][ie]->Divide(nearstd);
00507         if(plf==mnf && pln==mnn)
00508         {
00509           ND_NuMuCC_Minus1Sigma[systname][ie]->Scale(-1.);
00510         }
00511         
00512         FD_NuMuCC_Plus1Sigma[systname].push_back((TH2D*)farpl->Clone(Form("%s_%i_FD_NuMuCC_Plus1",systname.c_str(),ie+1)));
00513         FD_NuMuCC_Plus1Sigma[systname][ie]->Add(farstd,-1.0);
00514         FD_NuMuCC_Plus1Sigma[systname][ie]->Divide(farstd);
00515         FD_NuMuCC_Minus1Sigma[systname].push_back((TH2D*)farmn->Clone(Form("%s_%i_FD_NuMuCC_Minus1",systname.c_str(),ie+1)));
00516         FD_NuMuCC_Minus1Sigma[systname][ie]->Add(farstd,-1.0);
00517         FD_NuMuCC_Minus1Sigma[systname][ie]->Divide(farstd);
00518         if(plf==mnf && pln==mnn)
00519         {
00520           FD_NuMuCC_Minus1Sigma[systname][ie]->Scale(-1.);
00521         }
00522       }
00523       else if(bgs[j]==Background::kBNueCC)
00524       {
00525         ND_BNueCC_Plus1Sigma[systname].push_back((TH2D*)nearpl->Clone(Form("%s_%i_ND_BNueCC_Plus1",systname.c_str(),ie+1)));
00526         ND_BNueCC_Plus1Sigma[systname][ie]->Add(nearstd,-1.0);
00527         ND_BNueCC_Plus1Sigma[systname][ie]->Divide(nearstd);
00528         ND_BNueCC_Minus1Sigma[systname].push_back((TH2D*)nearmn->Clone(Form("%s_%i_ND_BNueCC_Minus1",systname.c_str(),ie+1)));
00529         ND_BNueCC_Minus1Sigma[systname][ie]->Add(nearstd,-1.0);
00530         ND_BNueCC_Minus1Sigma[systname][ie]->Divide(nearstd);
00531         if(plf==mnf && pln==mnn)
00532         {
00533           ND_BNueCC_Minus1Sigma[systname][ie]->Scale(-1.);
00534         }
00535         
00536         FD_BNueCC_Plus1Sigma[systname].push_back((TH2D*)farpl->Clone(Form("%s_%i_FD_BNueCC_Plus1",systname.c_str(),ie+1)));
00537         FD_BNueCC_Plus1Sigma[systname][ie]->Add(farstd,-1.0);
00538         FD_BNueCC_Plus1Sigma[systname][ie]->Divide(farstd);
00539         FD_BNueCC_Minus1Sigma[systname].push_back((TH2D*)farmn->Clone(Form("%s_%i_FD_BNueCC_Minus1",systname.c_str(),ie+1)));
00540         FD_BNueCC_Minus1Sigma[systname][ie]->Add(farstd,-1.0);
00541         FD_BNueCC_Minus1Sigma[systname][ie]->Divide(farstd);
00542         if(plf==mnf && pln==mnn)
00543         {
00544           FD_BNueCC_Minus1Sigma[systname][ie]->Scale(-1.);
00545         }
00546       }
00547       
00548       farstd->Divide(nearstd);
00549       farpl->Divide(nearpl);
00550       farmn->Divide(nearmn);
00551       
00552       fnpl = (TH2D*)farpl->Clone("fnpl");
00553       fnpl->Add(farstd,-1.);
00554       fnpl->Divide(farstd);
00555       
00556       fnmn = (TH2D*)farmn->Clone("fnmn");
00557       fnmn->Add(farstd,-1.);
00558       fnmn->Divide(farstd);
00559       
00560       if(plf==mnf && pln==mnn)
00561       {
00562         fnmn->Scale(-1.);
00563       }
00564       
00565       if(bgs[j]==Background::kNC)
00566       {
00567         FN_NC_Plus1Sigma[systname].push_back((TH2D*)fnpl->Clone(Form("%s_%i_FN_NC_Plus1",systname.c_str(),ie+1)));
00568         FN_NC_Minus1Sigma[systname].push_back((TH2D*)fnmn->Clone(Form("%s_%i_FN_NC_Minus1",systname.c_str(),ie+1)));
00569       }
00570       if(bgs[j]==Background::kNuMuCC)
00571       {
00572         FN_NuMuCC_Plus1Sigma[systname].push_back((TH2D*)fnpl->Clone(Form("%s_%i_FN_NuMuCC_Plus1",systname.c_str(),ie+1)));
00573         FN_NuMuCC_Minus1Sigma[systname].push_back((TH2D*)fnmn->Clone(Form("%s_%i_FN_NuMuCC_Minus1",systname.c_str(),ie+1)));
00574       }
00575       if(bgs[j]==Background::kBNueCC)
00576       {
00577         FN_BNueCC_Plus1Sigma[systname].push_back((TH2D*)fnpl->Clone(Form("%s_%i_FN_BNueCC_Plus1",systname.c_str(),ie+1)));
00578         FN_BNueCC_Minus1Sigma[systname].push_back((TH2D*)fnmn->Clone(Form("%s_%i_FN_BNueCC_Minus1",systname.c_str(),ie+1)));
00579       }
00580     }
00581   }
00582   
00583   return;
00584 }

void ErrorCalc::ReadSysFiles_FNExtrap_AllRuns ( int  n  )  [protected, virtual]

Reimplemented in ErrorCalc_Joint.

Definition at line 585 of file ErrorCalc.cxx.

References Background::AsString(), Extrap, FD_BNueCC_Minus1Sigma, FD_BNueCC_Plus1Sigma, FD_NC_Minus1Sigma, FD_NC_Plus1Sigma, FD_NuMuCC_Minus1Sigma, FD_NuMuCC_Plus1Sigma, FN_BNueCC_Minus1Sigma, FN_BNueCC_Plus1Sigma, FN_NC_Minus1Sigma, FN_NC_Plus1Sigma, FN_NuMuCC_Minus1Sigma, FN_NuMuCC_Plus1Sigma, FNExtrap_FarMinusTag, FNExtrap_FarPlusTag, FNExtrap_FileTag, FNExtrap_Flag, FNExtrap_NearMinusTag, FNExtrap_NearPlusTag, FNExtrap_StdTag, FNExtrap_SystName, Form(), gSystem(), Background::kBNueCC, Background::kNC, Background::kNuMuCC, MBH, ND_BNueCC_Minus1Sigma, ND_BNueCC_Plus1Sigma, ND_NC_Minus1Sigma, ND_NC_Plus1Sigma, ND_NuMuCC_Minus1Sigma, ND_NuMuCC_Plus1Sigma, MultiBinAnaHelper::Rebin2DHist(), and SysFileDir.

Referenced by InitializeSys().

00586 {
00587   int i;
00588   
00589   string plf,pln,mnf,mnn,std,systname,filetag;
00590   systname = FNExtrap_SystName.at(n);
00591   plf = FNExtrap_FarPlusTag.at(n);
00592   pln = FNExtrap_NearPlusTag.at(n);
00593   mnf = FNExtrap_FarMinusTag.at(n);
00594   mnn = FNExtrap_NearMinusTag.at(n);
00595   std = FNExtrap_StdTag.at(n);
00596   filetag=FNExtrap_FileTag.at(n);
00597   int flag = FNExtrap_Flag.at(n);
00598   
00599   vector<string> PreFiles;
00600   
00601   unsigned int ie;
00602   for(ie=0;ie<3;ie++)
00603   {
00604     PreFiles.push_back(SysFileDir+"/SysFile2D_"+filetag+"_Pre_"+Form("%i",ie+1)+".root");
00605     if(gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[ie].c_str())))
00606     {
00607       if(flag!=1) cout<<"Warning: "<<PreFiles[ie]<<" doesn't exist."<<endl;
00608     }
00609   }
00610   
00611   unsigned int j;
00612   vector<Background::Background_t> bgs;
00613   bgs.push_back(Background::kNC);
00614   bgs.push_back(Background::kNuMuCC);
00615   bgs.push_back(Background::kBNueCC);
00616   
00617   string name;
00618   int np,nr;
00619   double *r,*p;
00620   TFile *fpre;
00621   TH2D *farstd=0,*farpl=0,*farmn=0;
00622   TH2D *nearstd=0,*nearpl=0,*nearmn=0;
00623   TH2D *totalfarstd=0,*totalfarpl=0,*totalfarmn=0;
00624   TH2D *totalnearstd=0,*totalnearpl=0,*totalnearmn=0;
00625   TH2D *fnpl=0, *fnmn=0;
00626   double npot_pre,fpot_pre;
00627   TTree *treepre;
00628   double nPOTNear,nPOTFar;
00629   int nfiles=0;
00630   
00631   nr = Extrap[0]->GetNReco();
00632   np = Extrap[0]->GetNPID();
00633   r = new double[nr+1];
00634   p = new double[np+1];
00635   for(i=0;i<nr+1;i++)
00636   {
00637     r[i] = Extrap[0]->Pred_TotalBkgd->GetYaxis()->GetBinLowEdge(i+1);
00638   }
00639   for(i=0;i<np+1;i++)
00640   {
00641     p[i] = Extrap[0]->Pred_TotalBkgd->GetXaxis()->GetBinLowEdge(i+1);
00642   }
00643   
00644   //flag==1 means systematics related to numuCC extrapolation for signal/tau prediction - by definition 0 for these components
00645   if(flag==1 || (gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[0].c_str())) &&  gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[1].c_str())) && gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[2].c_str()))))//if all three Pre files do NOT exist
00646   {
00647     for(ie=0;ie<Extrap.size();ie++)
00648     {
00649       FN_NC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00650       FN_NC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00651       
00652       FN_NuMuCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NuMuCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00653       FN_NuMuCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NuMuCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00654       
00655       FN_BNueCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_BNueCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00656       FN_BNueCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_BNueCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00657       
00658       ND_NC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00659       ND_NC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00660       
00661       ND_NuMuCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NuMuCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00662       ND_NuMuCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NuMuCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00663       
00664       ND_BNueCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_BNueCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00665       ND_BNueCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_BNueCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00666       
00667       FD_NC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00668       FD_NC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00669       
00670       FD_NuMuCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NuMuCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00671       FD_NuMuCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NuMuCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00672       
00673       FD_BNueCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_BNueCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00674       FD_BNueCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_BNueCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00675     }
00676     
00677     return;
00678   }
00679   
00680   for(j=0;j<bgs.size();j++)
00681   {
00682     nfiles=0;
00683     for(ie=0;ie<3;ie++)
00684     {
00685       if(gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[ie].c_str())))//if file doesn't exist
00686       {
00687         continue;
00688       }
00689       fpre = new TFile(gSystem->ExpandPathName(PreFiles[ie].c_str()),"READ");
00690       
00691       treepre = (TTree*)fpre->Get("paramtree");
00692       treepre->SetBranchAddress("nearPOT",&npot_pre);
00693       treepre->SetBranchAddress("farPOT",&fpot_pre);
00694       treepre->GetEntry(0);
00695       
00696       nPOTFar=Extrap[ie]->GetFarPOT();
00697       nPOTNear=Extrap[ie]->GetNearPOT();
00698       
00699       name = string(Background::AsString(bgs[j])) + "_" + std + "_Presel/FD_RecoVs" + Extrap[0]->GetPID();
00700       farstd = (TH2D*)fpre->Get(name.c_str());
00701       farstd->Scale(nPOTFar/fpot_pre);
00702       MBH.Rebin2DHist(farstd,np,p,nr,r);
00703       if(nfiles==0) totalfarstd = (TH2D*)farstd->Clone();
00704       else totalfarstd->Add(farstd);
00705       
00706       name = string(Background::AsString(bgs[j])) + "_" + plf + "_Presel/FD_RecoVs" + Extrap[0]->GetPID();
00707       if(plf==std)
00708       {
00709         farpl = (TH2D*)farstd->Clone(name.c_str());
00710       }
00711       else
00712       {
00713         farpl = (TH2D*)fpre->Get(name.c_str());
00714         farpl->Scale(nPOTFar/fpot_pre);
00715         MBH.Rebin2DHist(farpl,np,p,nr,r);
00716       }
00717       if(nfiles==0) totalfarpl = (TH2D*)farpl->Clone();
00718       else totalfarpl->Add(farpl);
00719       
00720       name = string(Background::AsString(bgs[j])) + "_" + mnf + "_Presel/FD_RecoVs" + Extrap[0]->GetPID();
00721       if(mnf==std)
00722       {
00723         farmn = (TH2D*)farstd->Clone(name.c_str());
00724       }
00725       else if(mnf==plf)
00726       {
00727         farmn = (TH2D*)farpl->Clone(name.c_str());
00728       }
00729       else
00730       {
00731         farmn = (TH2D*)fpre->Get(name.c_str());
00732         farmn->Scale(nPOTFar/fpot_pre);
00733         MBH.Rebin2DHist(farmn,np,p,nr,r);
00734       }
00735       if(nfiles==0) totalfarmn = (TH2D*)farmn->Clone();
00736       else totalfarmn->Add(farmn);
00737       
00738       name = string(Background::AsString(bgs[j])) + "_" + std + "_Presel/ND_RecoVs" + Extrap[0]->GetPID();
00739       nearstd = (TH2D*)fpre->Get(name.c_str());
00740       nearstd->Scale(nPOTNear/npot_pre);
00741       MBH.Rebin2DHist(nearstd,np,p,nr,r);
00742       if(nfiles==0) totalnearstd = (TH2D*)nearstd->Clone();
00743       else totalnearstd->Add(nearstd);
00744       
00745       name = string(Background::AsString(bgs[j])) + "_" + pln + "_Presel/ND_RecoVs" + Extrap[0]->GetPID();
00746       if(pln==std)
00747       {
00748         nearpl = (TH2D*)nearstd->Clone(name.c_str());
00749       }
00750       else
00751       {
00752         nearpl = (TH2D*)fpre->Get(name.c_str());
00753         nearpl->Scale(nPOTNear/npot_pre);
00754         MBH.Rebin2DHist(nearpl,np,p,nr,r);
00755       }
00756       if(nfiles==0) totalnearpl = (TH2D*)nearpl->Clone();
00757       else totalnearpl->Add(nearpl);
00758       
00759       name = string(Background::AsString(bgs[j])) + "_" + mnn + "_Presel/ND_RecoVs" + Extrap[0]->GetPID();
00760       if(mnn==std)
00761       {
00762         nearmn=(TH2D*)nearstd->Clone(name.c_str());
00763       }
00764       else if(mnn==pln)
00765       {
00766         nearmn=(TH2D*)nearpl->Clone(name.c_str());
00767       }
00768       else
00769       {
00770         nearmn = (TH2D*)fpre->Get(name.c_str());
00771         nearmn->Scale(nPOTNear/npot_pre);
00772         MBH.Rebin2DHist(nearmn,np,p,nr,r);
00773       }
00774       if(nfiles==0) totalnearmn = (TH2D*)nearmn->Clone();
00775       else totalnearmn->Add(nearmn);
00776       
00777       nfiles++;
00778     }
00779     
00780     for(ie=0;ie<Extrap.size();ie++)
00781     {
00782       if(bgs[j]==Background::kNC)
00783       {
00784         ND_NC_Plus1Sigma[systname].push_back((TH2D*)totalnearpl->Clone(Form("%s_ND_NC_Plus1",systname.c_str())));
00785         ND_NC_Plus1Sigma[systname][ie]->Add(totalnearstd,-1.0);
00786         ND_NC_Plus1Sigma[systname][ie]->Divide(totalnearstd);
00787         ND_NC_Minus1Sigma[systname].push_back((TH2D*)totalnearmn->Clone(Form("%s_ND_NC_Minus1",systname.c_str())));
00788         ND_NC_Minus1Sigma[systname][ie]->Add(totalnearstd,-1.0);
00789         ND_NC_Minus1Sigma[systname][ie]->Divide(totalnearstd);
00790         if(plf==mnf && pln==mnn)
00791         {
00792           ND_NC_Minus1Sigma[systname][ie]->Scale(-1.);
00793         }
00794         
00795         FD_NC_Plus1Sigma[systname].push_back((TH2D*)totalfarpl->Clone(Form("%s_FD_NC_Plus1",systname.c_str())));
00796         FD_NC_Plus1Sigma[systname][ie]->Add(totalfarstd,-1.0);
00797         FD_NC_Plus1Sigma[systname][ie]->Divide(totalfarstd);
00798         FD_NC_Minus1Sigma[systname].push_back((TH2D*)totalfarmn->Clone(Form("%s_FD_NC_Minus1",systname.c_str())));
00799         FD_NC_Minus1Sigma[systname][ie]->Add(totalfarstd,-1.0);
00800         FD_NC_Minus1Sigma[systname][ie]->Divide(totalfarstd);
00801         if(plf==mnf && pln==mnn)
00802         {
00803           FD_NC_Minus1Sigma[systname][ie]->Scale(-1.);
00804         }
00805       }
00806       else if(bgs[j]==Background::kNuMuCC)
00807       {
00808         ND_NuMuCC_Plus1Sigma[systname].push_back((TH2D*)totalnearpl->Clone(Form("%s_ND_NuMuCC_Plus1",systname.c_str())));
00809         ND_NuMuCC_Plus1Sigma[systname][ie]->Add(totalnearstd,-1.0);
00810         ND_NuMuCC_Plus1Sigma[systname][ie]->Divide(totalnearstd);
00811         ND_NuMuCC_Minus1Sigma[systname].push_back((TH2D*)totalnearmn->Clone(Form("%s_ND_NuMuCC_Minus1",systname.c_str())));
00812         ND_NuMuCC_Minus1Sigma[systname][ie]->Add(totalnearstd,-1.0);
00813         ND_NuMuCC_Minus1Sigma[systname][ie]->Divide(totalnearstd);
00814         if(plf==mnf && pln==mnn)
00815         {
00816           ND_NuMuCC_Minus1Sigma[systname][ie]->Scale(-1.);
00817         }
00818         
00819         FD_NuMuCC_Plus1Sigma[systname].push_back((TH2D*)totalfarpl->Clone(Form("%s_FD_NuMuCC_Plus1",systname.c_str())));
00820         FD_NuMuCC_Plus1Sigma[systname][ie]->Add(totalfarstd,-1.0);
00821         FD_NuMuCC_Plus1Sigma[systname][ie]->Divide(totalfarstd);
00822         FD_NuMuCC_Minus1Sigma[systname].push_back((TH2D*)totalfarmn->Clone(Form("%s_FD_NuMuCC_Minus1",systname.c_str())));
00823         FD_NuMuCC_Minus1Sigma[systname][ie]->Add(totalfarstd,-1.0);
00824         FD_NuMuCC_Minus1Sigma[systname][ie]->Divide(totalfarstd);
00825         if(plf==mnf && pln==mnn)
00826         {
00827           FD_NuMuCC_Minus1Sigma[systname][ie]->Scale(-1.);
00828         }
00829       }
00830       else if(bgs[j]==Background::kBNueCC)
00831       {
00832         ND_BNueCC_Plus1Sigma[systname].push_back((TH2D*)totalnearpl->Clone(Form("%s_ND_BNueCC_Plus1",systname.c_str())));
00833         ND_BNueCC_Plus1Sigma[systname][ie]->Add(totalnearstd,-1.0);
00834         ND_BNueCC_Plus1Sigma[systname][ie]->Divide(totalnearstd);
00835         ND_BNueCC_Minus1Sigma[systname].push_back((TH2D*)totalnearmn->Clone(Form("%s_ND_BNueCC_Minus1",systname.c_str())));
00836         ND_BNueCC_Minus1Sigma[systname][ie]->Add(totalnearstd,-1.0);
00837         ND_BNueCC_Minus1Sigma[systname][ie]->Divide(totalnearstd);
00838         if(plf==mnf && pln==mnn)
00839         {
00840           ND_BNueCC_Minus1Sigma[systname][ie]->Scale(-1.);
00841         }
00842         
00843         FD_BNueCC_Plus1Sigma[systname].push_back((TH2D*)totalfarpl->Clone(Form("%s_FD_BNueCC_Plus1",systname.c_str())));
00844         FD_BNueCC_Plus1Sigma[systname][ie]->Add(totalfarstd,-1.0);
00845         FD_BNueCC_Plus1Sigma[systname][ie]->Divide(totalfarstd);
00846         FD_BNueCC_Minus1Sigma[systname].push_back((TH2D*)totalfarmn->Clone(Form("%s_FD_BNueCC_Minus1",systname.c_str())));
00847         FD_BNueCC_Minus1Sigma[systname][ie]->Add(totalfarstd,-1.0);
00848         FD_BNueCC_Minus1Sigma[systname][ie]->Divide(totalfarstd);
00849         if(plf==mnf && pln==mnn)
00850         {
00851           FD_BNueCC_Minus1Sigma[systname][ie]->Scale(-1.);
00852         }
00853       }
00854     }
00855     
00856     totalfarstd->Divide(totalnearstd);
00857     totalfarpl->Divide(totalnearpl);
00858     totalfarmn->Divide(totalnearmn);
00859     
00860     fnpl = (TH2D*)totalfarpl->Clone("fnpl");
00861     fnpl->Add(totalfarstd,-1.);
00862     fnpl->Divide(totalfarstd);
00863     
00864     fnmn = (TH2D*)totalfarmn->Clone("fnmn");
00865     fnmn->Add(totalfarstd,-1.);
00866     fnmn->Divide(totalfarstd);
00867     
00868     if(plf==mnf && pln==mnn)
00869     {
00870       fnmn->Scale(-1.);
00871     }
00872     
00873     for(ie=0;ie<Extrap.size();ie++)
00874     {
00875       if(bgs[j]==Background::kNC)
00876       {
00877         FN_NC_Plus1Sigma[systname].push_back((TH2D*)fnpl->Clone(Form("%s_FN_NC_Plus1",systname.c_str())));
00878         FN_NC_Minus1Sigma[systname].push_back((TH2D*)fnmn->Clone(Form("%s_FN_NC_Minus1",systname.c_str())));
00879       }
00880       if(bgs[j]==Background::kNuMuCC)
00881       {
00882         FN_NuMuCC_Plus1Sigma[systname].push_back((TH2D*)fnpl->Clone(Form("%s_FN_NuMuCC_Plus1",systname.c_str())));
00883         FN_NuMuCC_Minus1Sigma[systname].push_back((TH2D*)fnmn->Clone(Form("%s_FN_NuMuCC_Minus1",systname.c_str())));
00884       }
00885       if(bgs[j]==Background::kBNueCC)
00886       {
00887         FN_BNueCC_Plus1Sigma[systname].push_back((TH2D*)fnpl->Clone(Form("%s_FN_BNueCC_Plus1",systname.c_str())));
00888         FN_BNueCC_Minus1Sigma[systname].push_back((TH2D*)fnmn->Clone(Form("%s_FN_BNueCC_Minus1",systname.c_str())));
00889       }
00890     }
00891   }
00892   
00893   return;
00894 }

void ErrorCalc::SetGridPred ( int  nbins,
vector< vector< double > >  nc,
vector< vector< double > >  cc,
vector< vector< double > >  bnue,
vector< vector< double > >  tau,
vector< vector< double > >  sig 
) [virtual]

Reimplemented in ErrorCalc_Joint.

Definition at line 3860 of file ErrorCalc.cxx.

References NueConvention::bnue, Extrap, GridPred, Init, Initialize(), Background::kBNueCC, Background::kNC, Background::kNueCC, Background::kNuMuCC, and Background::kNuTauCC.

Referenced by NueFit2D::RunDataGrid(), NueFit2D_Joint::RunDataGrid(), NueFit2D::RunMultiBinFC(), NueFit2D::RunMultiBinFC_MHDeltaFit(), NueFit2D_Joint::RunMultiBinFC_MHDeltaFit(), NueFit2D_Joint::RunMultiBinPseudoExpts(), NueFit2D::RunMultiBinPseudoExpts(), NueFit2D_Joint::RunMultiBinPseudoExpts_MHDeltaFit(), and NueFit2D::RunMultiBinPseudoExpts_MHDeltaFit().

03861 {
03862   if(Extrap.size()==0)
03863   {
03864     cout<<"Error in ErrorCalc::SetGridPred(): You must add an Extrapolate2D object.  Quitting..."<<endl;
03865     return;
03866   }
03867   
03868   Initialize();
03869   if(!Init) return;
03870   
03871   int i;
03872   unsigned int ie;
03873   
03874   if(nbins!=GridPred[Background::kNC][0]->GetNbinsX())
03875   {
03876     cout<<"Warning in SetGridPred(): nbins != GridPred[]->GetNbinsX()!!!"<<endl;
03877   }
03878   
03879   for(ie=0;ie<Extrap.size();ie++)
03880   {
03881     GridPred[Background::kNC][ie]->Reset();
03882     GridPred[Background::kNuMuCC][ie]->Reset();
03883     GridPred[Background::kBNueCC][ie]->Reset();
03884     GridPred[Background::kNuTauCC][ie]->Reset();
03885     GridPred[Background::kNueCC][ie]->Reset();
03886     
03887     for(i=0;i<nbins;i++)
03888     {
03889       GridPred[Background::kNC][ie]->SetBinContent(i+1,nc[i][ie]);
03890       GridPred[Background::kNuMuCC][ie]->SetBinContent(i+1,cc[i][ie]);
03891       GridPred[Background::kBNueCC][ie]->SetBinContent(i+1,bnue[i][ie]);
03892       GridPred[Background::kNuTauCC][ie]->SetBinContent(i+1,tau[i][ie]);
03893       GridPred[Background::kNueCC][ie]->SetBinContent(i+1,sig[i][ie]);
03894     }
03895   }
03896   
03897   return;
03898 }

void ErrorCalc::SetNDDataFile ( string  s  )  [inline]

Definition at line 32 of file ErrorCalc.h.

References NDData_infile.

00032 { NDData_infile = s; };

void ErrorCalc::SetSysFileDir ( string  dir  )  [inline]

Definition at line 29 of file ErrorCalc.h.

References SysFileDir.

00029 { SysFileDir = dir; };

virtual void ErrorCalc::SetUseGrid ( bool  grid = false  )  [inline, virtual]

Member Data Documentation

std::map<string,TH2D*> ErrorCalc::ExtraCovariance [protected]
std::map<string,vector<TH2D*> > ErrorCalc::FD_BNueCC_Minus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::FD_BNueCC_Plus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::FD_NC_Minus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::FD_NC_Plus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::FD_NuMuCC_Minus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::FD_NuMuCC_Plus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::FN_BNueCC_Minus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::FN_BNueCC_Plus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::FN_NC_Minus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::FN_NC_Plus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::FN_NuMuCC_Minus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::FN_NuMuCC_Plus1Sigma [protected]
vector<string> ErrorCalc::FNExtrap_FileTag [protected]
vector<int> ErrorCalc::FNExtrap_Flag [protected]
vector<string> ErrorCalc::FNExtrap_StdTag [protected]
vector<string> ErrorCalc::FNExtrap_SystName [protected]

Definition at line 98 of file ErrorCalc.h.

Referenced by AddFNExtrapSyst(), and InitializeSys().

std::map<Background::Background_t, vector<TH1D*> > ErrorCalc::GridPred [protected]
Bool_t ErrorCalc::Init [protected]
Bool_t ErrorCalc::InitDecompSys [protected]
Bool_t ErrorCalc::InitSys [protected]

Definition at line 81 of file ErrorCalc.h.

Referenced by ErrorCalc(), ErrorCalc_Joint::ErrorCalc_Joint(), and InitializeSys().

std::map<string,vector<TH2D*> > ErrorCalc::ND_BNueCC_Minus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::ND_BNueCC_Plus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::ND_NC_Minus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::ND_NC_Plus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::ND_NuMuCC_Minus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::ND_NuMuCC_Plus1Sigma [protected]
TH2D* ErrorCalc::NDCovMatrix [protected]

Definition at line 88 of file ErrorCalc.h.

Referenced by CalculateHOOError(), and InitializeDecompSys().

string ErrorCalc::NDData_infile [protected]
std::map<string,vector<TH2D*> > ErrorCalc::NueCC_MC_Minus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::NueCC_MC_Plus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::NuTauCC_MC_Minus1Sigma [protected]
std::map<string,vector<TH2D*> > ErrorCalc::NuTauCC_MC_Plus1Sigma [protected]
std::map<string,vector<TH1D*> > ErrorCalc::Pred_CC_Fid_Minus1Sigma [protected]
std::map<string,vector<TH1D*> > ErrorCalc::Pred_CC_Fid_Plus1Sigma [protected]
vector<double> ErrorCalc::SpecialSyst_Value [protected]
std::vector<TH1D*> ErrorCalc::SysBkgd_Plus1Sigma

Definition at line 52 of file ErrorCalc.h.

Referenced by MakeFracError_Lists().

string ErrorCalc::SysFileDir [protected]
std::vector<TH1D*> ErrorCalc::SysSig_Plus1Sigma

Definition at line 53 of file ErrorCalc.h.

Referenced by MakeFracError_Lists().

bool ErrorCalc::UseGrid [protected]

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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1