ErrorCalc_Joint Class Reference

#include <ErrorCalc_Joint.h>

Inheritance diagram for ErrorCalc_Joint:
ErrorCalc

List of all members.

Public Member Functions

 ErrorCalc_Joint ()
virtual ~ErrorCalc_Joint ()
void AddExtrap (Extrapolate2D *E)
void AddExtrapFHC (Extrapolate2D *E)
void AddExtrapRHC (Extrapolate2D *E)
void AddCovarianceMatrix (string systname, string file, string histname, int flag, int isRHC)
void CalculateHOOError ()
void SetGridPred (int nbins, vector< vector< double > > nc, vector< vector< double > > cc, vector< vector< double > > bnue, vector< vector< double > > tau, vector< vector< double > > sig)

Private Member Functions

void Initialize ()
void ReadSysFiles_FNExtrap_AllRuns (int n)
void ReadSysFiles_Appearance_AllRuns (int n)
void ReadSpecialFiles (int n)
void ReadCovarianceFiles (int n)
void CalculateSystErrorMatrixExtrap ()
void CalculateSystErrorMatrixGrid ()

Private Attributes

vector< Extrapolate2D * > Extrap
vector< int > ExtrapType
std::map< string, TH2D * > ExtraCovariance
std::map< string, int > ExtraCovarianceTag
vector< int > ExtraCovariance_IsRHC

Detailed Description

Definition at line 22 of file ErrorCalc_Joint.h.


Constructor & Destructor Documentation

ErrorCalc_Joint::ErrorCalc_Joint (  ) 

Definition at line 5 of file ErrorCalc_Joint.cxx.

References ErrorCalc::Init, ErrorCalc::InitDecompSys, ErrorCalc::InitSys, and ErrorCalc::UseGrid.

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

ErrorCalc_Joint::~ErrorCalc_Joint (  )  [virtual]

Definition at line 14 of file ErrorCalc_Joint.cxx.

00015 {
00016 }


Member Function Documentation

void ErrorCalc_Joint::AddCovarianceMatrix ( string  systname,
string  file,
string  histname,
int  flag,
int  isRHC 
)

Definition at line 1274 of file ErrorCalc_Joint.cxx.

References ErrorCalc::ExtraCovariance_File, ErrorCalc::ExtraCovariance_Flag, ErrorCalc::ExtraCovariance_HistName, ExtraCovariance_IsRHC, and ErrorCalc::ExtraCovariance_SystName.

01275 {
01276   if(flag>7)
01277   {
01278     cout<<"Error in AddCovarianceMatrix(): 'flag' should be <=7.  Not adding this covariance matrix."<<endl;
01279     return;
01280   }
01281   
01282   ExtraCovariance_SystName.push_back(systname);
01283   ExtraCovariance_File.push_back(file);
01284   ExtraCovariance_HistName.push_back(histname);
01285   ExtraCovariance_Flag.push_back(flag);
01286   ExtraCovariance_IsRHC.push_back(isRHC);
01287   
01288   return;
01289 }

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

Reimplemented from ErrorCalc.

Definition at line 18 of file ErrorCalc_Joint.cxx.

00019 {
00020   cout << "For Joint Fit, use AddExtrapRHC or AddExtrapFHC instead.  Doing nothing..." << endl;
00021   return;
00022 }

void ErrorCalc_Joint::AddExtrapFHC ( Extrapolate2D E  ) 

Definition at line 24 of file ErrorCalc_Joint.cxx.

References Extrap, and ExtrapType.

00025 {
00026   Extrap.push_back(E);
00027   //Make sure to indicate the Extrapolation type (0 for FHC):
00028   ExtrapType.push_back(0);
00029   return;
00030 }

void ErrorCalc_Joint::AddExtrapRHC ( Extrapolate2D E  ) 

Definition at line 31 of file ErrorCalc_Joint.cxx.

References Extrap, and ExtrapType.

00032 {
00033   Extrap.push_back(E);
00034   //Make sure to indicate the Extrapolation type (1 for RHC):
00035   ExtrapType.push_back(1);
00036   return;
00037 }

void ErrorCalc_Joint::CalculateHOOError (  )  [virtual]

Reimplemented from ErrorCalc.

Definition at line 1851 of file ErrorCalc_Joint.cxx.

References ErrorCalc::CovMatrix_Decomp, ErrorCalc::Init, and Initialize().

01852 {
01853   //For the joint fit, we're just going to initialize this to 0.
01854   
01855   Initialize();
01856   if(!Init) return;
01857   
01858   CovMatrix_Decomp->Reset();
01859 
01860   
01861   return;
01862 }

void ErrorCalc_Joint::CalculateSystErrorMatrixExtrap (  )  [private, virtual]

Reimplemented from ErrorCalc.

Definition at line 1332 of file ErrorCalc_Joint.cxx.

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

01333 {
01334   //calculates the SYSTEMATIC part of the error matrix
01335   if(Extrap.size()==0)
01336   {
01337     cout<<"Error in ErrorCalc_Joint::CalculateSystErrorMatrix(): You must add an Extrapolate2D object.  Quitting..."<<endl;
01338     return;
01339   }
01340     
01341   Initialize();
01342 
01343   if(!Init) return;
01344   
01345   InitializeSys();
01346 
01347   //cout << "Calc matrix" << endl;
01348 
01349   string systname;
01350   
01351   unsigned int ie=0;
01352 
01353   int nPIDfhc=0;
01354   int nPIDrhc=0;
01355   int nRecofhc=0;
01356   int nRecorhc=0;
01357 
01358   for (ie=0; ie<Extrap.size(); ie++){
01359     if (ExtrapType[ie]==0) {
01360       nPIDfhc = Extrap[ie]->GetNPID();
01361       nRecofhc = Extrap[ie]->GetNReco();
01362     }
01363     if (ExtrapType[ie]==1) {
01364       nPIDrhc = Extrap[ie]->GetNPID();
01365       nRecorhc = Extrap[ie]->GetNReco();
01366     }
01367   }
01368 
01369   //int nPID = Extrap[0]->GetNPID();
01370   //int nbins = Extrap[0]->GetNReco()*Extrap[0]->GetNPID();
01371                                                                                                                    //int nt = Extrap[0]->GetNTrue();
01372   
01373   std::map<string, vector<TH2D*> >::iterator systiter = FN_NC_Plus1Sigma.begin();
01374   std::map<string, vector<TH2D*> >::iterator last  = FN_NC_Plus1Sigma.end();
01375   
01376   int i,j;
01377   int ip,ir,jp,jr;
01378   double element;
01379   double di,dj;
01380   bool addi, addj;
01381 
01382   CovMatrix->Reset();
01383   
01384   int nbinsFHC = nPIDfhc*nRecofhc;
01385   int nbinsRHC = nPIDrhc*nRecorhc;
01386 
01387   int nbins = nbinsFHC + nbinsRHC;
01388 
01389   for(i=0;i<nbins;i++)
01390   {
01391 
01392     if (i<nbinsFHC){
01393       ir = int(i/nPIDfhc);
01394       ip = i%nPIDfhc;
01395     }      
01396     if (i>=nbinsFHC){
01397       ir = int((i-nbinsFHC)/nPIDrhc);
01398       ip = (i-nbinsFHC)%nPIDrhc;
01399     }      
01400 
01401 
01402     for(j=0;j<nbins;j++)
01403       {
01404 
01405         if (j<nbinsFHC){
01406           jr = int(j/nPIDfhc);
01407           jp = j%nPIDfhc;
01408         }      
01409         if (j>=nbinsFHC){
01410           jr = int((j-nbinsFHC)/nPIDrhc);
01411           jp = (j-nbinsFHC)%nPIDrhc;
01412         }
01413         
01414         element = 0;
01415         systiter = FN_NC_Plus1Sigma.begin();
01416         while(systiter!=last)
01417           {
01418 
01419             systname = systiter->first;
01420             
01421             di=0;
01422             dj=0;
01423 
01424             for(ie=0;ie<Extrap.size();ie++)
01425               {
01426                 addi=false;
01427                 addj=false;
01428 
01429                 //Only add the element if bin is correct for that particular extrapolation
01430                 if ((ExtrapType[ie]==0&&i<nbinsFHC)||(ExtrapType[ie]==1&&i>=nbinsFHC)) addi = true;
01431                 if ((ExtrapType[ie]==0&&j<nbinsFHC)||(ExtrapType[ie]==1&&j>=nbinsFHC)) addj = true;
01432 
01433 
01434                 if (addi) 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);
01435                 
01436                 if (addj) 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);
01437                   
01438                 if(Extrap[ie]->GetFNforBeamNue())
01439                   {
01440                     if (addi) di += FN_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->Pred[Background::kBNueCC]->GetBinContent(ip+1,ir+1);
01441                     if (addj) dj += FN_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*Extrap[ie]->Pred[Background::kBNueCC]->GetBinContent(jp+1,jr+1);
01442                   }
01443                 else
01444                   {
01445                     if (addi) di += FD_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*Extrap[ie]->Pred[Background::kBNueCC]->GetBinContent(ip+1,ir+1);
01446                     if (addj) dj += FD_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*Extrap[ie]->Pred[Background::kBNueCC]->GetBinContent(jp+1,jr+1);
01447                   }
01448                 
01449               } 
01450             
01451             element += (di*dj);
01452             
01453             systiter++;
01454           }
01455         
01456 
01457         //Add an additional element (1% of stat error) to force matrix to be non-singular
01458         double totstat=0;
01459 
01460         for(ie=0;ie<Extrap.size();ie++)
01461           {
01462             addi=false;
01463             addj=false;
01464             
01465             //Only add the element if bin is correct for that particular extrapolation
01466             if ((ExtrapType[ie]==0&&i<nbinsFHC)||(ExtrapType[ie]==1&&i>=nbinsFHC)) addi = true;
01467             if ((ExtrapType[ie]==0&&j<nbinsFHC)||(ExtrapType[ie]==1&&j>=nbinsFHC)) addj = true;
01468             if ((i==j)&&addi&&addj) totstat += Extrap[ie]->Pred[Background::kNueCC]->GetBinContent(ip+1,ir+1) 
01469                                       + Extrap[ie]->Pred_TotalBkgd->GetBinContent(ip+1,ir+1);
01470           }
01471         
01472         element += (totstat*0.01);
01473         
01474         CovMatrix->SetBinContent(i+1,j+1,element);
01475 
01476       }
01477 
01478 
01479   }
01480   
01481   std::map<string,TH2D*>::iterator coviter = ExtraCovariance.begin();
01482   std::map<string,TH2D*>::iterator covlast  = ExtraCovariance.end();
01483   double ni,nj;
01484   int k=0;
01485   int flag;
01486   int ietype=-1;
01487 
01488   while(coviter!=covlast)
01489   {
01490     systname = coviter->first;
01491     flag = ExtraCovariance_Flag.at(k);
01492 
01493     ietype=ExtraCovarianceTag[systname];
01494 
01495     for(i=0;i<nbins;i++)
01496     {
01497       
01498       if (i<nbinsFHC){
01499         ir = int(i/nPIDfhc);
01500         ip = i%nPIDfhc;
01501       }      
01502       if (i>=nbinsFHC){
01503         ir = int((i-nbinsFHC)/nPIDrhc);
01504         ip = (i-nbinsFHC)%nPIDrhc;
01505       }      
01506 
01507       ni=0;
01508       for(ie=0;ie<Extrap.size();ie++)
01509       {
01510 
01511         //Only add the element if bin is correct for that particular extrapolation
01512         if ((ietype==0&&ExtrapType[ie]==0&&i<nbinsFHC)
01513             ||(ietype==1&&ExtrapType[ie]==1&&i>=nbinsFHC)) {
01514           
01515           if(flag==0 || flag==1 || flag==4 || flag==5)
01516             {
01517               ni+=Extrap[ie]->Pred[Background::kNC]->GetBinContent(ip+1,ir+1);
01518             }
01519           if(flag==0 || flag==1 || flag==4 || flag==6)
01520             {
01521               ni+=Extrap[ie]->Pred[Background::kNuMuCC]->GetBinContent(ip+1,ir+1);
01522             }
01523           if(flag==0 || flag==1 || flag==7)
01524             {
01525               ni+=Extrap[ie]->Pred[Background::kBNueCC]->GetBinContent(ip+1,ir+1);
01526             }
01527           if(flag==0 || flag==2)
01528             {
01529               ni+=Extrap[ie]->Pred[Background::kNuTauCC]->GetBinContent(ip+1,ir+1);
01530             }
01531           if(flag==0 || flag==3)
01532             {
01533               ni+=Extrap[ie]->Pred[Background::kNueCC]->GetBinContent(ip+1,ir+1);
01534             }
01535         }
01536       }
01537       
01538       for(j=0;j<nbins;j++)
01539       {
01540         
01541         if (j<nbinsFHC){
01542           jr = int(j/nPIDfhc);
01543           jp = j%nPIDfhc;
01544         }      
01545         if (j>=nbinsFHC){
01546           jr = int((j-nbinsFHC)/nPIDrhc);
01547           jp = (j-nbinsFHC)%nPIDrhc;
01548         }      
01549 
01550         nj=0;
01551           
01552         for(ie=0;ie<Extrap.size();ie++)
01553           {
01554             
01555             //Only add the element if bin is correct for that particular extrapolation
01556             if ((ietype==0&&ExtrapType[ie]==0&&j<nbinsFHC)
01557                 ||(ietype==1&&ExtrapType[ie]==1&&j>=nbinsFHC)) {
01558 
01559               if(flag==0 || flag==1 || flag==4 || flag==5)
01560                 {
01561                   nj+=Extrap[ie]->Pred[Background::kNC]->GetBinContent(jp+1,jr+1);
01562                 }
01563               if(flag==0 || flag==1 || flag==4 || flag==6)
01564                 {
01565                   nj+=Extrap[ie]->Pred[Background::kNuMuCC]->GetBinContent(jp+1,jr+1);
01566                 }
01567               if(flag==0 || flag==1 || flag==7)
01568                 {
01569                   nj+=Extrap[ie]->Pred[Background::kBNueCC]->GetBinContent(jp+1,jr+1);
01570                 }
01571               if(flag==0 || flag==2)
01572                 {
01573                   nj+=Extrap[ie]->Pred[Background::kNuTauCC]->GetBinContent(jp+1,jr+1);
01574                 }
01575               if(flag==0 || flag==3)
01576                 {
01577                   nj+=Extrap[ie]->Pred[Background::kNueCC]->GetBinContent(jp+1,jr+1);
01578                 }
01579             }
01580         }
01581 
01582         element = CovMatrix->GetBinContent(i+1,j+1);
01583         element += (ni*nj*ExtraCovariance[systname]->GetBinContent(i+1,j+1));
01584         CovMatrix->SetBinContent(i+1,j+1,element);
01585       }
01586     }
01587     
01588     coviter++;
01589     k++;
01590   }
01591  
01592   return;
01593 }

void ErrorCalc_Joint::CalculateSystErrorMatrixGrid (  )  [private, virtual]

Reimplemented from ErrorCalc.

Definition at line 1595 of file ErrorCalc_Joint.cxx.

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

01596 {
01597   //calculates the SYSTEMATIC part of the error matrix
01598   if(Extrap.size()==0)
01599   {
01600     cout<<"Error in ErrorCalc::CalculateSystErrorMatrixGrid(): You must add an Extrapolate2D object.  Quitting..."<<endl;
01601     return;
01602   }
01603 
01604   Initialize();
01605   if(!Init) return;
01606   
01607   InitializeSys();
01608 
01609   string systname;
01610   
01611   int nPIDfhc=0;
01612   int nPIDrhc=0;
01613   int nRecofhc=0;
01614   int nRecorhc=0;
01615   unsigned int ie=0;
01616 
01617   for (ie=0; ie<Extrap.size(); ie++){
01618 
01619     if (ExtrapType[ie]==0) {
01620       nPIDfhc = Extrap[ie]->GetNPID();
01621       nRecofhc = Extrap[ie]->GetNReco();
01622     }
01623     if (ExtrapType[ie]==1) {
01624       nPIDrhc = Extrap[ie]->GetNPID();
01625       nRecorhc = Extrap[ie]->GetNReco();
01626     }
01627   }
01628 
01629   int nbinsFHC = nPIDfhc*nRecofhc;
01630   int nbinsRHC = nPIDrhc*nRecorhc;
01631   
01632   int nbins = nbinsFHC + nbinsRHC;
01633 
01634   std::map<string, vector<TH2D*> >::iterator systiter = FN_NC_Plus1Sigma.begin();
01635   std::map<string, vector<TH2D*> >::iterator last  = FN_NC_Plus1Sigma.end();
01636   
01637 
01638   int i,j;
01639   int ii,jj;
01640   int ip,ir,jp,jr;
01641   double element;
01642   double di,dj;
01643   bool addi, addj;
01644   
01645 //  unsigned int ietmp=0;
01646 //  unsigned int ieFHC=0;
01647 //  unsigned int ieRHC=0;
01648   CovMatrix->Reset();
01649   
01650   for(ii=0;ii<nbins;ii++)
01651   {
01652 
01653     if (ii<nbinsFHC){
01654       i = ii;
01655       ir = int(ii/nPIDfhc);
01656       ip = ii%nPIDfhc;
01657     }      
01658     if (ii>=nbinsFHC){
01659       i = ii-nbinsFHC;
01660       ir = int((ii-nbinsFHC)/nPIDrhc);
01661       ip = (ii-nbinsFHC)%nPIDrhc;
01662     }      
01663 
01664     
01665     for(jj=0;jj<nbins;jj++)
01666     {
01667         if (jj<nbinsFHC){
01668           j = jj;
01669           jr = int(jj/nPIDfhc);
01670           jp = jj%nPIDfhc;
01671         }      
01672         if (jj>=nbinsFHC){
01673           j = jj-nbinsFHC;
01674           jr = int((jj-nbinsFHC)/nPIDrhc);
01675           jp = (jj-nbinsFHC)%nPIDrhc;
01676         }
01677       
01678       element = 0;
01679       systiter = FN_NC_Plus1Sigma.begin();
01680       while(systiter!=last)
01681       {
01682         systname = systiter->first;
01683         
01684         di=0;
01685         dj=0;
01686         for(ie=0;ie<Extrap.size();ie++)
01687         {
01688 
01689           addi=false;
01690           addj=false;
01691           
01692           //Only add the element if bin is correct for that particular extrapolation
01693           if ((ExtrapType[ie]==0&&ii<nbinsFHC)||(ExtrapType[ie]==1&&ii>=nbinsFHC)) addi = true;
01694           if ((ExtrapType[ie]==0&&jj<nbinsFHC)||(ExtrapType[ie]==1&&jj>=nbinsFHC)) addj = true;
01695 
01696 
01697           if (addi) 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);
01698           
01699           if (addj) 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);
01700           
01701           if(Extrap[ie]->GetFNforBeamNue())
01702           {
01703             if (addi) di += FN_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*GridPred[Background::kBNueCC][ie]->GetBinContent(i+1);
01704             
01705             if (addj) dj += FN_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*GridPred[Background::kBNueCC][ie]->GetBinContent(j+1);
01706             
01707           }
01708           else
01709           {
01710             if (addi) di += FD_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(ip+1,ir+1)*GridPred[Background::kBNueCC][ie]->GetBinContent(i+1);
01711             
01712             if (addj) dj += FD_BNueCC_Plus1Sigma[systname][ie]->GetBinContent(jp+1,jr+1)*GridPred[Background::kBNueCC][ie]->GetBinContent(j+1);
01713             
01714           }
01715         }
01716         
01717         element += (di*dj);
01718         
01719         systiter++;
01720       }
01721 
01722       //Add an additional element (1% of stat error) to force matrix to be non-singular
01723         double totstat=0;
01724 
01725         for(ie=0;ie<Extrap.size();ie++)
01726           {
01727             addi=false;
01728             addj=false;
01729             
01730             //Only add the element if bin is correct for that particular extrapolation
01731             if ((ExtrapType[ie]==0&&ii<nbinsFHC)||(ExtrapType[ie]==1&&ii>=nbinsFHC)) addi = true;
01732             if ((ExtrapType[ie]==0&&jj<nbinsFHC)||(ExtrapType[ie]==1&&jj>=nbinsFHC)) addj = true;
01733             if ((ii==jj)&&addi&&addj) totstat += GridPred[Background::kNC][ie]->GetBinContent(i+1) + GridPred[Background::kNuMuCC][ie]->GetBinContent(i+1) + GridPred[Background::kBNueCC][ie]->GetBinContent(i+1) + GridPred[Background::kNuTauCC][ie]->GetBinContent(i+1) + GridPred[Background::kNueCC][ie]->GetBinContent(i+1);
01734                                       
01735           }
01736         
01737         element += (totstat*0.01);
01738         
01739         CovMatrix->SetBinContent(ii+1,jj+1,element);
01740     }
01741   }
01742   
01743   std::map<string,TH2D*>::iterator coviter = ExtraCovariance.begin();
01744   std::map<string,TH2D*>::iterator covlast  = ExtraCovariance.end();
01745 
01746   double ni,nj;
01747   int k=0;
01748   int flag;
01749   int ietype=-1;
01750 
01751   k=0;
01752   while(coviter!=covlast)
01753   {
01754     systname = coviter->first;
01755     flag = ExtraCovariance_Flag.at(k);
01756 
01757     ietype=ExtraCovarianceTag[systname];
01758 
01759     for(ii=0;ii<nbins;ii++)
01760     {
01761       
01762       if (ii<nbinsFHC)  i = ii;
01763       if (ii>=nbinsFHC) i = ii-nbinsFHC;
01764       
01765       ni=0;
01766       for(ie=0;ie<Extrap.size();ie++)
01767       {
01768         //Only add the element if bin is correct for that particular extrapolation
01769         if ((ietype==0&&ExtrapType[ie]==0&&ii<nbinsFHC)
01770             ||(ietype==1&&ExtrapType[ie]==1&&ii>=nbinsFHC)) {
01771           
01772           if(flag==0 || flag==1 || flag==4 || flag==5)
01773             {
01774               ni+=GridPred[Background::kNC][ie]->GetBinContent(i+1);
01775             }
01776           if(flag==0 || flag==1 || flag==4 || flag==6)
01777             {
01778               ni+=GridPred[Background::kNuMuCC][ie]->GetBinContent(i+1);
01779             }
01780           if(flag==0 || flag==1 || flag==7)
01781             {
01782               ni+=GridPred[Background::kBNueCC][ie]->GetBinContent(i+1);
01783             }
01784           if(flag==0 || flag==2)
01785             {
01786               ni+=GridPred[Background::kNuTauCC][ie]->GetBinContent(i+1);
01787             }
01788           if(flag==0 || flag==3)
01789             {
01790               ni+=GridPred[Background::kNueCC][ie]->GetBinContent(i+1);
01791             }
01792         }
01793       }
01794       
01795       for(jj=0;jj<nbins;jj++)
01796         {
01797           
01798         if (jj<nbinsFHC)  j = jj;
01799         if (jj>=nbinsFHC) j = jj-nbinsFHC;
01800         
01801         nj=0;
01802         for(ie=0;ie<Extrap.size();ie++)
01803           {
01804             //Only add the element if bin is correct for that particular extrapolation
01805             if ((ietype==0&&ExtrapType[ie]==0&&jj<nbinsFHC)
01806                 ||(ietype==1&&ExtrapType[ie]==1&&jj>=nbinsFHC)) {
01807 
01808               if(flag==0 || flag==1 || flag==4 || flag==5)
01809                 {
01810                   nj+=GridPred[Background::kNC][ie]->GetBinContent(j+1);
01811                 }
01812               if(flag==0 || flag==1 || flag==4 || flag==6)
01813                 {
01814                   nj+=GridPred[Background::kNuMuCC][ie]->GetBinContent(j+1);
01815                 }
01816               if(flag==0 || flag==1 || flag==7)
01817                 {
01818                   nj+=GridPred[Background::kBNueCC][ie]->GetBinContent(j+1);
01819                 }
01820               if(flag==0 || flag==2)
01821                 {
01822                   nj+=GridPred[Background::kNuTauCC][ie]->GetBinContent(j+1);
01823                 }
01824               if(flag==0 || flag==3)
01825                 {
01826                   nj+=GridPred[Background::kNueCC][ie]->GetBinContent(j+1);
01827                 }
01828             }
01829           }
01830 
01831         //Added protection:
01832         if ((ietype==0&ii<nbinsFHC&&jj<nbinsFHC)
01833             ||(ietype==1&&ii>=nbinsFHC&&jj>=nbinsFHC)){
01834           element = CovMatrix->GetBinContent(ii+1,jj+1);
01835           element += (ni*nj*ExtraCovariance[systname]->GetBinContent(i+1,j+1));
01836           CovMatrix->SetBinContent(ii+1,jj+1,element);
01837         }
01838 
01839 
01840         }
01841     }
01842     
01843     coviter++;
01844     k++;
01845   }
01846   
01847   return;
01848 }

void ErrorCalc_Joint::Initialize (  )  [private, virtual]

Reimplemented from ErrorCalc.

Definition at line 40 of file ErrorCalc_Joint.cxx.

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

Referenced by CalculateHOOError(), CalculateSystErrorMatrixExtrap(), CalculateSystErrorMatrixGrid(), and SetGridPred().

00041 {
00042   if (Init) return;
00043   
00044   int nbinsFHC = 0;
00045   int nbinsRHC = 0;
00046 
00047   cout.precision(5);
00048   cout<<fixed;
00049   
00050   unsigned int ie;
00051   for(ie=0;ie<Extrap.size();ie++)
00052   {
00053     Extrap[ie]->GetPrediction();
00054     if(Extrap[ie]->GetOscFlag())
00055     {
00056       cout<<"Error in ErrorCalc_Joint::Initialize(): You have called Extrapolate2D::OscillatePrediction() before initializing the systematics!  Quitting..."<<endl;
00057       return;
00058     }
00059     //Get number of bins for each type:
00060     if (ExtrapType[ie]==0) nbinsFHC = Extrap[ie]->GetNReco()*Extrap[ie]->GetNPID();
00061     if (ExtrapType[ie]==1) nbinsRHC = Extrap[ie]->GetNReco()*Extrap[ie]->GetNPID();
00062   }
00063   
00064   gROOT->cd();
00065   int nbins = nbinsFHC + nbinsRHC;//Used to be Extrap[0]->GetNReco()*Extrap[0]->GetNPID();
00066   CovMatrix = new TH2D("CovMatrix","",nbins,-0.5,nbins-0.5,nbins,-0.5,nbins-0.5);
00067   CovMatrix_Decomp = new TH2D("CovMatrix_Decomp","",nbins,-0.5,nbins-0.5,nbins,-0.5,nbins-0.5);
00068   
00069 
00070   //Double-check compatibility with joint fit later, but this should be fine for now
00071   string st;
00072   for(ie=0;ie<Extrap.size();ie++)
00073   {
00074 
00075     if (ExtrapType[ie]==0) nbins = nbinsFHC;
00076     if (ExtrapType[ie]==1) nbins = nbinsRHC;
00077 
00078       st = string(Background::AsString(Background::kNC));  
00079       GridPred[Background::kNC].push_back(new TH1D(Form("GridPred_%s_%i",st.c_str(),ie),"",nbins,-0.5,nbins-0.5));
00080       st = string(Background::AsString(Background::kNuMuCC));  
00081       GridPred[Background::kNuMuCC].push_back(new TH1D(Form("GridPred_%s_%i",st.c_str(),ie),"",nbins,-0.5,nbins-0.5));
00082       st = string(Background::AsString(Background::kBNueCC));  
00083       GridPred[Background::kBNueCC].push_back(new TH1D(Form("GridPred_%s_%i",st.c_str(),ie),"",nbins,-0.5,nbins-0.5));
00084       st = string(Background::AsString(Background::kNuTauCC));  
00085       GridPred[Background::kNuTauCC].push_back(new TH1D(Form("GridPred_%s_%i",st.c_str(),ie),"",nbins,-0.5,nbins-0.5));
00086       st = string(Background::AsString(Background::kNueCC));  
00087       GridPred[Background::kNueCC].push_back(new TH1D(Form("GridPred_%s_%i",st.c_str(),ie),"",nbins,-0.5,nbins-0.5));
00088    }
00089 
00090   Init = true;
00091   
00092   return;
00093 }

void ErrorCalc_Joint::ReadCovarianceFiles ( int  n  )  [private, virtual]

Reimplemented from ErrorCalc.

Definition at line 1292 of file ErrorCalc_Joint.cxx.

References ExtraCovariance, ErrorCalc::ExtraCovariance_File, ErrorCalc::ExtraCovariance_HistName, ExtraCovariance_IsRHC, ErrorCalc::ExtraCovariance_SystName, ExtraCovarianceTag, Extrap, ExtrapType, and gSystem().

01293 {
01294 
01295   int matrixtype = ExtraCovariance_IsRHC.at(n);
01296 
01297   string systname = ExtraCovariance_SystName.at(n);
01298   string file = ExtraCovariance_File.at(n);
01299   string hist = ExtraCovariance_HistName.at(n);
01300 //   int flag = ExtraCovariance_Flag.at(n);
01301   
01302   if(gSystem->AccessPathName(gSystem->ExpandPathName(file.c_str())))
01303   {
01304     cout<<"Can't read "<<file<<endl;
01305     return;
01306   }
01307   
01308   TFile *f = new TFile(file.c_str(),"READ");
01309   TH2D *h = (TH2D*)f->Get(hist.c_str());
01310 
01311   unsigned int ie=0;
01312   int nbins=0;
01313 
01314   //Bins for different parts of joint fit:
01315   for (ie=0; ie<Extrap.size(); ie++){
01316     if (ExtrapType[ie]==0&&matrixtype==0) nbins = Extrap[ie]->GetNPID()*Extrap[ie]->GetNReco();
01317     if (ExtrapType[ie]==1&&matrixtype==1) nbins = Extrap[ie]->GetNPID()*Extrap[ie]->GetNReco();
01318   }
01319   
01320   if(h->GetNbinsX() != nbins)
01321   {
01322     cout<<"Error in ReadCovarianceFiles(): # bins from matrix in file doesn't match extrapolation bins!"<<endl;
01323     return;
01324   }
01325   
01326   ExtraCovariance[systname] = h;
01327   ExtraCovarianceTag[systname] = matrixtype;
01328   
01329   return;
01330 }

void ErrorCalc_Joint::ReadSpecialFiles ( int  n  )  [private, virtual]

Reimplemented from ErrorCalc.

Definition at line 1086 of file ErrorCalc_Joint.cxx.

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

01087 {
01088 
01089 
01090   string systname = SpecialSyst_SystName.at(n);
01091   TH2D *h;
01092   unsigned int ie;
01093   int i,nr,np,nt;
01094   double *r,*p,*t;
01095   
01096   //create empty histograms for all shifts
01097   for(ie=0;ie<Extrap.size();ie++)
01098   {
01099 
01100     nr = Extrap[ie]->GetNReco();
01101     np = Extrap[ie]->GetNPID();
01102     nt = Extrap[ie]->GetNTrue();
01103     r = new double[nr+1];
01104     p = new double[np+1];
01105     t = new double[nt+1];
01106     for(i=0;i<nr+1;i++)
01107       {
01108         r[i] = Extrap[ie]->Pred_TotalBkgd->GetYaxis()->GetBinLowEdge(i+1);
01109       }
01110     for(i=0;i<np+1;i++)
01111       {
01112         p[i] = Extrap[ie]->Pred_TotalBkgd->GetXaxis()->GetBinLowEdge(i+1);
01113       }
01114     for(i=0;i<nt+1;i++)
01115       {
01116         t[i] = Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNuTau]->GetZaxis()->GetBinLowEdge(i+1);
01117       }
01118 
01119     FN_NC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01120     FN_NC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01121     
01122     FN_NuMuCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NuMuCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01123     FN_NuMuCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NuMuCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01124     
01125     FN_BNueCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_BNueCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01126     FN_BNueCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_BNueCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01127     
01128     ND_NC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01129     ND_NC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01130     
01131     ND_NuMuCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NuMuCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01132     ND_NuMuCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NuMuCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01133     
01134     ND_BNueCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_BNueCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01135     ND_BNueCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_BNueCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01136     
01137     FD_NC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01138     FD_NC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01139     
01140     FD_NuMuCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NuMuCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01141     FD_NuMuCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NuMuCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01142     
01143     FD_BNueCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_BNueCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01144     FD_BNueCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_BNueCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01145     
01146     NueCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01147     NueCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01148     NuTauCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
01149     NuTauCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
01150     
01151     Pred_CC_Fid_Plus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Plus1",systname.c_str(),ie+1),"",nt,t));
01152     Pred_CC_Fid_Minus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Minus1",systname.c_str(),ie+1),"",nt,t));
01153   }
01154 
01155   int flag = SpecialSyst_Flag.at(n);
01156 
01157   if(SpecialSyst_Value.at(n)>0)//if value is set, use that number in all bins
01158   {
01159 
01160     for(unsigned int ie=0;ie<Extrap.size();ie++)
01161       {
01162         
01163         h = (TH2D*)Extrap[ie]->Pred_TotalBkgd->Clone("htemp");
01164         h->Reset();
01165         for(int ip=0;ip<Extrap[ie]->GetNPID();ip++)
01166           {
01167             for(int ir=0;ir<Extrap[ie]->GetNReco();ir++)
01168               {
01169                 h->SetBinContent(ip+1,ir+1,SpecialSyst_Value.at(n));
01170               }
01171           }
01172         
01173       if(flag==0 || flag==1)//for NC,NuMuCC,BNueCC
01174       {
01175         FN_NC_Plus1Sigma[systname][ie]->Add(h);
01176         FN_NC_Minus1Sigma[systname][ie]->Add(h,-1);
01177         
01178         FN_NuMuCC_Plus1Sigma[systname][ie]->Add(h);
01179         FN_NuMuCC_Minus1Sigma[systname][ie]->Add(h,-1);
01180         
01181         FN_BNueCC_Plus1Sigma[systname][ie]->Add(h);
01182         FN_BNueCC_Minus1Sigma[systname][ie]->Add(h,-1);
01183         
01184         FD_BNueCC_Plus1Sigma[systname][ie]->Add(h);
01185         FD_BNueCC_Minus1Sigma[systname][ie]->Add(h,-1);
01186       }
01187       if(flag==0 || flag==2)
01188       {
01189         NuTauCC_MC_Plus1Sigma[systname][ie]->Add(h);
01190         NuTauCC_MC_Minus1Sigma[systname][ie]->Add(h,-1);
01191       }
01192       if(flag==0 || flag==3)
01193       {
01194         NueCC_MC_Plus1Sigma[systname][ie]->Add(h);
01195         NueCC_MC_Minus1Sigma[systname][ie]->Add(h,-1);
01196       }
01197       if(flag==4)//for NC and NuMuCC (fake HOO)
01198       {
01199         FN_NC_Plus1Sigma[systname][ie]->Add(h);
01200         FN_NC_Minus1Sigma[systname][ie]->Add(h,-1);
01201         
01202         FN_NuMuCC_Plus1Sigma[systname][ie]->Add(h);
01203         FN_NuMuCC_Minus1Sigma[systname][ie]->Add(h,-1);
01204 
01205       }
01206 
01207 
01208     }
01209     
01210     return;
01211   }
01212   
01213   //if value is 0, read from a file
01214   string filetag,histname;
01215   filetag=SpecialSyst_FileTag.at(n);
01216   histname = SpecialSyst_HistName.at(n);
01217   
01218   TFile *f;
01219   
01220   for(unsigned int ie=0;ie<Extrap.size();ie++)
01221   {
01222     if(ExtrapType[ie]==0&&gSystem->AccessPathName(gSystem->ExpandPathName(Form("%s/%s_FHC%i.root",SysFileDir.c_str(),filetag.c_str(),Extrap[ie]->GetRunPeriod()))))
01223     {
01224       continue;
01225     }
01226     if(ExtrapType[ie]==1&&gSystem->AccessPathName(gSystem->ExpandPathName(Form("%s/%s_RHC%i.root",SysFileDir.c_str(),filetag.c_str(),Extrap[ie]->GetRunPeriod()))))
01227     {
01228       continue;
01229     }
01230     
01231     if (ExtrapType[ie]==0) f = new TFile(gSystem->ExpandPathName(Form("%s/%s_FHC%i.root",SysFileDir.c_str(),filetag.c_str(),Extrap[ie]->GetRunPeriod())),"READ");
01232     if (ExtrapType[ie]==1) f = new TFile(gSystem->ExpandPathName(Form("%s/%s_RHC%i.root",SysFileDir.c_str(),filetag.c_str(),Extrap[ie]->GetRunPeriod())),"READ");
01233   
01234     h = (TH2D*)f->Get(histname.c_str());
01235     
01236     if(flag==0 || flag==1)//for NC,NuMuCC,BNueCC
01237     {
01238       FN_NC_Plus1Sigma[systname][ie]->Add(h);
01239       FN_NC_Minus1Sigma[systname][ie]->Add(h,-1);
01240         
01241       FN_NuMuCC_Plus1Sigma[systname][ie]->Add(h);
01242       FN_NuMuCC_Minus1Sigma[systname][ie]->Add(h,-1);
01243         
01244       FN_BNueCC_Plus1Sigma[systname][ie]->Add(h);
01245       FN_BNueCC_Minus1Sigma[systname][ie]->Add(h,-1);
01246         
01247       FD_BNueCC_Plus1Sigma[systname][ie]->Add(h);
01248       FD_BNueCC_Minus1Sigma[systname][ie]->Add(h,-1);
01249     }
01250     if(flag==0 || flag==2)
01251     {
01252       NuTauCC_MC_Plus1Sigma[systname][ie]->Add(h);
01253       NuTauCC_MC_Minus1Sigma[systname][ie]->Add(h,-1);
01254     }
01255     if(flag==0 || flag==3)
01256     {
01257       NueCC_MC_Plus1Sigma[systname][ie]->Add(h);
01258       NueCC_MC_Minus1Sigma[systname][ie]->Add(h,-1);
01259     }
01260     if(flag==4)//for NC and NuMuCC (fake HOO)
01261       {
01262         FN_NC_Plus1Sigma[systname][ie]->Add(h);
01263         FN_NC_Minus1Sigma[systname][ie]->Add(h,-1);
01264 
01265         FN_NuMuCC_Plus1Sigma[systname][ie]->Add(h);
01266         FN_NuMuCC_Minus1Sigma[systname][ie]->Add(h,-1);
01267       }
01268 
01269   }
01270   
01271   return;
01272 }

void ErrorCalc_Joint::ReadSysFiles_Appearance_AllRuns ( int  n  )  [private, virtual]

Reimplemented from ErrorCalc.

Definition at line 451 of file ErrorCalc_Joint.cxx.

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

00452 {
00453 
00454   int i;
00455   int ietag=-1;
00456 
00457   string plf,pln,mnf,mnn,std,systname,filetag;
00458   systname = FNExtrap_SystName.at(n);
00459   plf = FNExtrap_FarPlusTag.at(n);
00460   pln = FNExtrap_NearPlusTag.at(n);
00461   mnf = FNExtrap_FarMinusTag.at(n);
00462   mnn = FNExtrap_NearMinusTag.at(n);
00463   std = FNExtrap_StdTag.at(n);
00464   filetag=FNExtrap_FileTag.at(n);
00465   int flag = FNExtrap_Flag.at(n);
00466   
00467   vector<string> FidFiles;
00468   vector<string> CClikeFiles;
00469   vector<string> PreFiles;
00470   
00471 
00472   unsigned int ie;
00473   for(ie=0;ie<Extrap.size();ie++)
00474   {
00475   
00476     if(flag==1)//systs related to NuMuCC-like prediction
00477     {
00478     if (ExtrapType[ie]==0) FidFiles.push_back(SysFileDir+"/SysFile2D_FHC"+filetag+"_Fid_1.root");
00479     if (ExtrapType[ie]==1) FidFiles.push_back(SysFileDir+"/SysFile2D_RHC"+filetag+"_Fid_1.root");
00480     if (ExtrapType[ie]==0) CClikeFiles.push_back(SysFileDir+"/SysFile2D_FHC"+filetag+"_CCLike_1.root");
00481     if (ExtrapType[ie]==1) CClikeFiles.push_back(SysFileDir+"/SysFile2D_RHC"+filetag+"_CCLike_1.root");
00482     }
00483     else if(flag==2)//uncertainty in selected tau or signal FD MC
00484     {
00485     if (ExtrapType[ie]==0) PreFiles.push_back(SysFileDir+"/SysFile2D_FHC"+filetag+"_Pre_1.root");
00486     if (ExtrapType[ie]==1) PreFiles.push_back(SysFileDir+"/SysFile2D_RHC"+filetag+"_Pre_1.root");
00487     
00488     }
00489   }
00490   
00491   unsigned int j;
00492   vector<Background::Background_t> app;
00493   app.push_back(Background::kNueCC);
00494   app.push_back(Background::kNuTauCC);
00495   
00496   string name;
00497   int np,nr,nt;
00498   double *r,*p,*t;
00499   TFile *fpre,*ffid,*fcc;
00500   double npot_pre,fpot_pre;
00501   double npot_fid,fpot_fid;
00502   double npot_cc,fpot_cc;
00503   TTree *treepre,*treefid,*treecc;
00504   double nPOTNear,nPOTFar;
00505   TH2D *prestd=0,*prepl=0,*premn=0;
00506   TH2D *totalprestd[2]={0},*totalprepl[2]={0},*totalpremn[2]={0};
00507   double sum,temp;
00508   int ir,it,ip;
00509   TH1D *ccstd[2]={0},*ccpl[2]={0},*ccmn[2]={0};
00510   TH1D *data=0;
00511   TH2D *farfidccstd=0,*farcclikestd=0,*nearcclikestd=0;
00512   TH2D *farfidccpl=0,*farcclikepl=0,*nearcclikepl=0;
00513   TH2D *farfidccmn=0,*farcclikemn=0,*nearcclikemn=0;
00514   TH1D *totaldata[2]={0};
00515   TH2D *totalfarfidccstd[2]={0},*totalfarcclikestd[2]={0},*totalnearcclikestd[2]={0};
00516   TH2D *totalfarfidccpl[2]={0},*totalfarcclikepl[2]={0},*totalnearcclikepl[2]={0};
00517   TH2D *totalfarfidccmn[2]={0},*totalfarcclikemn[2]={0},*totalnearcclikemn[2]={0};
00518   TH2D *h2;
00519   int nfiles[2]={0};
00520   double sum_plus_tau,sum_minus_tau,sum_plus_nue,sum_minus_nue;
00521   string ccname = "CC";
00522     
00523   if(flag==0)
00524   {
00525     for(ie=0;ie<Extrap.size();ie++)
00526     {
00527 
00528       nr = Extrap[ie]->GetNReco();
00529       np = Extrap[ie]->GetNPID();
00530       nt = Extrap[ie]->GetNTrue();
00531       r = new double[nr+1];
00532       p = new double[np+1];
00533       t = new double[nt+1];
00534       for(i=0;i<nr+1;i++)
00535         {
00536           r[i] = Extrap[ie]->Pred_TotalBkgd->GetYaxis()->GetBinLowEdge(i+1);
00537         }
00538       for(i=0;i<np+1;i++)
00539         {
00540           p[i] = Extrap[ie]->Pred_TotalBkgd->GetXaxis()->GetBinLowEdge(i+1);
00541         }
00542       for(i=0;i<nt+1;i++)
00543         {
00544           t[i] = Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNuTau]->GetZaxis()->GetBinLowEdge(i+1);
00545         }
00546 
00547 
00548        //create empty histograms for overall signal/tau shifts
00549       NueCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00550       NueCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00551       NuTauCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00552       NuTauCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00553       
00554        //creating empty numuCC shifts
00555       Pred_CC_Fid_Plus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Plus1",systname.c_str(),ie+1),"",nt,t));
00556       Pred_CC_Fid_Minus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Minus1",systname.c_str(),ie+1),"",nt,t));
00557 
00558     }
00559   }
00560   
00561   if(flag==1)
00562   {
00563     //create empty histograms for overall signal/tau shifts
00564     for(ie=0;ie<Extrap.size();ie++)
00565     {
00566 
00567       nr = Extrap[ie]->GetNReco();
00568       np = Extrap[ie]->GetNPID();
00569       r = new double[nr+1];
00570       p = new double[np+1];
00571 
00572       for(i=0;i<nr+1;i++)
00573         {
00574           r[i] = Extrap[ie]->Pred_TotalBkgd->GetYaxis()->GetBinLowEdge(i+1);
00575         }
00576       for(i=0;i<np+1;i++)
00577         {
00578           p[i] = Extrap[ie]->Pred_TotalBkgd->GetXaxis()->GetBinLowEdge(i+1);
00579         }
00580       NueCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00581       NueCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00582       NuTauCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00583       NuTauCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00584     }
00585     
00586 
00587     //RBTNOTE: DO THIS SEPARATELY FOR RHC AND FHC!
00588     bool preExist = false;
00589     for (ie=0; ie<FidFiles.size(); ie++){
00590       if (!gSystem->AccessPathName(gSystem->ExpandPathName(FidFiles[ie].c_str()))) preExist = true;
00591     }
00592     for (ie=0; ie<CClikeFiles.size(); ie++){
00593       if (!gSystem->AccessPathName(gSystem->ExpandPathName(CClikeFiles[ie].c_str()))) preExist = true;
00594     }
00595 
00596     //If none of the files exists
00597     if (!preExist)
00598     {
00599       for(ie=0;ie<Extrap.size();ie++)
00600       {
00601 
00602         nt = Extrap[ie]->GetNTrue();
00603         t = new double[nt+1];
00604         for(i=0;i<nt+1;i++)
00605           {
00606             t[i] = Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNuTau]->GetZaxis()->GetBinLowEdge(i+1);
00607           }
00608         
00609         Pred_CC_Fid_Plus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Plus1",systname.c_str(),ie+1),"",nt,t));
00610         Pred_CC_Fid_Minus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Minus1",systname.c_str(),ie+1),"",nt,t));
00611       }
00612       
00613       return;
00614     }
00615     nfiles[0]=0;
00616     nfiles[1]=0;
00617     for(ie=0;ie<Extrap.size();ie++)
00618     {
00619       if(gSystem->AccessPathName(gSystem->ExpandPathName(CClikeFiles[ie].c_str())) && gSystem->AccessPathName(gSystem->ExpandPathName(FidFiles[ie].c_str())))//if files don't exist
00620       { 
00621         continue;
00622       }
00623      
00624       ietag = ExtrapType[ie];
00625 
00626       //
00627       nPOTFar=Extrap[ie]->GetFarPOT();
00628       nPOTNear=Extrap[ie]->GetNearPOT();
00629 
00630       fcc = new TFile(gSystem->ExpandPathName(CClikeFiles[ie].c_str()),"READ");
00631       treecc = (TTree*)fcc->Get("paramtree");
00632       treecc->SetBranchAddress("nearPOT",&npot_cc);
00633       treecc->SetBranchAddress("farPOT",&fpot_cc);
00634       treecc->GetEntry(0);
00635       
00636       ffid = new TFile(gSystem->ExpandPathName(FidFiles[ie].c_str()),"READ");
00637       treefid = (TTree*)ffid->Get("paramtree");
00638       treefid->SetBranchAddress("nearPOT",&npot_fid);
00639       treefid->SetBranchAddress("farPOT",&fpot_fid);
00640       treefid->GetEntry(0);
00641       
00642 
00643       data = (TH1D*)Extrap[ie]->NDData_Reco_CClike->Clone(Form("data_%i",ie));
00644       //should scale data! don't need to?
00645       if(nfiles[ietag]==0) totaldata[ietag] = (TH1D*)data->Clone();
00646       else totaldata[ietag]->Add(data);
00647       
00648 
00649       name = "NuMuCC_" + std + "_Fid/FD_TrueVsReco";
00650       farfidccstd = (TH2D*)ffid->Get(name.c_str());
00651       farfidccstd->Scale(nPOTFar/fpot_fid);
00652       if(nfiles[ietag]==0) totalfarfidccstd[ietag] = (TH2D*)farfidccstd->Clone();
00653       else totalfarfidccstd[ietag]->Add(farfidccstd);
00654       
00655       name = "NuMuCC_" + std + "_" + ccname + "/FD_TrueVsReco";
00656       farcclikestd = (TH2D*)fcc->Get(name.c_str());
00657       name = "NC_" + std + "_" + ccname + "/FD_TrueVsReco";
00658       h2 = (TH2D*)fcc->Get(name.c_str());
00659       farcclikestd->Add(h2);
00660       name = "BNueCC_" + std + "_" + ccname + "/FD_TrueVsReco";
00661       h2 = (TH2D*)fcc->Get(name.c_str());
00662       farcclikestd->Add(h2);
00663       farcclikestd->Scale(nPOTFar/fpot_cc);
00664       if(nfiles[ietag]==0) totalfarcclikestd[ietag] = (TH2D*)farcclikestd->Clone();
00665       else totalfarcclikestd[ietag]->Add(farcclikestd);
00666       
00667       name = "NuMuCC_" + std + "_" + ccname + "/ND_TrueVsReco";
00668       nearcclikestd = (TH2D*)fcc->Get(name.c_str());
00669       name = "NC_" + std + "_" + ccname + "/ND_TrueVsReco";
00670       h2 = (TH2D*)fcc->Get(name.c_str());
00671       nearcclikestd->Add(h2);
00672       name = "BNueCC_" + std + "_" + ccname + "/ND_TrueVsReco";
00673       h2 = (TH2D*)fcc->Get(name.c_str());
00674       nearcclikestd->Add(h2);
00675       nearcclikestd->Scale(nPOTNear/npot_cc);
00676       if(nfiles[ietag]==0) totalnearcclikestd[ietag] = (TH2D*)nearcclikestd->Clone();
00677       else totalnearcclikestd[ietag]->Add(nearcclikestd);
00678       
00679 
00680       name = "NuMuCC_" + plf + "_Fid/FD_TrueVsReco";
00681       if(plf==std)
00682       {
00683         farfidccpl = (TH2D*)farfidccstd->Clone(name.c_str());
00684       }
00685       else
00686       {
00687         farfidccpl = (TH2D*)ffid->Get(name.c_str());
00688         farfidccpl->Scale(nPOTFar/fpot_fid);
00689       }
00690       if(nfiles[ietag]==0) totalfarfidccpl[ietag] = (TH2D*)farfidccpl->Clone();
00691       else totalfarfidccpl[ietag]->Add(farfidccpl);
00692       
00693       name = "NuMuCC_" + plf + "_" + ccname + "/FD_TrueVsReco";
00694       if(plf==std)
00695       {
00696         farcclikepl = (TH2D*)farcclikestd->Clone(name.c_str());
00697       }
00698       else
00699       {
00700         farcclikepl = (TH2D*)fcc->Get(name.c_str());
00701         name = "NC_" + plf + "_" + ccname + "/FD_TrueVsReco";
00702         h2 = (TH2D*)fcc->Get(name.c_str());
00703         farcclikepl->Add(h2);
00704         name = "BNueCC_" + plf + "_" + ccname + "/FD_TrueVsReco";
00705         h2 = (TH2D*)fcc->Get(name.c_str());
00706         farcclikepl->Add(h2);
00707         farcclikepl->Scale(nPOTFar/fpot_cc);
00708       }
00709       if(nfiles[ietag]==0) totalfarcclikepl[ietag] = (TH2D*)farcclikepl->Clone();
00710       else totalfarcclikepl[ietag]->Add(farcclikepl);
00711       
00712       name = "NuMuCC_" + pln + "_" + ccname + "/ND_TrueVsReco";
00713       if(pln==std)
00714       {
00715         nearcclikepl = (TH2D*)nearcclikestd->Clone(name.c_str());
00716       }
00717       else
00718       {
00719         nearcclikepl = (TH2D*)fcc->Get(name.c_str());
00720         name = "NC_" + pln + "_" + ccname + "/ND_TrueVsReco";
00721         h2 = (TH2D*)fcc->Get(name.c_str());
00722         nearcclikepl->Add(h2);
00723         name = "BNueCC_" + pln + "_" + ccname + "/ND_TrueVsReco";
00724         h2 = (TH2D*)fcc->Get(name.c_str());
00725         nearcclikepl->Add(h2);
00726         nearcclikepl->Scale(nPOTNear/npot_cc);
00727       }
00728       if(nfiles[ietag]==0) totalnearcclikepl[ietag] = (TH2D*)nearcclikepl->Clone();
00729       else totalnearcclikepl[ietag]->Add(nearcclikepl);
00730      
00731       name = "NuMuCC_" + mnf + "_Fid/FD_TrueVsReco";
00732       if(mnf==std)
00733       {
00734         farfidccmn = (TH2D*)farfidccstd->Clone(name.c_str());
00735       }
00736       else if(mnf==plf)
00737       {
00738         farfidccmn = (TH2D*)farfidccpl->Clone(name.c_str());
00739       }
00740       else
00741       {
00742         farfidccmn = (TH2D*)ffid->Get(name.c_str());
00743         farfidccmn->Scale(nPOTFar/fpot_fid);
00744       }
00745       if(nfiles[ietag]==0) totalfarfidccmn[ietag] = (TH2D*)farfidccmn->Clone();
00746       else totalfarfidccmn[ietag]->Add(farfidccmn);
00747       
00748       name = "NuMuCC_" + mnf + "_" + ccname + "/FD_TrueVsReco";
00749       if(mnf==std)
00750       {
00751         farcclikemn = (TH2D*)farcclikestd->Clone(name.c_str());
00752       }
00753       else if(mnf==plf)
00754       {
00755         farcclikemn = (TH2D*)farcclikepl->Clone(name.c_str());
00756       }
00757       else
00758       {
00759         farcclikemn = (TH2D*)fcc->Get(name.c_str());
00760         name = "NC_" + mnf + "_" + ccname + "/FD_TrueVsReco";
00761         h2 = (TH2D*)fcc->Get(name.c_str());
00762         farcclikemn->Add(h2);
00763         name = "BNueCC_" + mnf + "_" + ccname + "/FD_TrueVsReco";
00764         h2 = (TH2D*)fcc->Get(name.c_str());
00765         farcclikemn->Add(h2);
00766         farcclikemn->Scale(nPOTFar/fpot_cc);
00767       }
00768       if(nfiles[ietag]==0) totalfarcclikemn[ietag] = (TH2D*)farcclikemn->Clone();
00769       else totalfarcclikemn[ietag]->Add(farcclikemn);
00770       
00771       name = "NuMuCC_" + mnn + "_" + ccname + "/ND_TrueVsReco";
00772       if(mnn==std)
00773       {
00774         nearcclikemn = (TH2D*)nearcclikestd->Clone(name.c_str());
00775       }
00776       else if(mnn==pln)
00777       {
00778         nearcclikemn = (TH2D*)nearcclikepl->Clone(name.c_str());
00779       }
00780       else
00781       {
00782         nearcclikemn = (TH2D*)fcc->Get(name.c_str());
00783         name = "NC_" + mnn + "_" + ccname + "/ND_TrueVsReco";
00784         h2 = (TH2D*)fcc->Get(name.c_str());
00785         nearcclikemn->Add(h2);
00786         name = "BNueCC_" + mnn + "_" + ccname + "/ND_TrueVsReco";
00787         h2 = (TH2D*)fcc->Get(name.c_str());
00788         nearcclikemn->Add(h2);
00789         nearcclikemn->Scale(nPOTNear/npot_cc);
00790       }
00791       if(nfiles[ietag]==0) totalnearcclikemn[ietag] = (TH2D*)nearcclikemn->Clone();
00792       else totalnearcclikemn[ietag]->Add(nearcclikemn);
00793       
00794       nfiles[ietag]++;
00795     }
00796 
00797     //Make both types of file
00798     for (ietag=0; ietag<2; ietag++){
00799       if (nfiles[ietag]>0){
00800       for (ie=0; ie<Extrap.size(); ie++){
00801         if (ExtrapType[ie]==ietag) {
00802           nt = Extrap[ie]->GetNTrue();
00803           t = new double[nt+1];
00804           for(i=0;i<nt+1;i++) t[i] = Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNuTau]->GetZaxis()->GetBinLowEdge(i+1);
00805           break;
00806         }
00807       }
00808 
00809       ccstd[ietag] = new TH1D(Form("ccstd%i",ietag),"",nt,t);
00810       for(it=0;it<nt;it++)
00811         {
00812           sum=0;
00813           for(ir=0;ir<totalfarcclikestd[ietag]->GetNbinsX();ir++)
00814             {
00815               if(totalnearcclikestd[ietag]->Integral(ir+1,ir+1,1,nt)>0)
00816                 {
00817                   sum += (totaldata[ietag]->GetBinContent(ir+1)/totalnearcclikestd[ietag]->Integral(ir+1,ir+1,1,nt))*totalfarcclikestd[ietag]->GetBinContent(ir+1,it+1);
00818                 }
00819             }
00820           temp=0;
00821           if(totalfarcclikestd[ietag]->Integral(1,totalfarcclikestd[ietag]->GetNbinsX(),it+1,it+1)>0)
00822             {
00823               temp = totalfarfidccstd[ietag]->Integral(1,totalfarfidccstd[ietag]->GetNbinsX(),it+1,it+1)*(sum/totalfarcclikestd[ietag]->Integral(1,totalfarcclikestd[ietag]->GetNbinsX(),it+1,it+1));
00824             }
00825           ccstd[ietag]->SetBinContent(it+1,temp);
00826         }
00827       
00828       ccpl[ietag] = new TH1D(Form("ccpl%i",ietag),"",nt,t);
00829       for(it=0;it<nt;it++)
00830         {
00831           sum=0;
00832           for(ir=0;ir<totalfarcclikepl[ietag]->GetNbinsX();ir++)
00833             {
00834               if(totalnearcclikepl[ietag]->Integral(ir+1,ir+1,1,nt)>0)
00835                 {
00836                   sum += (totaldata[ietag]->GetBinContent(ir+1)/totalnearcclikepl[ietag]->Integral(ir+1,ir+1,1,nt))*totalfarcclikepl[ietag]->GetBinContent(ir+1,it+1);
00837                 }
00838             }
00839           temp=0;
00840           if(totalfarcclikepl[ietag]->Integral(1,totalfarcclikepl[ietag]->GetNbinsX(),it+1,it+1)>0)
00841             {
00842               temp = totalfarfidccpl[ietag]->Integral(1,totalfarfidccpl[ietag]->GetNbinsX(),it+1,it+1)*(sum/totalfarcclikepl[ietag]->Integral(1,totalfarcclikepl[ietag]->GetNbinsX(),it+1,it+1));
00843             }
00844           ccpl[ietag]->SetBinContent(it+1,temp);
00845         }
00846       ccmn[ietag] = new TH1D(Form("ccmn%i",ietag),"",nt,t);
00847       for(it=0;it<nt;it++)
00848         {
00849           sum=0;
00850           for(ir=0;ir<totalfarcclikemn[ietag]->GetNbinsX();ir++)
00851             {
00852               if(totalnearcclikemn[ietag]->Integral(ir+1,ir+1,1,nt)>0)
00853                 {
00854                   sum += (totaldata[ietag]->GetBinContent(ir+1)/totalnearcclikemn[ietag]->Integral(ir+1,ir+1,1,nt))*totalfarcclikemn[ietag]->GetBinContent(ir+1,it+1);
00855                 }
00856             }
00857           temp=0;
00858           if(totalfarcclikemn[ietag]->Integral(1,totalfarcclikemn[ietag]->GetNbinsX(),it+1,it+1)>0)
00859             {
00860               temp = totalfarfidccmn[ietag]->Integral(1,totalfarfidccmn[ietag]->GetNbinsX(),it+1,it+1)*(sum/totalfarcclikemn[ietag]->Integral(1,totalfarcclikemn[ietag]->GetNbinsX(),it+1,it+1));
00861             }
00862           ccmn[ietag]->SetBinContent(it+1,temp);
00863         }
00864 
00865       ccpl[ietag]->Add(ccstd[ietag],-1.);
00866       ccpl[ietag]->Divide(ccstd[ietag]);
00867       ccmn[ietag]->Add(ccstd[ietag],-1.);
00868       ccmn[ietag]->Divide(ccstd[ietag]);
00869 
00870       if(plf==mnf && pln==mnn)
00871         {
00872           ccmn[ietag]->Scale(-1.);
00873         }
00874       }
00875     }
00876 
00877     for(ie=0;ie<Extrap.size();ie++)
00878       {
00879         ietag = ExtrapType[ie];
00880 
00881         Pred_CC_Fid_Plus1Sigma[systname].push_back((TH1D*)ccpl[ietag]->Clone(Form("%s_%i_Pred_CC_Fid_Plus1",systname.c_str(),ie+1)));
00882         Pred_CC_Fid_Minus1Sigma[systname].push_back((TH1D*)ccmn[ietag]->Clone(Form("%s_%i_Pred_CC_Fid_Minus1",systname.c_str(),ie+1)));
00883       }
00884     
00885 
00886       for(ie=0;ie<Extrap.size();ie++)
00887         {
00888           nr = Extrap[ie]->GetNReco();
00889           np = Extrap[ie]->GetNPID();
00890           nt = Extrap[ie]->GetNTrue();
00891 
00892           for(ip=0;ip<np;ip++)
00893             {
00894               for(ir=0;ir<nr;ir++)
00895                 {
00896                   sum_plus_tau=0;
00897                   sum_minus_tau=0;
00898                   sum_plus_nue=0;
00899                   sum_minus_nue=0;
00900                   for(it=0;it<nt;it++)
00901                     {
00902 
00903                       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));
00904                       
00905                       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));
00906                       
00907                       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));
00908                       
00909                       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));
00910                     }
00911 
00912                   //make fractional errors again
00913                   sum_plus_tau = sum_plus_tau/Extrap[ie]->Pred[Background::kNuTauCC]->GetBinContent(ip+1,ir+1);
00914                   sum_minus_tau = sum_minus_tau/Extrap[ie]->Pred[Background::kNuTauCC]->GetBinContent(ip+1,ir+1);
00915                   sum_plus_nue = sum_plus_nue/Extrap[ie]->Pred[Background::kNueCC]->GetBinContent(ip+1,ir+1);
00916                   sum_minus_nue = sum_minus_nue/Extrap[ie]->Pred[Background::kNueCC]->GetBinContent(ip+1,ir+1);
00917 
00918                   NuTauCC_MC_Plus1Sigma[systname][ie]->SetBinContent(ip+1,ir+1,sum_plus_tau);
00919                   NuTauCC_MC_Minus1Sigma[systname][ie]->SetBinContent(ip+1,ir+1,sum_minus_tau);
00920                   NueCC_MC_Plus1Sigma[systname][ie]->SetBinContent(ip+1,ir+1,sum_plus_nue);
00921                   NueCC_MC_Minus1Sigma[systname][ie]->SetBinContent(ip+1,ir+1,sum_minus_nue);
00922 
00923                 }
00924             }
00925         }
00926 
00927   }
00928   
00929   if(flag==2)
00930     {
00931 
00932     //create empty histograms for numuCC shift
00933     for(ie=0;ie<Extrap.size();ie++)
00934     {
00935       nt = Extrap[ie]->GetNTrue();
00936       t = new double[nt+1];
00937       for(i=0;i<nt+1;i++)
00938         {
00939           t[i] = Extrap[ie]->Pred_3D[Extrapolate2D::qNuMuToNuTau]->GetZaxis()->GetBinLowEdge(i+1);
00940         }
00941       Pred_CC_Fid_Plus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Plus1",systname.c_str(),ie+1),"",nt,t));
00942       Pred_CC_Fid_Minus1Sigma[systname].push_back(new TH1D(Form("%s_%i_Pred_CC_Fid_Minus1",systname.c_str(),ie+1),"",nt,t));
00943 
00944     }
00945 
00946     //RBTNOTE: DO THIS SEPARATELY FOR FHC AND RHC
00947     bool preExist = false;
00948 
00949     for (ie=0; ie<PreFiles.size(); ie++){
00950       if (!gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[ie].c_str()))) preExist = true;
00951     }
00952 
00953     if (!preExist)
00954     {
00955       for(ie=0;ie<Extrap.size();ie++)
00956       {
00957 
00958         nr = Extrap[ie]->GetNReco();
00959         np = Extrap[ie]->GetNPID();
00960 
00961         NueCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00962         NueCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NueCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00963         NuTauCC_MC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00964         NuTauCC_MC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_NuTauCC_MC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00965       }
00966       
00967       return;
00968     }
00969 
00970     for(j=0;j<app.size();j++)
00971     {
00972 
00973       nfiles[0]=0;
00974       nfiles[1]=0;
00975 
00976       for(ie=0;ie<Extrap.size();ie++)
00977       {
00978         if(gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[ie].c_str())))//if file doesn't exist
00979         {
00980           continue;
00981         }
00982 
00983         nr = Extrap[ie]->GetNReco();
00984         np = Extrap[ie]->GetNPID();
00985         nt = Extrap[ie]->GetNTrue();
00986         r = new double[nr+1];
00987         p = new double[np+1];
00988         for(i=0;i<nr+1;i++)
00989           {
00990             r[i] = Extrap[ie]->Pred_TotalBkgd->GetYaxis()->GetBinLowEdge(i+1);
00991           }
00992         for(i=0;i<np+1;i++)
00993           {
00994             p[i] = Extrap[ie]->Pred_TotalBkgd->GetXaxis()->GetBinLowEdge(i+1);
00995           }
00996         
00997         ietag=ExtrapType[ie];
00998 
00999         nPOTFar=Extrap[ie]->GetFarPOT();
01000         nPOTNear=Extrap[ie]->GetNearPOT();
01001 
01002         fpre = new TFile(gSystem->ExpandPathName(PreFiles[ie].c_str()),"READ");
01003         treepre = (TTree*)fpre->Get("paramtree");
01004         treepre->SetBranchAddress("nearPOT",&npot_pre);
01005         treepre->SetBranchAddress("farPOT",&fpot_pre);
01006         treepre->GetEntry(0);
01007 
01008         name = string(Background::AsString(app[j])) + "_" + std + "_Presel/FD_RecoVs" + Extrap[ie]->GetPID();
01009 
01010         prestd = (TH2D*)fpre->Get(name.c_str());
01011         prestd->Scale(nPOTFar/fpot_pre);
01012         MBH.Rebin2DHist(prestd,np,p,nr,r);
01013         if(nfiles[ietag]==0) totalprestd[ietag] = (TH2D*)prestd->Clone();
01014         else totalprestd[ietag]->Add(prestd);
01015 
01016         name = string(Background::AsString(app[j])) + "_" + plf + "_Presel/FD_RecoVs" + Extrap[ie]->GetPID();
01017         if(plf==std)
01018         {
01019           prepl = (TH2D*)prestd->Clone(name.c_str());
01020         }
01021         else
01022         {
01023           prepl = (TH2D*)fpre->Get(name.c_str());
01024           prepl->Scale(nPOTFar/fpot_pre);
01025           MBH.Rebin2DHist(prepl,np,p,nr,r);
01026         }
01027         if(nfiles[ietag]==0) totalprepl[ietag] = (TH2D*)prepl->Clone();
01028         else totalprepl[ietag]->Add(prepl);
01029         
01030         name = string(Background::AsString(app[j])) + "_" + mnf + "_Presel/FD_RecoVs" + Extrap[ie]->GetPID();
01031         if(mnf==std)
01032         {
01033           premn = (TH2D*)prestd->Clone(name.c_str());
01034         }
01035         else if(mnf==plf)
01036         {
01037           premn = (TH2D*)prepl->Clone(name.c_str());
01038         }
01039         else
01040         {
01041           premn = (TH2D*)fpre->Get(name.c_str());
01042           premn->Scale(nPOTFar/fpot_pre);
01043           MBH.Rebin2DHist(premn,np,p,nr,r);
01044         }
01045         if(nfiles[ietag]==0) totalpremn[ietag] = (TH2D*)premn->Clone();
01046         else totalpremn[ietag]->Add(premn);
01047         
01048         nfiles[ietag]++;
01049       }
01050       
01051 
01052       for (ietag=0; ietag<2; ietag++){
01053         if (nfiles[ietag]>0){
01054         totalprepl[ietag]->Add(totalprestd[ietag],-1.);
01055         totalprepl[ietag]->Divide(totalprestd[ietag]);
01056         totalpremn[ietag]->Add(totalprestd[ietag],-1.);
01057         totalpremn[ietag]->Divide(totalprestd[ietag]);
01058         
01059         if(plf==mnf && pln==mnn)
01060           {
01061             totalpremn[ietag]->Scale(-1.);
01062           }
01063         }
01064       }
01065 
01066       for(ie=0;ie<Extrap.size();ie++)
01067       {
01068         ietag=ExtrapType[ie];
01069 
01070         if(app[j]==Background::kNueCC)
01071         {
01072           NueCC_MC_Plus1Sigma[systname].push_back((TH2D*)totalprepl[ietag]->Clone(Form("%s_%i_NueCC_MC_Plus1",systname.c_str(),ie+1)));
01073           NueCC_MC_Minus1Sigma[systname].push_back((TH2D*)totalpremn[ietag]->Clone(Form("%s_%i_NueCC_MC_Minus1",systname.c_str(),ie+1)));
01074         }
01075         else if(app[j]==Background::kNuTauCC)
01076         {
01077           NuTauCC_MC_Plus1Sigma[systname].push_back((TH2D*)totalprepl[ietag]->Clone(Form("%s_%i_NuTauCC_MC_Plus1",systname.c_str(),ie+1)));
01078           NuTauCC_MC_Minus1Sigma[systname].push_back((TH2D*)totalpremn[ietag]->Clone(Form("%s_%i_NuTauCC_MC_Minus1",systname.c_str(),ie+1)));
01079         }
01080       }
01081     }
01082   }
01083   
01084   return;
01085 }

void ErrorCalc_Joint::ReadSysFiles_FNExtrap_AllRuns ( int  n  )  [private, virtual]

---

Reimplemented from ErrorCalc.

Definition at line 95 of file ErrorCalc_Joint.cxx.

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

00096 {
00097 
00098   int i;
00099   
00100   string plf,pln,mnf,mnn,std,systname,filetag;
00101   systname = FNExtrap_SystName.at(n);
00102   plf = FNExtrap_FarPlusTag.at(n);
00103   pln = FNExtrap_NearPlusTag.at(n);
00104   mnf = FNExtrap_FarMinusTag.at(n);
00105   mnn = FNExtrap_NearMinusTag.at(n);
00106   std = FNExtrap_StdTag.at(n);
00107   filetag=FNExtrap_FileTag.at(n);
00108   int flag = FNExtrap_Flag.at(n);
00109   
00110   unsigned int ie;
00111 
00112   vector<string> PreFiles;
00113   
00114   for(ie=0;ie<Extrap.size();ie++)
00115   {
00116     if (ExtrapType[ie]==0) PreFiles.push_back(SysFileDir+"/SysFile2D_FHC"+filetag+"_Pre_1.root");
00117     if (ExtrapType[ie]==1) PreFiles.push_back(SysFileDir+"/SysFile2D_RHC"+filetag+"_Pre_1.root");
00118   }
00119 
00120   unsigned int j;
00121   vector<Background::Background_t> bgs;
00122   bgs.push_back(Background::kNC);
00123   bgs.push_back(Background::kNuMuCC);
00124   bgs.push_back(Background::kBNueCC);
00125   
00126   string name;
00127   int np,nr;
00128   double *r,*p;
00129   TFile *fpre;
00130   TH2D *farstd=0,*farpl=0,*farmn=0;
00131   TH2D *nearstd=0,*nearpl=0,*nearmn=0;
00132   TH2D *totalfarstd[2]={0},*totalfarpl[2]={0},*totalfarmn[2]={0};
00133   TH2D *totalnearstd[2]={0},*totalnearpl[2]={0},*totalnearmn[2]={0};
00134   TH2D *fnpl[2]={0}, *fnmn[2]={0};
00135   double npot_pre,fpot_pre;
00136   TTree *treepre;
00137   double nPOTNear,nPOTFar;
00138   int nfiles[2]={0};
00139    
00140   //RBTNOTE: THIS SHOULD BE DONE SEPARATELY FOR RHC AND FHC:
00141   //Addendum: Why not just do it separately for each extrapolation?
00142   bool preExist = false;
00143   for (ie=0; ie<PreFiles.size(); ie++){
00144     if (!gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[ie].c_str()))) preExist = true;
00145   }
00146 
00147   //flag==1 means systematics related to numuCC extrapolation for signal/tau prediction - by definition 0 for these components
00148   if(flag==1 || (!preExist))//if all three Pre files do NOT exist
00149   {
00150     for(ie=0;ie<Extrap.size();ie++)
00151     {
00152 
00153       nr = Extrap[ie]->GetNReco();
00154       np = Extrap[ie]->GetNPID();
00155       r = new double[nr+1];
00156       p = new double[np+1];
00157       for(i=0;i<nr+1;i++)
00158         {
00159           r[i] = Extrap[ie]->Pred_TotalBkgd->GetYaxis()->GetBinLowEdge(i+1);
00160         }
00161       for(i=0;i<np+1;i++)
00162         {
00163           p[i] = Extrap[ie]->Pred_TotalBkgd->GetXaxis()->GetBinLowEdge(i+1);
00164         }
00165 
00166       FN_NC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00167       FN_NC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00168       
00169       FN_NuMuCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NuMuCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00170       FN_NuMuCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_NuMuCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00171       
00172       FN_BNueCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_BNueCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00173       FN_BNueCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FN_BNueCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00174       
00175       ND_NC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00176       ND_NC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00177       
00178       ND_NuMuCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NuMuCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00179       ND_NuMuCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_NuMuCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00180       
00181       ND_BNueCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_BNueCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00182       ND_BNueCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_ND_BNueCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00183       
00184       FD_NC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00185       FD_NC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00186       
00187       FD_NuMuCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NuMuCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00188       FD_NuMuCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_NuMuCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00189       
00190       FD_BNueCC_Plus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_BNueCC_Plus1",systname.c_str(),ie+1),"",np,p,nr,r));
00191       FD_BNueCC_Minus1Sigma[systname].push_back(new TH2D(Form("%s_%i_FD_BNueCC_Minus1",systname.c_str(),ie+1),"",np,p,nr,r));
00192     }
00193     
00194     return;
00195   }
00196 
00197   int ietag=-1;
00198   
00199   for(j=0;j<bgs.size();j++)
00200   {
00201 
00202     //Set number of FHC and RHC files to 0
00203     nfiles[0]=0;
00204     nfiles[1]=0;
00205 
00206     for(ie=0;ie<Extrap.size();ie++)
00207     {
00208       if(gSystem->AccessPathName(gSystem->ExpandPathName(PreFiles[ie].c_str())))//if file doesn't exist
00209       {
00210         continue;
00211       }
00212 
00213       //Get a tag of 0 for FHC, 1 for RHC
00214       ietag = ExtrapType[ie];
00215 
00216       nr = Extrap[ie]->GetNReco();
00217       np = Extrap[ie]->GetNPID();
00218       r = new double[nr+1];
00219       p = new double[np+1];
00220       for(i=0;i<nr+1;i++)
00221         {
00222           r[i] = Extrap[ie]->Pred_TotalBkgd->GetYaxis()->GetBinLowEdge(i+1);
00223         }
00224       for(i=0;i<np+1;i++)
00225         {
00226           p[i] = Extrap[ie]->Pred_TotalBkgd->GetXaxis()->GetBinLowEdge(i+1);
00227         }
00228 
00229       //Scaling:
00230       nPOTFar=Extrap[ie]->GetFarPOT();
00231       nPOTNear=Extrap[ie]->GetNearPOT();
00233 
00234       fpre = new TFile(gSystem->ExpandPathName(PreFiles[ie].c_str()),"READ");
00235       
00236       treepre = (TTree*)fpre->Get("paramtree");
00237       treepre->SetBranchAddress("nearPOT",&npot_pre);
00238       treepre->SetBranchAddress("farPOT",&fpot_pre);
00239       treepre->GetEntry(0);
00240       
00241       name = string(Background::AsString(bgs[j])) + "_" + std + "_Presel/FD_RecoVs" + Extrap[ie]->GetPID();
00242       farstd = (TH2D*)fpre->Get(name.c_str());
00243       farstd->Scale(nPOTFar/fpot_pre);
00244       MBH.Rebin2DHist(farstd,np,p,nr,r);
00245       if(nfiles[ietag]==0) totalfarstd[ietag] = (TH2D*)farstd->Clone();
00246       else totalfarstd[ietag]->Add(farstd);
00247       
00248       name = string(Background::AsString(bgs[j])) + "_" + plf + "_Presel/FD_RecoVs" + Extrap[ie]->GetPID();
00249       if(plf==std)
00250       {
00251         farpl = (TH2D*)farstd->Clone(name.c_str());
00252       }
00253       else
00254       {
00255         farpl = (TH2D*)fpre->Get(name.c_str());
00256         farpl->Scale(nPOTFar/fpot_pre);
00257         MBH.Rebin2DHist(farpl,np,p,nr,r);
00258       }
00259       if(nfiles[ietag]==0) totalfarpl[ietag] = (TH2D*)farpl->Clone();
00260       else totalfarpl[ietag]->Add(farpl);
00261       
00262       name = string(Background::AsString(bgs[j])) + "_" + mnf + "_Presel/FD_RecoVs" + Extrap[ie]->GetPID();
00263       if(mnf==std)
00264       {
00265         farmn = (TH2D*)farstd->Clone(name.c_str());
00266       }
00267       else if(mnf==plf)
00268       {
00269         farmn = (TH2D*)farpl->Clone(name.c_str());
00270       }
00271       else
00272       {
00273         farmn = (TH2D*)fpre->Get(name.c_str());
00274         farmn->Scale(nPOTFar/fpot_pre);
00275         MBH.Rebin2DHist(farmn,np,p,nr,r);
00276       }
00277       if(nfiles[ietag]==0) totalfarmn[ietag] = (TH2D*)farmn->Clone();
00278       else totalfarmn[ietag]->Add(farmn);
00279       
00280       name = string(Background::AsString(bgs[j])) + "_" + std + "_Presel/ND_RecoVs" + Extrap[ie]->GetPID();
00281       nearstd = (TH2D*)fpre->Get(name.c_str());
00282       nearstd->Scale(nPOTNear/npot_pre);
00283       MBH.Rebin2DHist(nearstd,np,p,nr,r);
00284       if(nfiles[ietag]==0) totalnearstd[ietag] = (TH2D*)nearstd->Clone();
00285       else totalnearstd[ietag]->Add(nearstd);
00286       
00287       name = string(Background::AsString(bgs[j])) + "_" + pln + "_Presel/ND_RecoVs" + Extrap[ie]->GetPID();
00288       if(pln==std)
00289       {
00290         nearpl = (TH2D*)nearstd->Clone(name.c_str());
00291       }
00292       else
00293       {
00294         nearpl = (TH2D*)fpre->Get(name.c_str());
00295         nearpl->Scale(nPOTNear/npot_pre);
00296         MBH.Rebin2DHist(nearpl,np,p,nr,r);
00297       }
00298       if(nfiles[ietag]==0) totalnearpl[ietag] = (TH2D*)nearpl->Clone();
00299       else totalnearpl[ietag]->Add(nearpl);
00300       
00301       name = string(Background::AsString(bgs[j])) + "_" + mnn + "_Presel/ND_RecoVs" + Extrap[ie]->GetPID();
00302       if(mnn==std)
00303       {
00304         nearmn=(TH2D*)nearstd->Clone(name.c_str());
00305       }
00306       else if(mnn==pln)
00307       {
00308         nearmn=(TH2D*)nearpl->Clone(name.c_str());
00309       }
00310       else
00311       {
00312         nearmn = (TH2D*)fpre->Get(name.c_str());
00313         nearmn->Scale(nPOTNear/npot_pre);
00314         MBH.Rebin2DHist(nearmn,np,p,nr,r);
00315       }
00316       if(nfiles[ietag]==0) totalnearmn[ietag] = (TH2D*)nearmn->Clone();
00317       else totalnearmn[ietag]->Add(nearmn);
00318       
00319       nfiles[ietag]++;
00320     }
00321     
00322     for(ie=0;ie<Extrap.size();ie++)
00323     {
00324 
00325       //Get a tag of 0 for FHC, 1 for RHC
00326       ietag = ExtrapType[ie];
00327 
00328       if(bgs[j]==Background::kNC)
00329       {
00330         ND_NC_Plus1Sigma[systname].push_back((TH2D*)totalnearpl[ietag]->Clone(Form("%s_ND_NC_Plus1",systname.c_str())));
00331         ND_NC_Plus1Sigma[systname][ie]->Add(totalnearstd[ietag],-1.0);
00332         ND_NC_Plus1Sigma[systname][ie]->Divide(totalnearstd[ietag]);
00333         ND_NC_Minus1Sigma[systname].push_back((TH2D*)totalnearmn[ietag]->Clone(Form("%s_ND_NC_Minus1",systname.c_str())));
00334         ND_NC_Minus1Sigma[systname][ie]->Add(totalnearstd[ietag],-1.0);
00335         ND_NC_Minus1Sigma[systname][ie]->Divide(totalnearstd[ietag]);
00336         if(plf==mnf && pln==mnn)
00337         {
00338           ND_NC_Minus1Sigma[systname][ie]->Scale(-1.);
00339         }
00340         
00341         FD_NC_Plus1Sigma[systname].push_back((TH2D*)totalfarpl[ietag]->Clone(Form("%s_FD_NC_Plus1",systname.c_str())));
00342         FD_NC_Plus1Sigma[systname][ie]->Add(totalfarstd[ietag],-1.0);
00343         FD_NC_Plus1Sigma[systname][ie]->Divide(totalfarstd[ietag]);
00344         FD_NC_Minus1Sigma[systname].push_back((TH2D*)totalfarmn[ietag]->Clone(Form("%s_FD_NC_Minus1",systname.c_str())));
00345         FD_NC_Minus1Sigma[systname][ie]->Add(totalfarstd[ietag],-1.0);
00346         FD_NC_Minus1Sigma[systname][ie]->Divide(totalfarstd[ietag]);
00347         if(plf==mnf && pln==mnn)
00348         {
00349           FD_NC_Minus1Sigma[systname][ie]->Scale(-1.);
00350         }
00351       }
00352       else if(bgs[j]==Background::kNuMuCC)
00353       {
00354         ND_NuMuCC_Plus1Sigma[systname].push_back((TH2D*)totalnearpl[ietag]->Clone(Form("%s_ND_NuMuCC_Plus1",systname.c_str())));
00355         ND_NuMuCC_Plus1Sigma[systname][ie]->Add(totalnearstd[ietag],-1.0);
00356         ND_NuMuCC_Plus1Sigma[systname][ie]->Divide(totalnearstd[ietag]);
00357         ND_NuMuCC_Minus1Sigma[systname].push_back((TH2D*)totalnearmn[ietag]->Clone(Form("%s_ND_NuMuCC_Minus1",systname.c_str())));
00358         ND_NuMuCC_Minus1Sigma[systname][ie]->Add(totalnearstd[ietag],-1.0);
00359         ND_NuMuCC_Minus1Sigma[systname][ie]->Divide(totalnearstd[ietag]);
00360         if(plf==mnf && pln==mnn)
00361         {
00362           ND_NuMuCC_Minus1Sigma[systname][ie]->Scale(-1.);
00363         }
00364         
00365         FD_NuMuCC_Plus1Sigma[systname].push_back((TH2D*)totalfarpl[ietag]->Clone(Form("%s_FD_NuMuCC_Plus1",systname.c_str())));
00366         FD_NuMuCC_Plus1Sigma[systname][ie]->Add(totalfarstd[ietag],-1.0);
00367         FD_NuMuCC_Plus1Sigma[systname][ie]->Divide(totalfarstd[ietag]);
00368         FD_NuMuCC_Minus1Sigma[systname].push_back((TH2D*)totalfarmn[ietag]->Clone(Form("%s_FD_NuMuCC_Minus1",systname.c_str())));
00369         FD_NuMuCC_Minus1Sigma[systname][ie]->Add(totalfarstd[ietag],-1.0);
00370         FD_NuMuCC_Minus1Sigma[systname][ie]->Divide(totalfarstd[ietag]);
00371         if(plf==mnf && pln==mnn)
00372         {
00373           FD_NuMuCC_Minus1Sigma[systname][ie]->Scale(-1.);
00374         }
00375       }
00376       else if(bgs[j]==Background::kBNueCC)
00377       {
00378         ND_BNueCC_Plus1Sigma[systname].push_back((TH2D*)totalnearpl[ietag]->Clone(Form("%s_ND_BNueCC_Plus1",systname.c_str())));
00379         ND_BNueCC_Plus1Sigma[systname][ie]->Add(totalnearstd[ietag],-1.0);
00380         ND_BNueCC_Plus1Sigma[systname][ie]->Divide(totalnearstd[ietag]);
00381         ND_BNueCC_Minus1Sigma[systname].push_back((TH2D*)totalnearmn[ietag]->Clone(Form("%s_ND_BNueCC_Minus1",systname.c_str())));
00382         ND_BNueCC_Minus1Sigma[systname][ie]->Add(totalnearstd[ietag],-1.0);
00383         ND_BNueCC_Minus1Sigma[systname][ie]->Divide(totalnearstd[ietag]);
00384         if(plf==mnf && pln==mnn)
00385         {
00386           ND_BNueCC_Minus1Sigma[systname][ie]->Scale(-1.);
00387         }
00388         
00389         FD_BNueCC_Plus1Sigma[systname].push_back((TH2D*)totalfarpl[ietag]->Clone(Form("%s_FD_BNueCC_Plus1",systname.c_str())));
00390         FD_BNueCC_Plus1Sigma[systname][ie]->Add(totalfarstd[ietag],-1.0);
00391         FD_BNueCC_Plus1Sigma[systname][ie]->Divide(totalfarstd[ietag]);
00392         FD_BNueCC_Minus1Sigma[systname].push_back((TH2D*)totalfarmn[ietag]->Clone(Form("%s_FD_BNueCC_Minus1",systname.c_str())));
00393         FD_BNueCC_Minus1Sigma[systname][ie]->Add(totalfarstd[ietag],-1.0);
00394         FD_BNueCC_Minus1Sigma[systname][ie]->Divide(totalfarstd[ietag]);
00395         if(plf==mnf && pln==mnn)
00396         {
00397           FD_BNueCC_Minus1Sigma[systname][ie]->Scale(-1.);
00398         }
00399       }
00400     }
00401     
00402     //Make the total histograms for both FHC(itype=0) and RHC(itype=1)
00403     for (int itype=0; itype<2; itype++){
00404       if (nfiles[itype]>0){
00405         if (totalfarstd[itype]) totalfarstd[itype]->Divide(totalnearstd[itype]);
00406         if (totalfarpl[itype]) totalfarpl[itype]->Divide(totalnearpl[itype]);
00407         if (totalfarmn[itype]) totalfarmn[itype]->Divide(totalnearmn[itype]);
00408         
00409         fnpl[itype] = (TH2D*)totalfarpl[itype]->Clone("fnpl");
00410         fnpl[itype]->Add(totalfarstd[itype],-1.);
00411         fnpl[itype]->Divide(totalfarstd[itype]);
00412         
00413         fnmn[itype] = (TH2D*)totalfarmn[itype]->Clone("fnmn");
00414         fnmn[itype]->Add(totalfarstd[itype],-1.);
00415         fnmn[itype]->Divide(totalfarstd[itype]);
00416         
00417         if(plf==mnf && pln==mnn)
00418           {
00419             fnmn[itype]->Scale(-1.);
00420           }
00421       }
00422     }
00423 
00424     for(ie=0;ie<Extrap.size();ie++)
00425     {
00426 
00427       //Get a tag of 0 for FHC, 1 for RHC
00428       ietag = ExtrapType[ie];
00429 
00430       if(bgs[j]==Background::kNC)
00431       {
00432         FN_NC_Plus1Sigma[systname].push_back((TH2D*)fnpl[ietag]->Clone(Form("%s_FN_NC_Plus1",systname.c_str())));
00433         FN_NC_Minus1Sigma[systname].push_back((TH2D*)fnmn[ietag]->Clone(Form("%s_FN_NC_Minus1",systname.c_str())));
00434       }
00435       if(bgs[j]==Background::kNuMuCC)
00436       {
00437         FN_NuMuCC_Plus1Sigma[systname].push_back((TH2D*)fnpl[ietag]->Clone(Form("%s_FN_NuMuCC_Plus1",systname.c_str())));
00438         FN_NuMuCC_Minus1Sigma[systname].push_back((TH2D*)fnmn[ietag]->Clone(Form("%s_FN_NuMuCC_Minus1",systname.c_str())));
00439       }
00440       if(bgs[j]==Background::kBNueCC)
00441       {
00442         FN_BNueCC_Plus1Sigma[systname].push_back((TH2D*)fnpl[ietag]->Clone(Form("%s_FN_BNueCC_Plus1",systname.c_str())));
00443         FN_BNueCC_Minus1Sigma[systname].push_back((TH2D*)fnmn[ietag]->Clone(Form("%s_FN_BNueCC_Minus1",systname.c_str())));
00444       }
00445     }
00446   }
00447   
00448   return;
00449 }

void ErrorCalc_Joint::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 from ErrorCalc.

Definition at line 1864 of file ErrorCalc_Joint.cxx.

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

01865 {
01866   if(Extrap.size()==0)
01867   {
01868     cout<<"Error in ErrorCalc::SetGridPred(): You must add an Extrapolate2D object.  Quitting..."<<endl;
01869     return;
01870   }
01871 
01872   unsigned int ie;
01873   int nbinsFHC=0;
01874   int nbinsRHC=0;
01875 
01876   for (ie=0;ie<Extrap.size();ie++){
01877     if (ExtrapType[ie]==0) nbinsFHC=GridPred[Background::kNC][ie]->GetNbinsX();
01878     if (ExtrapType[ie]==1) nbinsRHC=GridPred[Background::kNC][ie]->GetNbinsX();
01879   }
01880 
01881   if (nbins != (nbinsFHC + nbinsRHC)) cout << "Warning!  SetGridPred mismatch!" << endl;
01882 
01883   Initialize();
01884   if(!Init) return;
01885   
01886   int i;
01887   int ii;  
01888 
01889   unsigned int ieFHC=0;
01890   unsigned int ieRHC=0;
01891   unsigned int ietmp=0;
01892 
01893   for(ie=0;ie<Extrap.size();ie++)
01894   {
01895     GridPred[Background::kNC][ie]->Reset();
01896     GridPred[Background::kNuMuCC][ie]->Reset();
01897     GridPred[Background::kBNueCC][ie]->Reset();
01898     GridPred[Background::kNuTauCC][ie]->Reset();
01899     GridPred[Background::kNueCC][ie]->Reset();
01900    
01901     for(i=0;i<nbins;i++)
01902     {
01903       if (ExtrapType[ie]==0) { ii=i; ietmp = ieFHC; }
01904       if (ExtrapType[ie]==1) { ii=i-nbinsFHC; ietmp = ieRHC; }
01905 
01906       if ((ExtrapType[ie]==0&&i<nbinsFHC)||
01907           (ExtrapType[ie]==1&&i>=nbinsFHC)){
01908 
01909         GridPred[Background::kNC][ie]->SetBinContent(ii+1,nc[i][ietmp]);
01910         GridPred[Background::kNuMuCC][ie]->SetBinContent(ii+1,cc[i][ietmp]);
01911         GridPred[Background::kBNueCC][ie]->SetBinContent(ii+1,bnue[i][ietmp]);
01912         GridPred[Background::kNuTauCC][ie]->SetBinContent(ii+1,tau[i][ietmp]);
01913         GridPred[Background::kNueCC][ie]->SetBinContent(ii+1,sig[i][ietmp]);
01914       }
01915     }
01916 
01917     if (ExtrapType[ie]==0) ieFHC++;
01918     if (ExtrapType[ie]==1) ieRHC++;
01919 
01920   }
01921   
01922   return;
01923 }


Member Data Documentation

std::map<string,TH2D*> ErrorCalc_Joint::ExtraCovariance [private]

Reimplemented from ErrorCalc.

Definition at line 51 of file ErrorCalc_Joint.h.

Referenced by CalculateSystErrorMatrixExtrap(), CalculateSystErrorMatrixGrid(), and ReadCovarianceFiles().

Definition at line 53 of file ErrorCalc_Joint.h.

Referenced by AddCovarianceMatrix(), and ReadCovarianceFiles().

std::map<string,int> ErrorCalc_Joint::ExtraCovarianceTag [private]

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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1