NueFit2D_Joint Class Reference

#include <NueFit2D_Joint.h>

Inheritance diagram for NueFit2D_Joint:
NueFit2D

List of all members.

Public Member Functions

 NueFit2D_Joint ()
virtual ~NueFit2D_Joint ()
void ScalePredictionFHC (double scale=1.0)
void ScalePredictionRHC (double scale=1.0)
void RunMultiBinPseudoExpts (bool Print=true)
void ReadGridFiles ()
void RunDataGrid (string filenorm, string fileinvt)
double GetSensitivityAt (double delta=0, bool normalhier=true)
void RunDeltaChi2Contour (int cl=0)
double GetLikelihood (double t12=0.6, double t23=0.785398, double t13=0, double dm2_32=2.32e-3, double dm2_21=7.59e-5, double delta=0)
void DoDeltaFit ()
void AddExtrap (Extrapolate2D *E)
void AddExtrapRHC (Extrapolate2D *E)
void AddExtrapFHC (Extrapolate2D *E)
void CombineFHCRHC ()
void RunMultiBinPseudoExpts_MHDeltaFit (bool Print=true)
void RunMultiBinFC_MHDeltaFit ()
void DoNSIFitQuick ()
void DoNSIFitForever ()
void GetNSIFFX2 ()
void RunNSIDeltaChi2Contour (int cl=0)
void DoIH (int doihinput=0)
void SetPhaseBroom (double pb=0.0)

Private Member Functions

void CalculateDlnLMatrix (TH2D *SystMatrix, TH2D *HOOHEMatrix)
void DefineStdDlnLMinuit ()
void DefineBinDlnLMinuit ()
double GetMinLikelihood (double delta=0, bool normalhier=true)
double GetMinLikelihood_Delta (bool normalhier=true)

Private Attributes

unsigned int nBinsFHC
unsigned int nBinsRHC
double PhaseBroom
int NHIH
double GridScale_Normal_FHC
double GridScale_Inverted_FHC
double GridScale_Normal_RHC
double GridScale_Inverted_RHC
bool printMatrix
double potScaleFHC
double potScaleRHC
vector< Extrapolate2D * > ExtrapRHC
vector< Extrapolate2D * > ExtrapFHC

Detailed Description

Definition at line 28 of file NueFit2D_Joint.h.


Constructor & Destructor Documentation

NueFit2D_Joint::NueFit2D_Joint (  ) 

Definition at line 6 of file NueFit2D_Joint.cxx.

References NueFit2D::delta, NueFit2D::Delta_emu, NueFit2D::Delta_etau, NueFit2D::Delta_etau_unc, NueFit2D::Delta_mutau, NueFit2D::DeltaMSq21, NueFit2D::DeltaMSq21Unc, NueFit2D::DeltaMSq32, NueFit2D::DeltaMSq32Unc, NueFit2D::Eps_ee, NueFit2D::Eps_emu, NueFit2D::Eps_etau, NueFit2D::Eps_etau_unc, NueFit2D::Eps_mumu, NueFit2D::Eps_mutau, NueFit2D::Eps_tautau, NueFit2D::ErrCalc, NueFit2D::ErrorMatrix, NueFit2D::ExternalErrorMatrix, NueFit2D::FormatChi2Hists(), NueFit2D::GridScale_Inverted, GridScale_Inverted_FHC, GridScale_Inverted_RHC, NueFit2D::GridScale_Normal, GridScale_Normal_FHC, GridScale_Normal_RHC, NueFit2D::IncludeOscParErrs, NueFit2D::InvErrorMatrix, NueFit2D::nBins, nBinsFHC, nBinsRHC, NueFit2D::NObs, potScaleFHC, potScaleRHC, printMatrix, NueFit2D::SetDeltaRange(), NueFit2D::SetepsetauRange(), NueFit2D::SetFitMethod(), NueFit2D::SetGridNorm(), NueFit2D::SetNDeltaSteps(), NueFit2D::SetNepsetauSteps(), NueFit2D::SetNExperiments(), NueFit2D::SetNSinSq2Th13Steps(), NueFit2D::SetOutputFile(), NueFit2D::SetSinSq2Th13Range(), NueFit2D::Theta12, NueFit2D::Theta12Unc, NueFit2D::Theta13, NueFit2D::Theta13Unc, NueFit2D::Theta23, and NueFit2D::Theta23Unc.

00007 {
00008 
00009   printMatrix=false;
00010 
00011   NObs=0;
00012   ErrCalc=0;
00013   ErrorMatrix=0;
00014   InvErrorMatrix=0;
00015   ExternalErrorMatrix=0;
00016 
00017   SetOutputFile();
00018   SetNExperiments();
00019   
00020   SetGridNorm();
00021   GridScale_Normal=1.;
00022   GridScale_Inverted=1.;
00023   GridScale_Normal_FHC=1.;
00024   GridScale_Inverted_FHC=1.;
00025   GridScale_Normal_RHC=1.;
00026   GridScale_Inverted_RHC=1.;
00027   
00028   nBins=0;
00029   nBinsFHC=0;
00030   nBinsRHC=0;
00031   
00032   IncludeOscParErrs = false;
00033   
00034   FormatChi2Hists();
00035   
00036   SetNDeltaSteps();
00037   SetNSinSq2Th13Steps();
00038   SetDeltaRange();
00039   SetSinSq2Th13Range();
00040   SetepsetauRange();
00041   SetNepsetauSteps();
00042   
00043   SetFitMethod();
00044   
00045   potScaleFHC=1.0;
00046   potScaleRHC=1.0;
00047 
00048   Theta12 = 0.59365;
00049   Theta12Unc = 0.041;
00050   Theta23 = 0.7854;
00051   Theta23Unc = 0.122;
00052   Theta13 = 0.1594;
00053   Theta13Unc = 0.011;
00054   DeltaMSq32 = 2.32e-3;
00055   DeltaMSq32Unc = 0.11e-3;
00056   DeltaMSq21 = 8.0e-5;
00057   DeltaMSq21Unc = 0.6e-5;
00058   delta = 0;
00059   
00060   //NSI parameters:
00061   Eps_ee = 0.0;
00062   Eps_emu = 0.0;
00063   Eps_etau = 0.0;
00064   Eps_etau_unc = 0.01;
00065   Eps_mumu = 0.0;
00066   Eps_mutau = 0.0; 
00067   Eps_tautau = 0.0; 
00068   Delta_emu = 0.0;
00069   Delta_etau = 0.0;
00070   Delta_etau_unc = 0.01;
00071   Delta_mutau = 0.0;
00072 
00073   return;
00074 }

NueFit2D_Joint::~NueFit2D_Joint (  )  [virtual]

Definition at line 75 of file NueFit2D_Joint.cxx.

00076 {
00077 }


Member Function Documentation

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

Reimplemented from NueFit2D.

Definition at line 174 of file NueFit2D_Joint.cxx.

00175 {
00176 
00177   cout << "Can't add generic extrapolation for a joint fit.  Use FHC/RHC options instead.  No extrapolation was added." << endl;
00178   return;
00179 
00180 }

void NueFit2D_Joint::AddExtrapFHC ( Extrapolate2D E  ) 

Definition at line 182 of file NueFit2D_Joint.cxx.

References ExtrapFHC.

00183 {
00184   ExtrapFHC.push_back(E);
00185   return;
00186 
00187 }

void NueFit2D_Joint::AddExtrapRHC ( Extrapolate2D E  ) 

Definition at line 189 of file NueFit2D_Joint.cxx.

References ExtrapRHC.

00190 {
00191   ExtrapRHC.push_back(E);
00192   return;
00193 
00194 }

void NueFit2D_Joint::CalculateDlnLMatrix ( TH2D *  SystMatrix = 0,
TH2D *  HOOHEMatrix = 0 
) [private, virtual]

Reimplemented from NueFit2D.

Definition at line 80 of file NueFit2D_Joint.cxx.

References ExtrapFHC, ExtrapRHC, NueFit2D::InvErrorMatrix, NueFit2D::nBins, potScaleFHC, potScaleRHC, and printMatrix.

00080                                                                               {
00081   //Calculate the covariance matrix for the "bin by bin" method
00082 
00083   //Make new inverted matrix if not done already
00084   if (InvErrorMatrix==0){
00085     if (SysMatrix){
00086       InvErrorMatrix = (TH2D*)SysMatrix->Clone("InvErrorMatrix");
00087     } else if (HOOHEMatrix) {
00088       InvErrorMatrix = (TH2D*)HOOHEMatrix->Clone("InvErrorMatrix");
00089     } else {
00090       cout << "Didn't give me any matrices to invert!" << endl;
00091       return;
00092     }
00093   }
00094   //Reset everything to zero:
00095   InvErrorMatrix->Reset();
00096   int totbins = nBins;
00097   int i=0;
00098   int j=0;
00099   int k=0;
00100   double scalei=1.0;
00101   double scalej=1.0;
00102   int nFHC=0;
00103   int nRHC=0;
00104 
00105   if (ExtrapFHC.size()>0) nFHC = ExtrapFHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
00106   if (ExtrapRHC.size()>0) nRHC = ExtrapRHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
00107   if (printMatrix) cout << "PRE-INVERSION" << endl;
00108   //Add together error elements into a single array
00109   double *Varray = new double[totbins*totbins];
00110   for(i=0;i<totbins;i++)
00111   {
00112 
00113     if (i<nFHC) { scalei=potScaleFHC; }
00114     else { scalei=potScaleRHC; }
00115 
00116     for(j=0;j<totbins;j++)
00117     {
00118 
00119     if (j<nFHC) { scalej=potScaleFHC; }
00120     else { scalej=potScaleRHC; }
00121 
00122       k = i*totbins + j;
00123       Varray[k]=0;
00124       if (SysMatrix!=0) {
00125         Varray[k] += scalei*scalej*SysMatrix->GetBinContent(i+1,j+1);
00126       }
00127 
00128       if (HOOHEMatrix!=0) {
00129         Varray[k] += HOOHEMatrix->GetBinContent(i+1,j+1);
00130       }
00131 
00132       Varray[k] = Varray[k]*1.0e0;
00133 
00134     }
00135 
00136  }
00137 
00138   //Hand elements to a TMatrix
00139   TMatrixD *Vmatrix = new TMatrixD(totbins,totbins,Varray);
00140 
00141   if (printMatrix) Vmatrix->Print();
00142   //Vmatrix->Print();
00143 
00144   cout.precision(20);
00145   //Insert it (determ found for debugging purposes)!
00146   double determ;
00147   TMatrixD Vinvmatrix = Vmatrix->Invert(&determ);
00148   if (printMatrix) cout << "DETERM=" << determ << endl;
00149   cout.precision(5);
00150 
00151   if (printMatrix) Vinvmatrix.Print();
00152 
00153   //Extract the array of the inverted matrix:
00154   Double_t *Vinvarray = Vinvmatrix.GetMatrixArray();
00155 
00156   //Turn it into the inverted covariance matrix
00157   if (printMatrix) cout << "INVERTED" << endl;
00158   //make Vinv out of array
00159   for(i=0;i<totbins;i++)
00160   {
00161     for(j=0;j<totbins;j++)
00162     {
00163       InvErrorMatrix->SetBinContent(i+1,j+1,Vinvarray[i*totbins + j]*1.0e0);
00164     }
00165   }
00166 
00167   delete [] Varray;
00168   
00169   if (printMatrix) cout << "-------------" << endl;
00170   return;
00171 
00172 }

void NueFit2D_Joint::CombineFHCRHC (  ) 

Definition at line 196 of file NueFit2D_Joint.cxx.

References NueFit2D::Bkgd, ExtrapFHC, ExtrapRHC, potScaleFHC, potScaleRHC, and NueFit2D::Sig.

Referenced by DoDeltaFit(), DoNSIFitForever(), DoNSIFitQuick(), GetLikelihood(), GetMinLikelihood(), GetMinLikelihood_Delta(), GetNSIFFX2(), GetSensitivityAt(), RunDeltaChi2Contour(), and RunNSIDeltaChi2Contour().

00197 {
00198 
00199   Bkgd->Reset();
00200   Sig->Reset();
00201 
00202   int nFHC=0;
00203   int nRHC=0;
00204  
00205   if (ExtrapFHC.size()>0) nFHC = ExtrapFHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
00206   if (ExtrapRHC.size()>0) nRHC = ExtrapRHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
00207   unsigned int ie=0;
00208   TH1D* tmpsig = (TH1D*)Sig->Clone("tmpsig");
00209   tmpsig->Reset();
00210   TH1D* tmpbkg = (TH1D*)Bkgd->Clone("tmpbkg");
00211   tmpbkg->Reset();
00212 
00213   for(ie=0;ie<ExtrapFHC.size();ie++)
00214     {
00215       tmpbkg->Reset();
00216       tmpsig->Reset();
00217       for (int nb = 0; nb<nFHC; nb++){
00218         tmpbkg->SetBinContent(nb+1,potScaleFHC*(ExtrapFHC[ie]->Pred_TotalBkgd_VsBinNumber->GetBinContent(nb+1)));
00219         tmpsig->SetBinContent(nb+1,potScaleFHC*(ExtrapFHC[ie]->Pred_Signal_VsBinNumber->GetBinContent(nb+1)));
00220       }
00221       Bkgd->Add(tmpbkg);
00222       Sig->Add(tmpsig);
00223     }
00224 
00225   for(ie=0;ie<ExtrapRHC.size();ie++)
00226     {
00227       tmpbkg->Reset();
00228       tmpsig->Reset();
00229       for (int nb = 0; nb<nRHC; nb++){
00230         tmpbkg->SetBinContent(nb+1+nFHC,potScaleRHC*(ExtrapRHC[ie]->Pred_TotalBkgd_VsBinNumber->GetBinContent(nb+1)));
00231         tmpsig->SetBinContent(nb+1+nFHC,potScaleRHC*(ExtrapRHC[ie]->Pred_Signal_VsBinNumber->GetBinContent(nb+1)));
00232       }
00233       Bkgd->Add(tmpbkg);
00234       Sig->Add(tmpsig);
00235 
00236     }
00237 
00238 }

void NueFit2D_Joint::DefineBinDlnLMinuit (  )  [private, virtual]

Reimplemented from NueFit2D.

Definition at line 414 of file NueFit2D_Joint.cxx.

References ExtrapFHC, ExtrapRHC, NueFit2D::minuit, NueFit2D::nBins, and NueFit2D::NObs.

Referenced by DoDeltaFit(), DoNSIFitForever(), DoNSIFitQuick(), GetLikelihood(), GetMinLikelihood(), GetMinLikelihood_Delta(), GetNSIFFX2(), GetSensitivityAt(), RunDataGrid(), RunDeltaChi2Contour(), RunMultiBinFC_MHDeltaFit(), RunMultiBinPseudoExpts(), RunMultiBinPseudoExpts_MHDeltaFit(), and RunNSIDeltaChi2Contour().

00414                                         {
00415   //Function sets the maximum size of Minuit.
00416 
00417   //Clear things first:
00418 
00419   int npar = 0;
00420 
00421   if (nBins != 0){
00422     npar = nBins;
00423   } else if (NObs != 0){
00424       npar = NObs->GetNbinsX();
00425   } else if ((ExtrapFHC.size()+ExtrapRHC.size())!=0){
00426     npar = ExtrapFHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX()
00427      + ExtrapRHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
00428 
00429   } else {
00430     cout << "ERROR: Add extrapolation or NObs before initializing Minuit size"
00431          << endl;
00432     return;
00433   }
00434 
00435   //Make new minuit:
00436   minuit = new TMinuit(npar);
00437   minuit->SetPrintLevel(-1);
00438 
00439   double arglist[1];
00440   int ierflg=0;
00441 
00442   arglist[0] = 1.e-5;
00443   minuit->mnexcm("SET EPS",arglist,1,ierflg);
00444 
00445 
00446 }

void NueFit2D_Joint::DefineStdDlnLMinuit (  )  [private, virtual]

Reimplemented from NueFit2D.

Definition at line 448 of file NueFit2D_Joint.cxx.

References ExtrapFHC, ExtrapRHC, NueFit2D::FracErr_Bkgd_List, NueFit2D::minuit, NueFit2D::nBins, and NueFit2D::NObs.

Referenced by DoDeltaFit(), DoNSIFitForever(), DoNSIFitQuick(), GetLikelihood(), GetMinLikelihood(), GetMinLikelihood_Delta(), GetNSIFFX2(), GetSensitivityAt(), RunDataGrid(), RunDeltaChi2Contour(), RunMultiBinFC_MHDeltaFit(), RunMultiBinPseudoExpts(), RunMultiBinPseudoExpts_MHDeltaFit(), and RunNSIDeltaChi2Contour().

00448                                         {
00449   //Function sets the maximum size of Minuit.
00450 
00451   //Clear things first:
00452   //minuit->mncler();
00453 
00454   int npar = 0;
00455 
00456   //Number of systematics inputted:
00457   if (FracErr_Bkgd_List.size()!=0){
00458     npar = FracErr_Bkgd_List.size();
00459   } 
00460 
00461   int nb = 0;
00462   //Number of bins (for HOOHE):
00463   if (nBins != 0){
00464     nb = nBins;
00465   } else if (NObs != 0){
00466     nb = NObs->GetNbinsX();
00467   } else if ((ExtrapFHC.size()+ExtrapRHC.size())!=0){
00468     nb = ExtrapFHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX()
00469      + ExtrapRHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
00470   } else {
00471     cout << "ERROR: Add extrapolation or NObs before initializing Minuit size"
00472          << endl;
00473     return;
00474   }
00475 
00476   npar += nb;
00477 
00478   //make new minuit
00479   minuit = new TMinuit(npar+nb);
00480   minuit->SetPrintLevel(-1);
00481 
00482 }

void NueFit2D_Joint::DoDeltaFit (  ) 

Definition at line 3608 of file NueFit2D_Joint.cxx.

References MuELoss::a, NueFit2D::BinLikelihood(), NueFit2D::Bkgd, CombineFHCRHC(), DefineBinDlnLMinuit(), DefineStdDlnLMinuit(), NueFit2D::delta, NueFit2D::DeltaMSq21, NueFit2D::DeltaMSq32, NueFit2D::DeltaMSq32Unc, NueFit2D::ErrorMatrix, ExtrapFHC, ExtrapRHC, NueFit2D::FitMethod, NueFit2D::FracErr_Bkgd, NueFit2D::FracErr_Sig, Munits::g, NueFit2D::InvErrorMatrix, OscPar::kDelta, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kTh12, OscPar::kTh13, OscPar::kTh23, Munits::m, NueFit2D::nBins, NueFit2D::NExp, NueFit2D::NObs, NueFit2D::PoissonChi2(), NueFit2D::ScaledChi2(), NueFit2D::Sig, NueFit2D::StandardChi2(), NueFit2D::StandardLikelihood(), NueFit2D::Theta12, NueFit2D::Theta13, NueFit2D::Theta13Unc, NueFit2D::Theta23, and NueFit2D::Theta23Unc.

03609 {
03610   //probably don't want to set central values in the function call
03611   //make separate functions that also set the uncertainties, ie
03612   //NueFit2D::SetTheta13(theta13 central value, +1 sigma error, -1sigam error)
03613     
03614   if(NObs==0)
03615     {
03616       cout<<"NObs not set.  Quitting..."<<endl;
03617       return;
03618     }
03619   if(FitMethod==1 && (FracErr_Bkgd==0 || FracErr_Sig==0))
03620     {
03621       cout<<"FracErr_Bkgd and FracErr_Sig need to be set for ScaledChi2.  Quitting..."<<endl;
03622       return;
03623     }
03624   if((ExtrapFHC.size()+ExtrapRHC.size())==0)
03625     {
03626       cout<<"No Extrapolate2D input.  Quitting..."<<endl;
03627       return;
03628     }
03629   
03630   if(FitMethod==3) DefineStdDlnLMinuit();
03631   if(FitMethod==4) DefineBinDlnLMinuit();
03632   
03633   Bkgd = (TH1D*)NObs->Clone("Bkgd");
03634   Bkgd->Reset();
03635   Sig = (TH1D*)NObs->Clone("Sig");
03636   Sig->Reset();
03637   TH1D *NExp = (TH1D*)NObs->Clone("NExp");
03638   NExp->Reset();
03639   
03640   unsigned int ie;
03641   for(ie=0;ie<ExtrapFHC.size();ie++)
03642     {
03643       ExtrapFHC[ie]->GetPrediction();
03644       ExtrapFHC[ie]->SetOscPar(OscPar::kTh12,Theta12);
03645       ExtrapFHC[ie]->SetOscPar(OscPar::kTh23,Theta23);
03646       ExtrapFHC[ie]->SetOscPar(OscPar::kTh13,Theta13);
03647       ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM23,DeltaMSq32);
03648       ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM12,DeltaMSq21);
03649       ExtrapFHC[ie]->SetOscPar(OscPar::kDelta,delta);
03650       ExtrapFHC[ie]->OscillatePrediction();
03651     }
03652   for(ie=0;ie<ExtrapRHC.size();ie++)
03653     {
03654       ExtrapRHC[ie]->GetPrediction();
03655       ExtrapRHC[ie]->SetOscPar(OscPar::kTh12,Theta12);
03656       ExtrapRHC[ie]->SetOscPar(OscPar::kTh23,Theta23);
03657       ExtrapRHC[ie]->SetOscPar(OscPar::kTh13,Theta13);
03658       ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM23,DeltaMSq32);
03659       ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM12,DeltaMSq21);
03660       ExtrapRHC[ie]->SetOscPar(OscPar::kDelta,delta);
03661       ExtrapRHC[ie]->OscillatePrediction();
03662     }
03663   
03664   Int_t i;  
03665   if(ErrorMatrix==0) ErrorMatrix = new TH2D("ErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
03666   if(InvErrorMatrix==0) InvErrorMatrix = new TH2D("InvErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
03667   
03668   //Double_t ss2th13 = 0;
03669   Double_t delchi2 = 0;
03670   Double_t best = 0;
03671   Double_t ttbest = 0;
03672   Double_t otbest = 0;
03673   Double_t dmsbest = 0;
03674   //
03675   double chi2[3],val[3]; 
03676   TGraph *g3;
03677   TF1* fit;
03678   double minchi2;
03679   double global_minchi2;
03680   
03681   //find global minimum  
03682   //find minimum in delta
03683   best=-1.;
03684   minchi2=100;
03685   for(i=0;i<3;i++)
03686     {
03687       chi2[i]=-1.;
03688       val[i]=-1.;
03689     }
03690   
03691   double inc = 100.0;
03692   double Deltaincrement = (1/inc)*TMath::Pi(); 
03693   double db;
03694   for(i=0; i<2*inc+1 ; i++)//i scans over all delta values so we can find the best fit delta
03695     {
03696       chi2[0] = chi2[1];
03697       chi2[1] = chi2[2];
03698       val[0] = val[1];
03699       val[1] = val[2];
03700       
03701       db = i*Deltaincrement + 0.5;
03702       cout << "It: " << i << " with delta: " << db << endl;
03703       for(ie=0;ie<ExtrapFHC.size();ie++)
03704         {
03705           ExtrapFHC[ie]->SetOscPar(OscPar::kDelta,db);
03706           ExtrapFHC[ie]->OscillatePrediction();
03707         }
03708       
03709       for(ie=0;ie<ExtrapRHC.size();ie++)
03710         {
03711           ExtrapRHC[ie]->SetOscPar(OscPar::kDelta,db);
03712           ExtrapRHC[ie]->OscillatePrediction();
03713         }
03714       
03715       Bkgd->Reset();
03716       Sig->Reset();
03717       NExp->Reset();
03718       
03719       CombineFHCRHC();
03720       NExp->Add(Bkgd);
03721       NExp->Add(Sig);
03722       
03723       chi2[2] = 1e10;
03724       if(FitMethod==0)
03725         {
03726           chi2[2] = PoissonChi2(NExp);
03727         }
03728       else if(FitMethod==1)
03729         {
03730           chi2[2] = ScaledChi2(Bkgd,Sig);
03731         }
03732       else if(FitMethod==2)
03733         {
03734           chi2[2] = StandardChi2(NExp);
03735         }
03736       else if(FitMethod==3)
03737         {
03738           //Likelihood: "Standard" (N syst, N nuisance)
03739           //Calculate the likelihood (x2 for chi)
03740           chi2[2] = StandardLikelihood();
03741         }
03742       else if(FitMethod==4)
03743         {
03744           //Likelihood: Bin by Bin Calculation of Systematics
03745           //Calculate the likelihood (x2 for chi)
03746           chi2[2] = BinLikelihood();
03747           //cout << "ss2th13 = " << ss2th13 << " Chi2=" << chi2[2] << endl;
03748         }
03749       else
03750         {
03751           cout<<"Error in GetLikelihood(): Unknown 'FitMethod'."<<endl;
03752         }
03753       
03754       val[2] = db;
03755       
03756       if(i<2) continue;
03757       
03758       if(i<3 && chi2[2]>chi2[1] && chi2[1]>chi2[0])//first three points are increasing, first point is minimum.
03759         {
03760           best = val[0];
03761           minchi2 = chi2[0];
03762           cout<<"minimum at 1st point: "<<minchi2<<" at "<<best<<endl;
03763           break;
03764         }
03765       
03766       if(chi2[2]>chi2[1] && chi2[0]>chi2[1])//found minimum
03767         {
03768           g3 = new TGraph(3, val, chi2);
03769           fit = new TF1("pol2", "pol2");
03770           g3->Fit(fit, "Q");//fit to second order polynominal
03771           if(fit->GetParameter(2) > 0)//if the x^2 term is nonzero
03772             {
03773               best = -fit->GetParameter(1)/(2*fit->GetParameter(2));//the location of the minimum is -p1/(2*p2)
03774               minchi2 = fit->GetParameter(0) + fit->GetParameter(1)*best + fit->GetParameter(2)*best*best;
03775               cout<<"minimum with fit: "<<minchi2<<" at "<<best<<endl;
03776             }
03777           else//if the x^2 term is zero, then just use the minimum you got by scanning
03778             {
03779               best = val[1];
03780               minchi2 = chi2[1];
03781               cout<<"minimum with scan: "<<minchi2<<" at "<<best<<endl;
03782             }
03783           break;
03784         }
03785     }
03786   
03787   double f;
03788   minchi2=1e10;
03789   inc = 2.0;//20.0;
03790   int d,a,r,m;
03791   Deltaincrement = 0.01; //*TMath::Pi();
03792   double Aincrement =2*Theta23Unc/inc;
03793   double Rincrement =2*Theta13Unc/inc;
03794   double Mincrement = 0.01e-3;
03795   double das, aas, ras, mas;
03796   double dbest; 
03797   double penpar;
03798   
03799   double ass, rss;
03800   //NEW
03801   das = 0.0;
03802   ass = 0.3;
03803   rss = 0.0167;
03804   mas = -2.70e-3;
03805   
03806   //Delta Increment
03807   for(d=0;d<=10;d++){
03808     das = best - (5)*Deltaincrement + d*Deltaincrement;
03809     
03810     //Delta increment for Andy's combined surface
03811     //das = d*(2*TMath::Pi())/32.0;
03812     
03813     //Theta 23 increment
03814     for(a=0;a<=inc;a++){
03815       aas = Theta23 - (inc/2)*Aincrement + a*Aincrement;
03816       
03817       //Theta23 increment used for Andy's combined surface
03818       //for(a=0;a<8;a++){
03819       //ass = 0.3 + a*0.06;
03820       //aas = TMath::ASin(TMath::Sqrt(ass)); 
03821       
03822       //Theta13 increment
03823       for(r=0;r<=inc;r++){
03824         ras = Theta13 - (inc/2)*Rincrement+ (r)*Rincrement;
03825         
03826         //Theta13 increment used for Andy's combined surface  
03827         //for(r=0;r<7;r++){
03828         //rss = 0.0167 + r*(0.0317-0.0167)/6;
03829         // ras = TMath::ASin(TMath::Sqrt(rss));
03830         
03831         //Mass increment
03832         for(m=0;m<=10;m++){
03833           mas = DeltaMSq32 - (5)*Mincrement + (m)*Mincrement;
03834           
03835           //Mass increment for Andy's combined surface
03836           //for(m=0;m<7;m++){
03837           //mas = -2.70e-3 + m*0.1e-3;
03838           
03839           cout << "Finding minchi2 for: Delta = " << das << " Theta23= " << aas << " Theta13 = " << ras << " with dms= " << mas << endl;
03840           for(ie=0;ie<ExtrapFHC.size();ie++)
03841             {
03842               ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM23,mas);
03843               ExtrapFHC[ie]->SetOscPar(OscPar::kTh13,ras);
03844               ExtrapFHC[ie]->SetOscPar(OscPar::kTh23,aas); 
03845               ExtrapFHC[ie]->SetOscPar(OscPar::kDelta,das);
03846               //Th12 and DeltaM12 set for Andy's surface
03847               //ExtrapFHC[ie]->SetOscPar(OscPar::kTh12,0.587);
03848               //ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM12,7.54e-5);
03849               ExtrapFHC[ie]->OscillatePrediction();
03850             }
03851           
03852           for(ie=0;ie<ExtrapRHC.size();ie++)
03853             {
03854               ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM23,mas);
03855               ExtrapRHC[ie]->SetOscPar(OscPar::kTh13,ras);
03856               ExtrapRHC[ie]->SetOscPar(OscPar::kTh23,aas);
03857               ExtrapRHC[ie]->SetOscPar(OscPar::kDelta,das);
03858               //Th12 and DeltaM12 set for Andy's surface
03859               //ExtrapRHC[ie]->SetOscPar(OscPar::kTh12,0.587);
03860               //ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM12,7.54e-5);
03861               ExtrapRHC[ie]->OscillatePrediction();
03862             }
03863           
03864           Bkgd->Reset();
03865           Sig->Reset();
03866           NExp->Reset();
03867           
03868           CombineFHCRHC();
03869           NExp->Add(Bkgd);
03870           NExp->Add(Sig);
03871           
03872           penpar = ((TMath::Sin(2*aas))*(TMath::Sin(2*aas)) - ((TMath::Sin(2*Theta23))*(TMath::Sin(2*Theta23))))*((TMath::Sin(2*aas))*(TMath::Sin(2*aas)) - ((TMath::Sin(2*Theta23))*(TMath::Sin(2*Theta23))))/(0.036*0.036);
03873           
03874           f = BinLikelihood() + penpar + (ras - Theta13)*(ras-Theta13)/(Theta13Unc*Theta13Unc) + (mas - DeltaMSq32)*(mas-DeltaMSq32)/(DeltaMSq32Unc*DeltaMSq32Unc);
03875           
03876           //outfile for Andy's surface - to swap back, only use BinLikelihood above and remove penalty terms
03877           //outfile << mas << " " << ass << " " << rss << " " << das << " " << f << endl;    
03878           if(f<minchi2){
03879             minchi2=f;
03880             dbest = das;
03881             ttbest = aas;
03882             otbest = ras;
03883             dmsbest = mas;
03884           }
03885         }//m
03886       }//r
03887     }//a
03888   }//d
03889   
03890   // gm2 for -2LogL   
03891   global_minchi2 = 0; 
03892   
03893   cout << "Final Points - Theta13: " << otbest << " Theta23: " << ttbest << " with DMS: " << dmsbest << endl;  
03894 
03895   double plotinc = 20;
03896   int nn = int(2*plotinc+1); // number of values in vector deltaval. We'' use the same number to fill TH1* histogram.
03897   double deltaval[41];
03898   for(int i=0;i<nn;i++){
03899     deltaval[i]=i*(1/plotinc)*TMath::Pi();
03900   } 
03901   
03902   double deltatwologl[41];
03903   
03904   double dp;
03905   for(i=0; i<2*plotinc+1; i++)//looping over delta again, but now only the points we want to display on the delta plot
03906     {
03907       dp=deltaval[i];
03908       cout << "It: " << i << " with  delta: " << dp << endl;
03909       for(ie=0;ie<ExtrapFHC.size();ie++)
03910         {
03911           ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM23,dmsbest);
03912           ExtrapFHC[ie]->SetOscPar(OscPar::kTh13,otbest);
03913           ExtrapFHC[ie]->SetOscPar(OscPar::kTh23,ttbest);
03914           ExtrapFHC[ie]->SetOscPar(OscPar::kDelta,dp);
03915           ExtrapFHC[ie]->OscillatePrediction();
03916         }
03917       for(ie=0;ie<ExtrapRHC.size();ie++)
03918         {
03919           ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM23,dmsbest);
03920           ExtrapRHC[ie]->SetOscPar(OscPar::kTh13,otbest);
03921           ExtrapRHC[ie]->SetOscPar(OscPar::kTh23,ttbest);
03922           ExtrapRHC[ie]->SetOscPar(OscPar::kDelta,dp);
03923           ExtrapRHC[ie]->OscillatePrediction();
03924         }
03925       
03926       Bkgd->Reset();
03927       Sig->Reset();
03928       NExp->Reset();
03929       
03930       CombineFHCRHC();
03931       NExp->Add(Bkgd);
03932       NExp->Add(Sig);
03933       
03934       delchi2 = 1e10;
03935       if(FitMethod==0)
03936         {
03937           delchi2 = PoissonChi2(NExp) - global_minchi2;
03938         }
03939       else if(FitMethod==1)
03940         {
03941           delchi2 = ScaledChi2(Bkgd,Sig) - global_minchi2;
03942         }
03943       else if(FitMethod==2)
03944         {
03945           delchi2 = StandardChi2(NExp) - global_minchi2;
03946         }
03947       else if(FitMethod==3)
03948         {
03949           //Likelihood: "Standard" (N syst, N nuisance)
03950           //Calculate the likelihood (x2 for chi)
03951           delchi2 = StandardLikelihood() - global_minchi2;
03952         }
03953       else if(FitMethod==4)
03954         {
03955           //Likelihood: Bin by Bin Calculation of Systematics
03956           //Calculate the likelihood (x2 for chi)
03957           delchi2 = BinLikelihood() - global_minchi2;
03958         }
03959       else
03960         {
03961           cout<<"Error in GetLikelihood(): Unknown 'FitMethod'."<<endl;
03962         }
03963       deltatwologl[i] = delchi2;
03964     }//end loop over delta values
03965   
03966   //make a TGraph
03967   TFile *w = new TFile("DeltaOutput.root","RECREATE");
03968   TGraph *g = new TGraph(nn,deltaval,deltatwologl);
03969   //g->Draw("al");
03970   g->SetName("deltasweep");
03971   g->Write();
03972   w->Close();
03973   //outfile for surface 
03974   //outfile.close();
03975   return;
03976 }

void NueFit2D_Joint::DoIH ( int  doihinput = 0  )  [inline]

Definition at line 152 of file NueFit2D_Joint.h.

References NHIH.

00152                               {
00153       NHIH = doihinput;
00154     }

void NueFit2D_Joint::DoNSIFitForever (  ) 

Definition at line 4299 of file NueFit2D_Joint.cxx.

References NueFit2D::BinLikelihood(), NueFit2D::Bkgd, CombineFHCRHC(), DefineBinDlnLMinuit(), DefineStdDlnLMinuit(), NueFit2D::delta, NueFit2D::DeltaHigh, NueFit2D::DeltaLow, NueFit2D::epsetauHigh, NueFit2D::epsetauLow, NueFit2D::ErrorMatrix, ExtrapFHC, ExtrapRHC, NueFit2D::FitMethod, NueFit2D::FracErr_Bkgd, NueFit2D::FracErr_Sig, Munits::g, gSystem(), NueFit2D::InvErrorMatrix, NueFit2D::nBins, NueFit2D::nDeltaSteps, NueFit2D::nepsetauSteps, NueFit2D::NExp, NHIH, NueFit2D::NObs, NueFit2D::outFileName, NueFit2D::PoissonChi2(), NueFit2D::ScaledChi2(), NueFit2D::Sig, NueFit2D::StandardChi2(), and NueFit2D::StandardLikelihood().

04300 {
04301   //probably don't want to set central values in the function call
04302   //make separate functions that also set the uncertainties, ie
04303   //NueFit2D::SetTheta13(theta13 central value, +1 sigma error, -1sigam error)
04304     
04305   if(NObs==0)
04306     {
04307       cout<<"NObs not set.  Quitting..."<<endl;
04308       return;
04309     }
04310   if(FitMethod==1 && (FracErr_Bkgd==0 || FracErr_Sig==0))
04311     {
04312       cout<<"FracErr_Bkgd and FracErr_Sig need to be set for ScaledChi2.  Quitting..."<<endl;
04313       return;
04314     }
04315   if((ExtrapFHC.size()+ExtrapRHC.size())==0)
04316     {
04317       cout<<"No Extrapolate2D input.  Quitting..."<<endl;
04318       return;
04319     }
04320   
04321   if(FitMethod==3) DefineStdDlnLMinuit();
04322   if(FitMethod==4) DefineBinDlnLMinuit();
04323   
04324   Bkgd = (TH1D*)NObs->Clone("Bkgd");
04325   Bkgd->Reset();
04326   Sig = (TH1D*)NObs->Clone("Sig");
04327   Sig->Reset();
04328   TH1D *NExp = (TH1D*)NObs->Clone("NExp");
04329   NExp->Reset();
04330   
04331   unsigned int ie;
04332   for(ie=0;ie<ExtrapFHC.size();ie++)
04333     {
04334       ExtrapFHC[ie]->GetPrediction();
04335     }
04336   for(ie=0;ie<ExtrapRHC.size();ie++)
04337     {
04338       ExtrapRHC[ie]->GetPrediction();
04339     }
04340   
04341   
04342   if(ErrorMatrix==0) ErrorMatrix = new TH2D("ErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
04343   if(InvErrorMatrix==0) InvErrorMatrix = new TH2D("InvErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
04344   
04345   Int_t i,j,k,q,s;
04346   i = 1;
04347   Int_t idx = 0;
04348   Double_t epsetau;
04349   Double_t delta = 0;
04350   Double_t Deltaincrement = 0;
04351   if(nDeltaSteps>0) Deltaincrement = (DeltaHigh - DeltaLow)/(nDeltaSteps);
04352   Double_t epsetauincrement = 0;
04353   if(nepsetauSteps>0) epsetauincrement = (epsetauHigh - epsetauLow)/(nepsetauSteps);
04354   int arraysize = (nDeltaSteps+1)*(nepsetauSteps+1);  
04355   cout << arraysize << endl;
04356   double x[18045],y[18045],chi[18045]; 
04357   //double x[53489],y[53489],chi[53489];  
04358   TF1* fit;
04359   double minchi2;
04360   double mc2;
04361   
04362   if(NHIH==1){
04363   for(ie=0;ie<ExtrapFHC.size();ie++) ExtrapFHC[ie]->InvertMassHierarchy();
04364   for(ie=0;ie<ExtrapRHC.size();ie++) ExtrapRHC[ie]->InvertMassHierarchy();
04365   }//DoIH function sets inverted hierarchy
04366   
04367  for(j=0;j<nDeltaSteps+1;j++)
04368   {
04369     delta = j*Deltaincrement*TMath::Pi() + DeltaLow*TMath::Pi();
04370     for(ie=0;ie<ExtrapFHC.size();ie++){
04371       ExtrapFHC[ie]->SetDelta_etau(delta);
04372     }
04373     for(ie=0;ie<ExtrapRHC.size();ie++){
04374       ExtrapRHC[ie]->SetDelta_etau(delta);
04375     }
04376     
04377     minchi2=100;
04378     mc2=1000;
04379 
04380     for(q=0;q<=nepsetauSteps;q++){
04381       if(idx>arraysize){
04382         cout << "The code is not doing what you think it is doing, check arraysize and idx." << endl;
04383        }
04384       epsetau = epsetauLow+q*epsetauincrement;
04385       for(ie=0;ie<ExtrapFHC.size();ie++)
04386       {
04387         ExtrapFHC[ie]->SetEps_etau(epsetau);
04388         ExtrapFHC[ie]->OscillatePrediction();
04389       }
04390       for(ie=0;ie<ExtrapRHC.size();ie++)
04391       {
04392         ExtrapRHC[ie]->SetEps_etau(epsetau);
04393         ExtrapRHC[ie]->OscillatePrediction();
04394       }
04395 
04396       Bkgd->Reset();
04397       Sig->Reset();
04398       NExp->Reset();
04399 
04400       CombineFHCRHC();
04401       NExp->Add(Bkgd);
04402       NExp->Add(Sig);
04403       
04404       if(FitMethod==0)
04405       {
04406         x[idx] = epsetau;
04407         y[idx] = delta/TMath::Pi();
04408         chi[idx] = PoissonChi2(NExp);
04409       }
04410       else if(FitMethod==1)
04411       {
04412         x[idx] = epsetau;
04413         y[idx] = delta/TMath::Pi();
04414         chi[idx] = ScaledChi2(Bkgd,Sig);
04415       }
04416       else if(FitMethod==2)
04417       {
04418         x[idx] = epsetau;
04419         y[idx] = delta/TMath::Pi();
04420         chi[idx] = StandardChi2(NExp);
04421       }
04422       else if(FitMethod==3)
04423       {
04424         //Likelihood: "Standard" (N syst, N nuisance)
04425         //Calculate the likelihood (x2 for chi)
04426         x[idx] = epsetau;
04427         y[idx] = delta/TMath::Pi();
04428         chi[idx] = StandardLikelihood();
04429       }
04430       else if(FitMethod==4)
04431 
04432       {
04433         //Likelihood: Bin by Bin Calculation of Systematics
04434         //Calculate the likelihood (x2 for chi)
04435         //cout << idx << " " << epsetau << " " << delta << endl;
04436         x[idx] = epsetau;
04437         y[idx] = delta/TMath::Pi();
04438         chi[idx] = BinLikelihood();
04439       }
04440       else
04441       {
04442         cout<<"Error in RunNSIDeltaChi2Contour(): Unknown 'FitMethod'."<<endl;
04443       }
04444     
04445     idx++;    
04446   
04447     }//epsetauloop 
04448  
04449   }//deltaloop
04450  
04451  double dbest[1];
04452  double epsbest[1];
04453  for(k=0;k<arraysize;k++){
04454    if(chi[k]<mc2){
04455      mc2=chi[k];
04456      dbest[0]=y[k];
04457      epsbest[0]=x[k];
04458    }
04459  }
04460  cout << "Subtracting minchi value of: " << mc2 << " from coordinate ("<<epsbest[0]<<","<<dbest[0]<<"Pi)"<< endl;
04461 
04462  for(s=0;s<arraysize;s++){
04463    chi[s]=chi[s]-mc2;
04464    if(abs(chi[s])<0.00001){
04465      cout <<"I have a zero!"<<endl;
04466    }
04467    
04468  }
04469 
04470   //make a TGraph
04471   TFile *w = new TFile(gSystem->ExpandPathName(outFileName.c_str()),"RECREATE");
04472   TGraph *bp = new TGraph(i,epsbest,dbest);
04473   TGraph2D *g = new TGraph2D(arraysize,x,y,chi);
04474   g->GetYaxis()->SetTitle("#delta_{e#tau}/#pi");
04475   g->GetXaxis()->SetTitle("#epsilon_{e#tau}");
04476   g->GetZaxis()->SetTitle("-2#DeltalnL");
04477   bp->SetName("bp");
04478   g->SetName("ets");
04479   g->Write();
04480   bp->Write();
04481   w->Close();
04482   return;
04483 }

void NueFit2D_Joint::DoNSIFitQuick (  ) 

Definition at line 3980 of file NueFit2D_Joint.cxx.

References NueFit2D::BinLikelihood(), NueFit2D::Bkgd, CombineFHCRHC(), DefineBinDlnLMinuit(), DefineStdDlnLMinuit(), NueFit2D::ErrorMatrix, ExtrapFHC, ExtrapRHC, NueFit2D::FitMethod, NueFit2D::FracErr_Bkgd, NueFit2D::FracErr_Sig, Munits::g, NueFit2D::InvErrorMatrix, OscPar::kEps_etau, NueFit2D::nBins, NueFit2D::NExp, NueFit2D::NObs, NueFit2D::PoissonChi2(), NueFit2D::ScaledChi2(), NueFit2D::Sig, NueFit2D::StandardChi2(), and NueFit2D::StandardLikelihood().

03981 {
03982   //probably don't want to set central values in the function call
03983   //make separate functions that also set the uncertainties, ie
03984   //NueFit2D::SetTheta13(theta13 central value, +1 sigma error, -1sigam error)
03985   if(NObs==0)
03986     {
03987       cout<<"NObs not set.  Quitting..."<<endl;
03988       return;
03989     }
03990   if(FitMethod==1 && (FracErr_Bkgd==0 || FracErr_Sig==0))
03991     {
03992       cout<<"FracErr_Bkgd and FracErr_Sig need to be set for ScaledChi2.  Quitting..."<<endl;
03993       return;
03994     }
03995   if((ExtrapFHC.size()+ExtrapRHC.size())==0)
03996     {
03997       cout<<"No Extrapolate2D input.  Quitting..."<<endl;
03998       return;
03999     }
04000   
04001   if(FitMethod==3) DefineStdDlnLMinuit();
04002   if(FitMethod==4) DefineBinDlnLMinuit();
04003   
04004   Bkgd = (TH1D*)NObs->Clone("Bkgd");
04005   Bkgd->Reset();
04006   Sig = (TH1D*)NObs->Clone("Sig");
04007   Sig->Reset();
04008   TH1D *NExp = (TH1D*)NObs->Clone("NExp");
04009   NExp->Reset();
04010   
04011   unsigned int ie;
04012   for(ie=0;ie<ExtrapFHC.size();ie++)
04013     {
04014       ExtrapFHC[ie]->GetPrediction();
04015       /*
04016       ExtrapFHC[ie]->SetOscPar(OscPar::kTh12,Theta12);
04017       ExtrapFHC[ie]->SetOscPar(OscPar::kTh23,Theta23);
04018       ExtrapFHC[ie]->SetOscPar(OscPar::kTh13,Theta13);
04019       ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM23,DeltaMSq32);
04020       ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM12,DeltaMSq21);
04021       ExtrapFHC[ie]->SetOscPar(OscPar::kDelta,delta);
04022 
04023       // Add NSI parameters:
04024       ExtrapFHC[ie]->SetOscPar(OscPar::kEps_ee, Eps_ee);
04025       ExtrapFHC[ie]->SetOscPar(OscPar::kEps_emu, Eps_emu);
04026       ExtrapFHC[ie]->SetOscPar(OscPar::kEps_etau, Eps_etau);
04027       ExtrapFHC[ie]->SetOscPar(OscPar::kEps_mumu, Eps_mumu);
04028       ExtrapFHC[ie]->SetOscPar(OscPar::kEps_mutau, Eps_mutau); 
04029       ExtrapFHC[ie]->SetOscPar(OscPar::kEps_tautau, Eps_tautau); 
04030       ExtrapFHC[ie]->SetOscPar(OscPar::kDelta_emu, Delta_emu);
04031       ExtrapFHC[ie]->SetOscPar(OscPar::kDelta_etau, Delta_etau);
04032       ExtrapFHC[ie]->SetOscPar(OscPar::kDelta_mutau, Delta_mutau);
04033       //Oscillate prediction
04034       ExtrapFHC[ie]->OscillatePrediction();
04035       */
04036     }
04037   for(ie=0;ie<ExtrapRHC.size();ie++)
04038     {
04039       ExtrapRHC[ie]->GetPrediction();
04040       /*
04041       ExtrapRHC[ie]->SetOscPar(OscPar::kTh12,Theta12);
04042       ExtrapRHC[ie]->SetOscPar(OscPar::kTh23,Theta23);
04043       ExtrapRHC[ie]->SetOscPar(OscPar::kTh13,Theta13);
04044       ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM23,DeltaMSq32);
04045       ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM12,DeltaMSq21);
04046       ExtrapRHC[ie]->SetOscPar(OscPar::kDelta,delta);      
04047       // Add NSI parameters:
04048       ExtrapRHC[ie]->SetOscPar(OscPar::kEps_ee, Eps_ee);
04049       ExtrapRHC[ie]->SetOscPar(OscPar::kEps_emu, Eps_emu);
04050       ExtrapRHC[ie]->SetOscPar(OscPar::kEps_etau, Eps_etau);
04051       ExtrapRHC[ie]->SetOscPar(OscPar::kEps_mumu, Eps_mumu);
04052       ExtrapRHC[ie]->SetOscPar(OscPar::kEps_mutau, Eps_mutau); 
04053       ExtrapRHC[ie]->SetOscPar(OscPar::kEps_tautau, Eps_tautau); 
04054       ExtrapRHC[ie]->SetOscPar(OscPar::kDelta_emu, Delta_emu);
04055       ExtrapRHC[ie]->SetOscPar(OscPar::kDelta_etau, Delta_etau);
04056       ExtrapRHC[ie]->SetOscPar(OscPar::kDelta_mutau, Delta_mutau);
04057       // Oscillate Prediction
04058       ExtrapRHC[ie]->OscillatePrediction();
04059       */
04060     }
04061   
04062   Int_t i;  
04063   if(ErrorMatrix==0) ErrorMatrix = new TH2D("ErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
04064   if(InvErrorMatrix==0) InvErrorMatrix = new TH2D("InvErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
04065   
04066   //Double_t ss2th13 = 0;
04067   Double_t delchi2 = 0;
04068   Double_t best = 0;
04069   //
04070   double chi2[601],val[601]; 
04071   TGraph *g3;
04072   TF1* fit;
04073   double minchi2;
04074   double global_minchi2;
04075   
04076   //find global minimum  
04077   //find minimum in delta
04078   best=-1.;
04079   minchi2=1.0e6;
04080 
04081   double inc = 600.0;
04082   for(i=0;i<=inc;i++)
04083     {
04084       chi2[i]=-1.;
04085       val[i]=-1.;
04086     }
04087   
04088   double epsetauincrement = (12.0/inc);
04089   double nn = 601; 
04090   double epsetaulog[601];
04091   double epsetauval[601];
04092   for(int thing=0;thing<=inc;thing++){
04093     epsetauval[thing]=thing*epsetauincrement-6.0;
04094   }
04095 
04096   double db;
04097   for(i=0; i<inc+1 ; i++)//i scans over all delta values so we can find the best fit delta
04098     {
04099       /*
04100       chi2[0] = chi2[1];
04101       chi2[1] = chi2[2];
04102       val[0] = val[1];
04103       val[1] = val[2];
04104       */
04105       db = i*epsetauincrement - 6.0;
04106   
04107       cout << "It: " << i << " with eps_etau: " << db << endl;
04108       for(ie=0;ie<ExtrapFHC.size();ie++)
04109         {
04110           ExtrapFHC[ie]->SetOscPar(OscPar::kEps_etau,db);
04111           ExtrapFHC[ie]->OscillatePrediction();
04112         }
04113       
04114       for(ie=0;ie<ExtrapRHC.size();ie++)
04115         {
04116           ExtrapRHC[ie]->SetOscPar(OscPar::kEps_etau,db);
04117           ExtrapRHC[ie]->OscillatePrediction();
04118         }
04119      
04120       Bkgd->Reset();
04121       Sig->Reset();
04122       NExp->Reset();
04123       
04124       CombineFHCRHC();
04125       NExp->Add(Bkgd);
04126       NExp->Add(Sig);
04127       
04128 
04129       if(FitMethod==0)
04130         {
04131           chi2[i] = PoissonChi2(NExp);
04132         }
04133       else if(FitMethod==1)
04134         {
04135           chi2[i] = ScaledChi2(Bkgd,Sig);
04136         }
04137       else if(FitMethod==2)
04138         {
04139           chi2[i] = StandardChi2(NExp);
04140         }
04141       else if(FitMethod==3)
04142         {
04143           //Likelihood: "Standard" (N syst, N nuisance)
04144           //Calculate the likelihood (x2 for chi)
04145           chi2[i] = StandardLikelihood();
04146         }
04147       else if(FitMethod==4)
04148         {
04149           //Likelihood: Bin by Bin Calculation of Systematics
04150           //Calculate the likelihood (x2 for chi)
04151           chi2[i] = BinLikelihood();
04152           //cout << "ss2th13 = " << ss2th13 << " Chi2=" << chi2[2] << endl;
04153         }
04154       else
04155         {
04156           cout<<"Error in GetLikelihood(): Unknown 'FitMethod'."<<endl;
04157         }
04158       
04159       val[i] = db;
04160 
04161 
04162       /*      
04163       if(i<2) continue;
04164       
04165       if(i<3 && chi2[2]>chi2[1] && chi2[1]>chi2[0])//first three points are increasing, first point is minimum.
04166         {
04167           best = val[0];
04168           minchi2 = chi2[0];
04169           cout<<"minimum at 1st point: "<<minchi2<<" at "<<best<<endl;
04170           break;
04171         }
04172       
04173       if(chi2[2]>chi2[1] && chi2[0]>chi2[1])//found minimum
04174         {
04175           g3 = new TGraph(3, val, chi2);
04176           fit = new TF1("pol2", "pol2");
04177           g3->Fit(fit, "Q");//fit to second order polynominal
04178           if(fit->GetParameter(2) > 0)//if the x^2 term is nonzero
04179             {
04180               best = -fit->GetParameter(1)/(2*fit->GetParameter(2));//the location of the minimum is -p1/(2*p2)
04181               minchi2 = fit->GetParameter(0) + fit->GetParameter(1)*best + fit->GetParameter(2)*best*best;
04182               cout<<"minimum with fit: "<<minchi2<<" at "<<best<<endl;
04183             }
04184           else//if the x^2 term is zero, then just use the minimum you got by scanning
04185             {
04186               best = val[1];
04187               minchi2 = chi2[1];
04188               cout<<"minimum with scan: "<<minchi2<<" at "<<best<<endl;
04189             }
04190           break;
04191         }
04192       */
04193     }
04194   int idx;
04195  
04196   for(int bst=0; bst<=inc; bst++){
04197     if(chi2[bst]<minchi2){
04198       cout << "Chi2 of " <<chi2[bst]<< " at " <<val[bst] <<endl;
04199       minchi2 = chi2[bst];
04200       cout << "Minchi2 is " <<minchi2<<endl;
04201       best = val[bst];
04202       idx = bst;
04203     }
04204   }
04205 
04206   if(idx>1 && idx<inc-1){
04207       double va[3]; va[0]=val[idx-2]; va[1]=val[idx-1]; va[2]=val[idx];
04208       double ch[3]; ch[0]=chi2[idx-2]; ch[1]=chi2[idx-1]; ch[2]=chi2[idx];
04209       g3 = new TGraph(3, va, ch);
04210       fit = new TF1("pol2", "pol2");
04211       g3->Fit(fit, "Q");//fit to second order polynominal
04212           if(fit->GetParameter(2) > 0)//if the x^2 term is nonzero
04213             {
04214               best = -fit->GetParameter(1)/(2*fit->GetParameter(2));//the location of the minimum is -p1/(2*p2)
04215               minchi2 = fit->GetParameter(0) + fit->GetParameter(1)*best + fit->GetParameter(2)*best*best;
04216               cout<<"minimum with fit: "<<minchi2<<" at "<<best<<endl;
04217             }
04218           else//if the x^2 term is zero, then just use the minimum you got by scanning
04219             {
04220               cout<<"Keeping scan best of: "<<minchi2<<" at "<<best<<endl;
04221             }
04222          
04223     }
04224  
04225 
04226   cout << "Fit 1 Complete" << endl;
04227   minchi2=1e10;
04228  global_minchi2=0;
04229  
04230 
04231   double dp;
04232   for(i=0; i<inc+1; i++){
04233     //looping over eps_Etau again, but now only the points we want to display on the eps plot
04234     
04235       dp=epsetauval[i];
04236       cout << "It: " << i << " with  eps_etau: " << dp << endl;
04237       for(ie=0;ie<ExtrapFHC.size();ie++)
04238         {
04239           ExtrapFHC[ie]->SetOscPar(OscPar::kEps_etau,dp);
04240           ExtrapFHC[ie]->OscillatePrediction();
04241         }
04242       for(ie=0;ie<ExtrapRHC.size();ie++)
04243         {
04244           ExtrapRHC[ie]->SetOscPar(OscPar::kEps_etau,dp);
04245           ExtrapRHC[ie]->OscillatePrediction();
04246         }
04247       Bkgd->Reset();
04248       Sig->Reset();
04249       NExp->Reset();
04250       
04251       CombineFHCRHC();
04252       NExp->Add(Bkgd);
04253       NExp->Add(Sig);
04254       
04255       delchi2 = 1e10;
04256       if(FitMethod==0)
04257         {
04258           delchi2 = PoissonChi2(NExp) - global_minchi2;
04259         }
04260       else if(FitMethod==1)
04261         {
04262           delchi2 = ScaledChi2(Bkgd,Sig) - global_minchi2;
04263         }
04264       else if(FitMethod==2)
04265         {
04266           delchi2 = StandardChi2(NExp) - global_minchi2;
04267         }
04268       else if(FitMethod==3)
04269         {
04270           //Likelihood: "Standard" (N syst, N nuisance)
04271           //Calculate the likelihood (x2 for chi)
04272           delchi2 = StandardLikelihood() - global_minchi2;
04273         }
04274       else if(FitMethod==4)
04275         {
04276           //Likelihood: Bin by Bin Calculation of Systematics
04277           //Calculate the likelihood (x2 for chi)
04278           delchi2 = BinLikelihood() - global_minchi2;
04279         }
04280       else
04281         {
04282           cout<<"Error in GetLikelihood(): Unknown 'FitMethod'."<<endl;
04283         }
04284       epsetaulog[i] = delchi2;
04285     }//end loop over delta values
04286   cout <<"Fit 2 Complete" << endl;
04287   //make a TGraph
04288   TFile *w = new TFile("NSIQuick.root","RECREATE");
04289   cout <<"Outfile Opened" << endl;
04290   TGraph *g = new TGraph(nn,epsetauval,epsetaulog);
04291   g->SetName("etausweep");
04292   g->GetXaxis()->SetTitle("#epsilon_{e#tau}");
04293   g->GetYaxis()->SetTitle("-2lnL");
04294   g->Write();
04295   w->Close();
04296    return;
04297 }

double NueFit2D_Joint::GetLikelihood ( double  t12 = 0.6,
double  t23 = 0.785398,
double  t13 = 0,
double  dm2_32 = 2.32e-3,
double  dm2_21 = 7.59e-5,
double  delta = 0 
) [virtual]

Reimplemented from NueFit2D.

Definition at line 1981 of file NueFit2D_Joint.cxx.

References NueFit2D::BinLikelihood(), NueFit2D::Bkgd, CombineFHCRHC(), DefineBinDlnLMinuit(), DefineStdDlnLMinuit(), NueFit2D::ErrorMatrix, ExtrapFHC, ExtrapRHC, NueFit2D::FitMethod, NueFit2D::FracErr_Bkgd, NueFit2D::FracErr_Sig, NueFit2D::InvErrorMatrix, OscPar::kDelta, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kTh12, OscPar::kTh13, OscPar::kTh23, NueFit2D::nBins, NueFit2D::NExp, NueFit2D::NObs, NueFit2D::nSinSq2Th13Steps, NueFit2D::PoissonChi2(), NueFit2D::ScaledChi2(), NueFit2D::Sig, NueFit2D::SinSq2Th13High, NueFit2D::SinSq2Th13Low, NueFit2D::StandardChi2(), and NueFit2D::StandardLikelihood().

01982 {
01983   if(NObs==0)
01984   {
01985     cout<<"NObs not set.  Quitting..."<<endl;
01986     return 0;
01987   }
01988   if(FitMethod==1 && (FracErr_Bkgd==0 || FracErr_Sig==0))
01989   {
01990     cout<<"FracErr_Bkgd and FracErr_Sig need to be set for ScaledChi2.  Quitting..."<<endl;
01991     return 0;
01992   }
01993   if((ExtrapFHC.size()+ExtrapRHC.size())==0)
01994   {
01995     cout<<"No Extrapolate2D input.  Quitting..."<<endl;
01996     return 0;
01997   }
01998   
01999   if(FitMethod==3) DefineStdDlnLMinuit();
02000   if(FitMethod==4) DefineBinDlnLMinuit();
02001   
02002   Bkgd = (TH1D*)NObs->Clone("Bkgd");
02003   Bkgd->Reset();
02004   Sig = (TH1D*)NObs->Clone("Sig");
02005   Sig->Reset();
02006   TH1D *NExp = (TH1D*)NObs->Clone("NExp");
02007   NExp->Reset();
02008   
02009   unsigned int ie;
02010   for(ie=0;ie<ExtrapFHC.size();ie++)
02011   {
02012     ExtrapFHC[ie]->GetPrediction();
02013     ExtrapFHC[ie]->SetOscPar(OscPar::kTh12,t12);
02014     ExtrapFHC[ie]->SetOscPar(OscPar::kTh23,t23);
02015     ExtrapFHC[ie]->SetOscPar(OscPar::kTh13,t13);
02016     ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM23,dm2_32);
02017     ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM12,dm2_21);
02018     ExtrapFHC[ie]->SetOscPar(OscPar::kDelta,delta);
02019     ExtrapFHC[ie]->OscillatePrediction();
02020   }
02021   for(ie=0;ie<ExtrapRHC.size();ie++)
02022   {
02023     ExtrapRHC[ie]->GetPrediction();
02024     ExtrapRHC[ie]->SetOscPar(OscPar::kTh12,t12);
02025     ExtrapRHC[ie]->SetOscPar(OscPar::kTh23,t23);
02026     ExtrapRHC[ie]->SetOscPar(OscPar::kTh13,t13);
02027     ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM23,dm2_32);
02028     ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM12,dm2_21);
02029     ExtrapRHC[ie]->SetOscPar(OscPar::kDelta,delta);
02030     ExtrapRHC[ie]->OscillatePrediction();
02031   }
02032   
02033   Int_t i;
02034   
02035   if(ErrorMatrix==0) ErrorMatrix = new TH2D("ErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
02036   if(InvErrorMatrix==0) InvErrorMatrix = new TH2D("InvErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
02037   
02038   Double_t ss2th13 = 0;
02039   Double_t delchi2 = 0;
02040   Double_t best = 0;
02041   Double_t Th13increment = 0;
02042   if(nSinSq2Th13Steps>0) Th13increment = (SinSq2Th13High - SinSq2Th13Low)/(nSinSq2Th13Steps);
02043   double chi2[3],val[3];
02044   TGraph *g3;
02045   TF1* fit;
02046   double minchi2;
02047   
02048   best=-1.;
02049   minchi2=100;
02050   for(i=0;i<3;i++)
02051   {
02052     chi2[i]=-1.;
02053     val[i]=-1.;
02054   }
02055   for(i=0;i<nSinSq2Th13Steps+1;i++)
02056   {
02057     chi2[0] = chi2[1];
02058     chi2[1] = chi2[2];
02059     val[0] = val[1];
02060     val[1] = val[2];
02061     ss2th13 = i*Th13increment + SinSq2Th13Low;
02062 
02063     for(ie=0;ie<ExtrapFHC.size();ie++)
02064     {
02065       ExtrapFHC[ie]->SetSinSq2Th13(ss2th13);
02066       ExtrapFHC[ie]->OscillatePrediction();
02067     }
02068 
02069     for(ie=0;ie<ExtrapRHC.size();ie++)
02070     {
02071       ExtrapRHC[ie]->SetSinSq2Th13(ss2th13);
02072       ExtrapRHC[ie]->OscillatePrediction();
02073     }
02074 
02075     Bkgd->Reset();
02076     Sig->Reset();
02077     NExp->Reset();
02078 
02079     CombineFHCRHC();
02080     NExp->Add(Bkgd);
02081     NExp->Add(Sig);
02082     
02083     chi2[2] = 1e10;
02084     if(FitMethod==0)
02085     {
02086       chi2[2] = PoissonChi2(NExp);
02087     }
02088     else if(FitMethod==1)
02089     {
02090       chi2[2] = ScaledChi2(Bkgd,Sig);
02091     }
02092     else if(FitMethod==2)
02093     {
02094       chi2[2] = StandardChi2(NExp);
02095     }
02096     else if(FitMethod==3)
02097     {
02098       //Likelihood: "Standard" (N syst, N nuisance)
02099       //Calculate the likelihood (x2 for chi)
02100       chi2[2] = StandardLikelihood();
02101     }
02102     else if(FitMethod==4)
02103     {
02104       //Likelihood: Bin by Bin Calculation of Systematics
02105       //Calculate the likelihood (x2 for chi)
02106       chi2[2] = BinLikelihood();
02107       //cout << "ss2th13 = " << ss2th13 << " Chi2=" << chi2[2] << endl;
02108     }
02109     else
02110     {
02111       cout<<"Error in GetLikelihood(): Unknown 'FitMethod'."<<endl;
02112     }
02113     
02114     val[2] = ss2th13;
02115     
02116     if(i<2) continue;
02117     
02118     if(i<3 && chi2[2]>chi2[1] && chi2[1]>chi2[0])//first three points are increasing, first point is minimum.
02119     {
02120       best = val[0];
02121       minchi2 = chi2[0];
02122       cout<<"minimum at 1st point: "<<minchi2<<" at "<<best<<endl;
02123       break;
02124     }
02125     
02126     if(chi2[2]>chi2[1] && chi2[0]>chi2[1])//found minimum
02127     {
02128       g3 = new TGraph(3, val, chi2);
02129       fit = new TF1("pol2", "pol2");
02130       g3->Fit(fit, "Q");//fit to second order polynominal
02131       if(fit->GetParameter(2) > 0)//if the x^2 term is nonzero
02132       {
02133         best = -fit->GetParameter(1)/(2*fit->GetParameter(2));//the location of the minimum is -p1/(2*p2)
02134         minchi2 = fit->GetParameter(0) + fit->GetParameter(1)*best + fit->GetParameter(2)*best*best;
02135         cout<<"minimum with fit: "<<minchi2<<" at "<<best<<endl;
02136       }
02137       else//if the x^2 term is zero, then just use the minimum you got by scanning
02138       {
02139         best = val[1];
02140         minchi2 = chi2[1];
02141         cout<<"minimum with scan: "<<minchi2<<" at "<<best<<endl;
02142       }
02143       break;
02144     }
02145   }
02146   
02147   for(ie=0;ie<ExtrapFHC.size();ie++)
02148   {
02149     ExtrapFHC[ie]->SetOscPar(OscPar::kTh13,t13);
02150     ExtrapFHC[ie]->OscillatePrediction();
02151   }
02152   for(ie=0;ie<ExtrapRHC.size();ie++)
02153   {
02154     ExtrapRHC[ie]->SetOscPar(OscPar::kTh13,t13);
02155     ExtrapRHC[ie]->OscillatePrediction();
02156   }
02157 
02158   Bkgd->Reset();
02159   Sig->Reset();
02160   NExp->Reset();
02161 
02162   CombineFHCRHC();
02163   NExp->Add(Bkgd);
02164   NExp->Add(Sig);
02165   
02166   delchi2 = 1e10;
02167   if(FitMethod==0)
02168   {
02169     delchi2 = PoissonChi2(NExp) - minchi2;
02170   }
02171   else if(FitMethod==1)
02172   {
02173     delchi2 = ScaledChi2(Bkgd,Sig) - minchi2;
02174   }
02175   else if(FitMethod==2)
02176   {
02177     delchi2 = StandardChi2(NExp) - minchi2;
02178   }
02179   else if(FitMethod==3)
02180   {
02181     //Likelihood: "Standard" (N syst, N nuisance)
02182     //Calculate the likelihood (x2 for chi)
02183     delchi2 = StandardLikelihood() - minchi2;
02184   }
02185   else if(FitMethod==4)
02186   {
02187     //Likelihood: Bin by Bin Calculation of Systematics
02188     //Calculate the likelihood (x2 for chi)
02189     delchi2 = BinLikelihood() - minchi2;
02190   }
02191   else
02192   {
02193     cout<<"Error in GetLikelihood(): Unknown 'FitMethod'."<<endl;
02194   }
02195   
02196   cout<<"delchi2 = "<<delchi2<<endl;
02197   
02198   return delchi2;
02199 }

double NueFit2D_Joint::GetMinLikelihood ( double  delta = 0,
bool  normalhier = true 
) [private, virtual]

Reimplemented from NueFit2D.

Definition at line 2201 of file NueFit2D_Joint.cxx.

References NueFit2D::BinLikelihood(), NueFit2D::Bkgd, CombineFHCRHC(), DefineBinDlnLMinuit(), DefineStdDlnLMinuit(), NueFit2D::ErrCalc, NueFit2D::ErrorMatrix, ExtrapFHC, ExtrapRHC, NueFit2D::FitMethod, NueFit2D::FracErr_Bkgd, NueFit2D::FracErr_Sig, Munits::g, NueFit2D::grid_i_dm2_21, NueFit2D::grid_i_dm2_32, NueFit2D::grid_i_th12, NueFit2D::grid_i_th23, NueFit2D::grid_n_dm2_21, NueFit2D::grid_n_dm2_32, NueFit2D::grid_n_th12, NueFit2D::grid_n_th23, NueFit2D::InvErrorMatrix, OscPar::kDelta, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kTh12, OscPar::kTh23, max, NueFit2D::nBins, nBinsFHC, nBinsRHC, NueFit2D::NExp, NueFit2D::NObs, Munits::ns, NueFit2D::PoissonChi2(), NueFit2D::ScaledChi2(), ErrorCalc::SetUseGrid(), NueFit2D::Sig, NueFit2D::StandardChi2(), and NueFit2D::StandardLikelihood().

Referenced by RunDataGrid(), and RunMultiBinPseudoExpts().

02202 {
02203   if(NObs==0)
02204   {
02205     cout<<"NObs not set.  Quitting..."<<endl;
02206     return 0;
02207   }
02208   if(FitMethod==1 && (FracErr_Bkgd==0 || FracErr_Sig==0))
02209   {
02210     cout<<"FracErr_Bkgd and FracErr_Sig need to be set for ScaledChi2.  Quitting..."<<endl;
02211     return 0;
02212   }
02213   if((ExtrapFHC.size()+ExtrapRHC.size())==0)
02214   {
02215     cout<<"No Extrapolate2D input.  Quitting..."<<endl;
02216     return 0;
02217   }
02218  
02219   if(FitMethod==3) DefineStdDlnLMinuit();
02220   if(FitMethod==4) DefineBinDlnLMinuit();
02221   if(FitMethod==3 || FitMethod==4)
02222   {
02223     if(ErrCalc!=0) ErrCalc->SetUseGrid(false);
02224   }
02225   
02226   Bkgd = (TH1D*)NObs->Clone("Bkgd");
02227   Bkgd->Reset();
02228   Sig = (TH1D*)NObs->Clone("Sig");
02229   Sig->Reset();
02230   TH1D *NExp = (TH1D*)NObs->Clone("NExp");
02231   NExp->Reset();
02232   
02233   int nBinsFHC=0;
02234   int nBinsRHC=0;
02235   if (ExtrapFHC.size()>0) nBinsFHC = ExtrapFHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
02236   if (ExtrapRHC.size()>0) nBinsRHC = ExtrapRHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
02237   nBins = nBinsFHC + nBinsRHC;  
02238   
02239   unsigned int ie;
02240   for(ie=0;ie<ExtrapFHC.size();ie++)
02241   {
02242 
02243     ExtrapFHC[ie]->GetPrediction();
02244     if(normalhier)
02245     {
02246       ExtrapFHC[ie]->SetOscPar(OscPar::kTh12,grid_n_th12);
02247       ExtrapFHC[ie]->SetOscPar(OscPar::kTh23,grid_n_th23);
02248       ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM12,grid_n_dm2_21);
02249       ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM23,grid_n_dm2_32);
02250     }
02251     else
02252     {
02253       ExtrapFHC[ie]->SetOscPar(OscPar::kTh12,grid_i_th12);
02254       ExtrapFHC[ie]->SetOscPar(OscPar::kTh23,grid_i_th23);
02255       ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM12,grid_i_dm2_21);
02256       ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM23,grid_i_dm2_32);
02257     }
02258     ExtrapFHC[ie]->SetOscPar(OscPar::kDelta,delta);
02259     ExtrapFHC[ie]->OscillatePrediction();
02260   }
02261 
02262   for(ie=0;ie<ExtrapRHC.size();ie++)
02263   {
02264 
02265     ExtrapRHC[ie]->GetPrediction();
02266     if(normalhier)
02267     {
02268       ExtrapRHC[ie]->SetOscPar(OscPar::kTh12,grid_n_th12);
02269       ExtrapRHC[ie]->SetOscPar(OscPar::kTh23,grid_n_th23);
02270       ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM12,grid_n_dm2_21);
02271       ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM23,grid_n_dm2_32);
02272     }
02273     else
02274     {
02275       ExtrapRHC[ie]->SetOscPar(OscPar::kTh12,grid_i_th12);
02276       ExtrapRHC[ie]->SetOscPar(OscPar::kTh23,grid_i_th23);
02277       ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM12,grid_i_dm2_21);
02278       ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM23,grid_i_dm2_32);
02279     }
02280     ExtrapRHC[ie]->SetOscPar(OscPar::kDelta,delta);
02281     ExtrapRHC[ie]->OscillatePrediction();
02282   }
02283   Int_t i;
02284   
02285   if(ErrorMatrix==0) ErrorMatrix = new TH2D("ErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
02286   if(InvErrorMatrix==0) InvErrorMatrix = new TH2D("InvErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
02287   
02288   Double_t ss2th13 = 0;
02289   TF1* fit;
02290   double minchi2;
02291   
02292   const int ns=30;
02293   double max=0.6;
02294   double width = max/ns;
02295   double th13[ns],c2[ns];
02296   for(i=0;i<ns;i++)
02297   {
02298 
02299     ss2th13 = i*width;
02300     th13[i] = ss2th13;
02301 
02302     for(ie=0;ie<ExtrapFHC.size();ie++)
02303       {
02304         ExtrapFHC[ie]->SetSinSq2Th13(ss2th13);
02305         ExtrapFHC[ie]->OscillatePrediction();
02306       }
02307     for(ie=0;ie<ExtrapRHC.size();ie++)
02308       {
02309         ExtrapRHC[ie]->SetSinSq2Th13(ss2th13);
02310         ExtrapRHC[ie]->OscillatePrediction();
02311       }
02312     Bkgd->Reset();
02313     Sig->Reset();
02314     NExp->Reset();
02315  
02316   //Turn extrapolations into a prediction:
02317     CombineFHCRHC();
02318     NExp->Add(Bkgd);
02319     NExp->Add(Sig);
02320  
02321     c2[i] = 1e10;
02322     if(FitMethod==0)
02323     {
02324       c2[i] = PoissonChi2(NExp);
02325     }
02326     else if(FitMethod==1)
02327     {
02328       c2[i] = ScaledChi2(Bkgd,Sig);
02329     }
02330     else if(FitMethod==2)
02331     {
02332       c2[i] = StandardChi2(NExp);
02333     }
02334     else if(FitMethod==3)
02335     {
02336       //Likelihood: "Standard" (N syst, N nuisance)
02337       //Calculate the likelihood (x2 for chi)
02338       c2[i] = StandardLikelihood();
02339     }
02340     else if(FitMethod==4)
02341     {
02342       //Likelihood: Bin by Bin Calculation of Systematics
02343       //Calculate the likelihood (x2 for chi)
02344       c2[i] = BinLikelihood();
02345     }
02346     else
02347     {
02348       cout<<"Error in GetMinLikelihood(): Unknown 'FitMethod'."<<endl;
02349     }
02350   }
02351   
02352   TGraph *g = new TGraph(ns,th13,c2);
02353   fit = new TF1("pol6", "pol6");
02354   g->Fit(fit, "Q");//fit to second order polynominal
02355   minchi2 = fit->GetMinimum(0,max);
02356   
02357   return minchi2;
02358 }

double NueFit2D_Joint::GetMinLikelihood_Delta ( bool  normalhier = true  )  [private, virtual]

Reimplemented from NueFit2D.

Definition at line 2359 of file NueFit2D_Joint.cxx.

References NueFit2D::BinLikelihood(), NueFit2D::Bkgd, CombineFHCRHC(), DefineBinDlnLMinuit(), DefineStdDlnLMinuit(), NueFit2D::delta, NueFit2D::ErrCalc, NueFit2D::ErrorMatrix, ExtrapFHC, ExtrapRHC, NueFit2D::FitMethod, NueFit2D::FracErr_Bkgd, NueFit2D::FracErr_Sig, NueFit2D::grid_i_dm2_21, NueFit2D::grid_i_dm2_32, NueFit2D::grid_i_th12, NueFit2D::grid_i_th13, NueFit2D::grid_i_th23, NueFit2D::grid_n_dm2_21, NueFit2D::grid_n_dm2_32, NueFit2D::grid_n_th12, NueFit2D::grid_n_th13, NueFit2D::grid_n_th23, NueFit2D::InvErrorMatrix, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kTh12, OscPar::kTh13, OscPar::kTh23, NueFit2D::nBins, nBinsFHC, nBinsRHC, NueFit2D::nDeltaSteps, NueFit2D::NExp, NueFit2D::NObs, NueFit2D::PoissonChi2(), NueFit2D::ScaledChi2(), ErrorCalc::SetUseGrid(), NueFit2D::Sig, NueFit2D::StandardChi2(), and NueFit2D::StandardLikelihood().

Referenced by RunMultiBinFC_MHDeltaFit(), and RunMultiBinPseudoExpts_MHDeltaFit().

02360 {
02361   if(NObs==0)
02362   {
02363     cout<<"NObs not set.  Quitting..."<<endl;
02364     return 0;
02365   }
02366   if(FitMethod==1 && (FracErr_Bkgd==0 || FracErr_Sig==0))
02367   {
02368     cout<<"FracErr_Bkgd and FracErr_Sig need to be set for ScaledChi2.  Quitting..."<<endl;
02369     return 0;
02370   }
02371   if((ExtrapFHC.size()+ExtrapRHC.size())==0)
02372   {
02373     cout<<"No Extrapolate2D input.  Quitting..."<<endl;
02374     return 0;
02375   }
02376   
02377   if(FitMethod==3) DefineStdDlnLMinuit();
02378   if(FitMethod==4) DefineBinDlnLMinuit();
02379   if(FitMethod==3 || FitMethod==4)
02380   {
02381     if(ErrCalc!=0) ErrCalc->SetUseGrid(false);
02382   }
02383   
02384   Bkgd = (TH1D*)NObs->Clone("Bkgd");
02385   Bkgd->Reset();
02386   Sig = (TH1D*)NObs->Clone("Sig");
02387   Sig->Reset();
02388   TH1D *NExp = (TH1D*)NObs->Clone("NExp");
02389   NExp->Reset();
02390   
02391   int nBinsFHC=0;
02392   int nBinsRHC=0;
02393   if (ExtrapFHC.size()>0) nBinsFHC = ExtrapFHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
02394   if (ExtrapRHC.size()>0) nBinsRHC = ExtrapRHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
02395   nBins = nBinsFHC + nBinsRHC;  
02396   
02397   unsigned int ie;
02398   for(ie=0;ie<ExtrapFHC.size();ie++)
02399   {
02400     
02401     ExtrapFHC[ie]->GetPrediction();
02402     if(normalhier)
02403     {
02404       ExtrapFHC[ie]->SetOscPar(OscPar::kTh13,grid_n_th13);
02405       ExtrapFHC[ie]->SetOscPar(OscPar::kTh12,grid_n_th12);
02406       ExtrapFHC[ie]->SetOscPar(OscPar::kTh23,grid_n_th23);
02407       ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM12,grid_n_dm2_21);
02408       ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM23,grid_n_dm2_32);
02409     }
02410     else
02411     {
02412       ExtrapFHC[ie]->SetOscPar(OscPar::kTh13,grid_i_th13);
02413       ExtrapFHC[ie]->SetOscPar(OscPar::kTh12,grid_i_th12);
02414       ExtrapFHC[ie]->SetOscPar(OscPar::kTh23,grid_i_th23);
02415       ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM12,grid_i_dm2_21);
02416       ExtrapFHC[ie]->SetOscPar(OscPar::kDeltaM23,grid_i_dm2_32);
02417     }
02418     ExtrapFHC[ie]->OscillatePrediction();
02419   }
02420   
02421   for(ie=0;ie<ExtrapRHC.size();ie++)
02422   {
02423     
02424     ExtrapRHC[ie]->GetPrediction();
02425     if(normalhier)
02426     {
02427       ExtrapRHC[ie]->SetOscPar(OscPar::kTh13,grid_n_th13);
02428       ExtrapRHC[ie]->SetOscPar(OscPar::kTh12,grid_n_th12);
02429       ExtrapRHC[ie]->SetOscPar(OscPar::kTh23,grid_n_th23);
02430       ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM12,grid_n_dm2_21);
02431       ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM23,grid_n_dm2_32);
02432     }
02433     else
02434     {
02435       ExtrapRHC[ie]->SetOscPar(OscPar::kTh13,grid_i_th13);
02436       ExtrapRHC[ie]->SetOscPar(OscPar::kTh12,grid_i_th12);
02437       ExtrapRHC[ie]->SetOscPar(OscPar::kTh23,grid_i_th23);
02438       ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM12,grid_i_dm2_21);
02439       ExtrapRHC[ie]->SetOscPar(OscPar::kDeltaM23,grid_i_dm2_32);
02440     }
02441     ExtrapRHC[ie]->OscillatePrediction();
02442   }
02443   Int_t i;
02444   
02445   if(ErrorMatrix==0) ErrorMatrix = new TH2D("ErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
02446   if(InvErrorMatrix==0) InvErrorMatrix = new TH2D("InvErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
02447   
02448   Double_t delta = 0;
02449   //TF1* fit;
02450   double minchi2,mindelta;
02451   Double_t increment = 0;
02452   if(nDeltaSteps>0) increment = 2*TMath::Pi()/(nDeltaSteps);
02453   const int nd=nDeltaSteps+1;
02454   double *d = new double[nDeltaSteps+1];
02455   double *c2 = new double[nDeltaSteps+1];
02456   for(i=0;i<nd;i++)
02457   {
02458     delta = i*increment;
02459     d[i] = delta;
02460     
02461     for(ie=0;ie<ExtrapFHC.size();ie++)
02462     {
02463       ExtrapFHC[ie]->SetDeltaCP(delta);
02464       ExtrapFHC[ie]->OscillatePrediction();
02465     }
02466     for(ie=0;ie<ExtrapRHC.size();ie++)
02467     {
02468       ExtrapRHC[ie]->SetDeltaCP(delta);
02469       ExtrapRHC[ie]->OscillatePrediction();
02470     }
02471     Bkgd->Reset();
02472     Sig->Reset();
02473     NExp->Reset();
02474     
02475     //Turn extrapolations into a prediction:
02476     CombineFHCRHC();
02477     NExp->Add(Bkgd);
02478     NExp->Add(Sig);
02479     
02480     c2[i] = 1e10;
02481     if(FitMethod==0)
02482     {
02483       c2[i] = PoissonChi2(NExp);
02484     }
02485     else if(FitMethod==1)
02486     {
02487       c2[i] = ScaledChi2(Bkgd,Sig);
02488     }
02489     else if(FitMethod==2)
02490     {
02491       c2[i] = StandardChi2(NExp);
02492     }
02493     else if(FitMethod==3)
02494     {
02495       //Likelihood: "Standard" (N syst, N nuisance)
02496       //Calculate the likelihood (x2 for chi)
02497       c2[i] = StandardLikelihood();
02498     }
02499     else if(FitMethod==4)
02500     {
02501       //Likelihood: Bin by Bin Calculation of Systematics
02502       //Calculate the likelihood (x2 for chi)
02503       c2[i] = BinLikelihood();
02504     }
02505     else
02506     {
02507       cout<<"Error in GetMinLikelihood_Delta(): Unknown 'FitMethod'."<<endl;
02508     }
02509   }
02510   
02511   //can we actually fit this to a nice function?
02512   //   TGraph *g = new TGraph(ns,th13,c2);
02513   //   fit = new TF1("pol6", "pol6");
02514   //   g->Fit(fit, "Q");//fit to second order polynominal
02515   //   minchi2 = fit->GetMinimum(0,max);
02516   
02517   minchi2=1000;
02518   mindelta=0;
02519   for(i=0;i<nd;i++)
02520   {
02521     if(c2[i]<minchi2)
02522     {
02523       minchi2 = c2[i];
02524       mindelta = d[i];
02525     }
02526   }
02527   
02528   cout<<"min delta = "<<mindelta<<" , min chi2 = "<<minchi2<<endl;
02529   
02530   delete [] d;
02531   delete [] c2;
02532   
02533   return minchi2;
02534 }

void NueFit2D_Joint::GetNSIFFX2 (  ) 

Definition at line 4485 of file NueFit2D_Joint.cxx.

References NueFit2D::BinLikelihood(), NueFit2D::Bkgd, CombineFHCRHC(), DefineBinDlnLMinuit(), DefineStdDlnLMinuit(), NueFit2D::delta, NueFit2D::DeltaHigh, NueFit2D::DeltaLow, NueFit2D::epsetauHigh, NueFit2D::epsetauLow, NueFit2D::ErrorMatrix, ExtrapFHC, ExtrapRHC, NueFit2D::FitMethod, NueFit2D::FracErr_Bkgd, NueFit2D::FracErr_Sig, NueFit2D::InvErrorMatrix, NueFit2D::nBins, NueFit2D::nDeltaSteps, NueFit2D::nepsetauSteps, NueFit2D::NExp, NHIH, NueFit2D::NObs, NueFit2D::PoissonChi2(), NueFit2D::ScaledChi2(), NueFit2D::Sig, NueFit2D::StandardChi2(), and NueFit2D::StandardLikelihood().

04486 {
04487   //probably don't want to set central values in the function call
04488   //make separate functions that also set the uncertainties, ie
04489   //NueFit2D::SetTheta13(theta13 central value, +1 sigma error, -1sigam error)
04490     
04491   if(NObs==0)
04492     {
04493       cout<<"NObs not set.  Quitting..."<<endl;
04494       return;
04495     }
04496   if(FitMethod==1 && (FracErr_Bkgd==0 || FracErr_Sig==0))
04497     {
04498       cout<<"FracErr_Bkgd and FracErr_Sig need to be set for ScaledChi2.  Quitting..."<<endl;
04499       return;
04500     }
04501   if((ExtrapFHC.size()+ExtrapRHC.size())==0)
04502     {
04503       cout<<"No Extrapolate2D input.  Quitting..."<<endl;
04504       return;
04505     }
04506   
04507   if(FitMethod==3) DefineStdDlnLMinuit();
04508   if(FitMethod==4) DefineBinDlnLMinuit();
04509   
04510   Bkgd = (TH1D*)NObs->Clone("Bkgd");
04511   Bkgd->Reset();
04512   Sig = (TH1D*)NObs->Clone("Sig");
04513   Sig->Reset();
04514   TH1D *NExp = (TH1D*)NObs->Clone("NExp");
04515   NExp->Reset();
04516   
04517   unsigned int ie;
04518   for(ie=0;ie<ExtrapFHC.size();ie++)
04519     {
04520       ExtrapFHC[ie]->GetPrediction();
04521     }
04522   for(ie=0;ie<ExtrapRHC.size();ie++)
04523     {
04524       ExtrapRHC[ie]->GetPrediction();
04525     }
04526   
04527   
04528   if(ErrorMatrix==0) ErrorMatrix = new TH2D("ErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
04529   if(InvErrorMatrix==0) InvErrorMatrix = new TH2D("InvErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
04530   
04531   Int_t i,j,k,q,s;
04532   i = 1;
04533   Int_t idx = 0;
04534   Double_t epsetau;
04535   Double_t delta = 0;
04536   Double_t Deltaincrement = 0;
04537   if(nDeltaSteps>0) Deltaincrement = (DeltaHigh - DeltaLow)/(nDeltaSteps);
04538   Double_t epsetauincrement = 0;
04539   if(nepsetauSteps>0) epsetauincrement = (epsetauHigh - epsetauLow)/(nepsetauSteps);
04540   int arraysize = 1;  
04541   cout << arraysize << endl;
04542   double x[1],y[1],chi[1]; 
04543    TF1* fit;
04544   double minchi2;
04545   double mc2;
04546   
04547   if(NHIH==1){
04548   for(ie=0;ie<ExtrapFHC.size();ie++) ExtrapFHC[ie]->InvertMassHierarchy();
04549   for(ie=0;ie<ExtrapRHC.size();ie++) ExtrapRHC[ie]->InvertMassHierarchy();
04550   }//DoIH function sets inverted hierarchy
04551   
04552  for(j=0;j<nDeltaSteps+1;j++)
04553   {
04554     delta = j*Deltaincrement*TMath::Pi() + DeltaLow*TMath::Pi();
04555     for(ie=0;ie<ExtrapFHC.size();ie++){
04556       ExtrapFHC[ie]->SetDelta_etau(delta);
04557     }
04558     for(ie=0;ie<ExtrapRHC.size();ie++){
04559       ExtrapRHC[ie]->SetDelta_etau(delta);
04560     }
04561     
04562     minchi2=100;
04563     mc2=1000;
04564 
04565     for(q=0;q<=nepsetauSteps;q++){
04566       if(idx>arraysize){
04567         cout << "The code is not doing what you think it is doing, check arraysize and idx." << endl;
04568        }
04569       epsetau = epsetauLow+q*epsetauincrement;
04570       for(ie=0;ie<ExtrapFHC.size();ie++)
04571       {
04572         ExtrapFHC[ie]->SetEps_etau(epsetau);
04573         ExtrapFHC[ie]->OscillatePrediction();
04574       }
04575       for(ie=0;ie<ExtrapRHC.size();ie++)
04576       {
04577         ExtrapRHC[ie]->SetEps_etau(epsetau);
04578         ExtrapRHC[ie]->OscillatePrediction();
04579       }
04580 
04581       Bkgd->Reset();
04582       Sig->Reset();
04583       NExp->Reset();
04584 
04585       CombineFHCRHC();
04586       NExp->Add(Bkgd);
04587       NExp->Add(Sig);
04588       
04589       if(FitMethod==0)
04590       {
04591         x[idx] = epsetau;
04592         y[idx] = delta/TMath::Pi();
04593         chi[idx] = PoissonChi2(NExp);
04594       }
04595       else if(FitMethod==1)
04596       {
04597         x[idx] = epsetau;
04598         y[idx] = delta/TMath::Pi();
04599         chi[idx] = ScaledChi2(Bkgd,Sig);
04600       }
04601       else if(FitMethod==2)
04602       {
04603         x[idx] = epsetau;
04604         y[idx] = delta/TMath::Pi();
04605         chi[idx] = StandardChi2(NExp);
04606       }
04607       else if(FitMethod==3)
04608       {
04609         //Likelihood: "Standard" (N syst, N nuisance)
04610         //Calculate the likelihood (x2 for chi)
04611         x[idx] = epsetau;
04612         y[idx] = delta/TMath::Pi();
04613         chi[idx] = StandardLikelihood();
04614       }
04615       else if(FitMethod==4)
04616 
04617       {
04618         //Likelihood: Bin by Bin Calculation of Systematics
04619         //Calculate the likelihood (x2 for chi)
04620         //cout << idx << " " << epsetau << " " << delta << endl;
04621         x[idx] = epsetau;
04622         y[idx] = delta/TMath::Pi();
04623         chi[idx] = BinLikelihood();
04624       }
04625       else
04626       {
04627         cout<<"Error in RunNSIDeltaChi2Contour(): Unknown 'FitMethod'."<<endl;
04628       }
04629     
04630     idx++;    
04631   
04632     }//epsetauloop 
04633  
04634   }//deltaloop
04635  
04636  double dbest[1];
04637  double epsbest[1];
04638  for(k=0;k<arraysize;k++){
04639    if(chi[k]<mc2){
04640      mc2=chi[k];
04641      dbest[0]=y[k];
04642      epsbest[0]=x[k];
04643    }
04644  }
04645  cout << "Subtracting minchi value of: " << mc2 << " from coordinate ("<<epsbest[0]<<","<<dbest[0]<<"Pi)"<< endl;
04646 
04647  for(s=0;s<arraysize;s++){
04648    chi[s]=chi[s]-mc2;
04649    if(abs(chi[s])<0.00001){
04650      cout <<"I have a zero!"<<endl;
04651    }
04652    
04653  }
04654 
04655 }

double NueFit2D_Joint::GetSensitivityAt ( double  delta = 0,
bool  normalhier = true 
) [virtual]

Reimplemented from NueFit2D.

Definition at line 1215 of file NueFit2D_Joint.cxx.

References NueFit2D::BinLikelihood(), NueFit2D::Bkgd, CombineFHCRHC(), DefineBinDlnLMinuit(), DefineStdDlnLMinuit(), NueFit2D::ErrorMatrix, ExtrapFHC, ExtrapRHC, NueFit2D::FitMethod, NueFit2D::FracErr_Bkgd, NueFit2D::FracErr_Sig, NueFit2D::InvErrorMatrix, NueFit2D::nBins, NueFit2D::NExp, NueFit2D::NObs, NueFit2D::nSinSq2Th13Steps, NueFit2D::PoissonChi2(), NueFit2D::ScaledChi2(), NueFit2D::Sig, NueFit2D::SinSq2Th13High, NueFit2D::SinSq2Th13Low, NueFit2D::StandardChi2(), and NueFit2D::StandardLikelihood().

01216 {
01217 
01218   if(NObs==0)
01219   {
01220     cout<<"NObs not set.  Quitting..."<<endl;
01221     return 0;
01222   }
01223   if(FitMethod==1 && (FracErr_Bkgd==0 || FracErr_Sig==0))
01224   {
01225     cout<<"FracErr_Bkgd and FracErr_Sig need to be set for ScaledChi2.  Quitting..."<<endl;
01226     return 0;
01227   }
01228   if((ExtrapFHC.size()+ExtrapRHC.size())==0)
01229   {
01230     cout<<"No Extrapolate2D input.  Quitting..."<<endl;
01231     return 0;
01232   }
01233   
01234   if(FitMethod==3) DefineStdDlnLMinuit();
01235   if(FitMethod==4) DefineBinDlnLMinuit();
01236   
01237   Bkgd = (TH1D*)NObs->Clone("Bkgd");
01238   Bkgd->Reset();
01239   Sig = (TH1D*)NObs->Clone("Sig");
01240   Sig->Reset();
01241   TH1D *NExp = (TH1D*)NObs->Clone("NExp");
01242   NExp->Reset();
01243   
01244   unsigned int ie;
01245   for(ie=0;ie<ExtrapFHC.size();ie++) ExtrapFHC[ie]->GetPrediction();
01246   for(ie=0;ie<ExtrapRHC.size();ie++) ExtrapRHC[ie]->GetPrediction();
01247 
01248   Int_t i;
01249   
01250   if(ErrorMatrix==0) ErrorMatrix = new TH2D("ErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
01251   if(InvErrorMatrix==0) InvErrorMatrix = new TH2D("InvErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
01252   
01253   Double_t ss2th13 = 0;
01254   Double_t delchi2 = 0,delchi2prev = 0;
01255   Double_t contourlvl = 2.71;
01256   Double_t prev = 0;
01257   Double_t best = 0;
01258   Double_t Th13increment = 0;
01259   if(nSinSq2Th13Steps>0) Th13increment = (SinSq2Th13High - SinSq2Th13Low)/(nSinSq2Th13Steps);
01260   double sens=-1;
01261   double chi2[3],val[3];
01262   TGraph *g3;
01263   TF1* fit;
01264   double minchi2;
01265   double limit = 0;
01266   
01267   for(ie=0;ie<ExtrapRHC.size();ie++)
01268   {
01269     ExtrapRHC[ie]->SetDeltaCP(delta);
01270     if(!normalhier) ExtrapRHC[ie]->InvertMassHierarchy();
01271   }
01272 
01273   for(ie=0;ie<ExtrapFHC.size();ie++)
01274   {
01275     ExtrapFHC[ie]->SetDeltaCP(delta);
01276     if(!normalhier) ExtrapFHC[ie]->InvertMassHierarchy();
01277   }
01278   
01279   contourlvl = 2.71;
01280   
01281   best=-1.;
01282   minchi2=100;
01283   for(i=0;i<3;i++)
01284   {
01285     chi2[i]=-1.;
01286     val[i]=-1.;
01287   }
01288   for(i=0;i<nSinSq2Th13Steps+1;i++)
01289   {
01290     chi2[0] = chi2[1];
01291     chi2[1] = chi2[2];
01292     val[0] = val[1];
01293     val[1] = val[2];
01294     ss2th13 = i*Th13increment + SinSq2Th13Low;
01295 
01296     for(ie=0;ie<ExtrapFHC.size();ie++)
01297     {
01298       ExtrapFHC[ie]->SetSinSq2Th13(ss2th13);
01299       ExtrapFHC[ie]->OscillatePrediction();
01300     }
01301 
01302     for(ie=0;ie<ExtrapRHC.size();ie++)
01303     {
01304       ExtrapRHC[ie]->SetSinSq2Th13(ss2th13);
01305       ExtrapRHC[ie]->OscillatePrediction();
01306     }
01307 
01308     Bkgd->Reset();
01309     Sig->Reset();
01310     NExp->Reset();
01311 
01312     //Turn extrapolations into a prediction:
01313     CombineFHCRHC();
01314     NExp->Add(Bkgd);
01315     NExp->Add(Sig);
01316     
01317     chi2[2] = 1e10;
01318     if(FitMethod==0)
01319     {
01320       chi2[2] = PoissonChi2(NExp);
01321     }
01322     else if(FitMethod==1)
01323     {
01324       chi2[2] = ScaledChi2(Bkgd,Sig);
01325     }
01326     else if(FitMethod==2)
01327     {
01328       chi2[2] = StandardChi2(NExp);
01329     }
01330     else if(FitMethod==3)
01331     {
01332       //Likelihood: "Standard" (N syst, N nuisance)
01333       //Calculate the likelihood (x2 for chi)
01334       chi2[2] = StandardLikelihood();
01335     }
01336     else if(FitMethod==4)
01337     {
01338       //Likelihood: Bin by Bin Calculation of Systematics
01339       //Calculate the likelihood (x2 for chi)
01340       chi2[2] = BinLikelihood();
01341       //cout << ss2th13 << " -> " << chi2[2] << endl;
01342     }
01343     else
01344     {
01345       cout<<"Error in GetSensitivityAt(): Unknown 'FitMethod'."<<endl;
01346     }
01347     
01348     val[2] = ss2th13;
01349     
01350     if(i<2) continue;
01351     
01352     if(i<3 && chi2[2]>chi2[1] && chi2[1]>chi2[0])//first three points are increasing, first point is minimum.
01353     {
01354       best = val[0];
01355       minchi2 = chi2[0];
01356       cout<<"minimum at 1st point: "<<minchi2<<" at "<<best<<endl;
01357       break; //CHANGE BACK
01358     }
01359     
01360     if(chi2[2]>chi2[1] && chi2[0]>chi2[1])//found minimum
01361     {
01362       g3 = new TGraph(3, val, chi2);
01363       fit = new TF1("pol2", "pol2");
01364       g3->Fit(fit, "Q");//fit to second order polynominal
01365       if(fit->GetParameter(2) > 0)//if the x^2 term is nonzero
01366       {
01367         best = -fit->GetParameter(1)/(2*fit->GetParameter(2));//the location of the minimum is -p1/(2*p2)
01368         minchi2 = fit->GetParameter(0) + fit->GetParameter(1)*best + fit->GetParameter(2)*best*best;
01369         cout<<"minimum with fit: "<<minchi2<<" at "<<best<<endl;
01370       }
01371       else//if the x^2 term is zero, then just use the minimum you got by scanning
01372       {
01373         best = val[1];
01374         minchi2 = chi2[1];
01375         cout<<"minimum with scan: "<<minchi2<<" at "<<best<<endl;
01376       }
01377       break; //CHANGE BACK
01378     }
01379 
01380     //if (ss2th13>0.11) break; //DELETE ME
01381   }
01382   
01383   limit = -1;
01384   delchi2prev = 1000;
01385   prev = 0;
01386   for(i=0;i<nSinSq2Th13Steps+1;i++)
01387   {
01388     ss2th13 = i*Th13increment + SinSq2Th13Low;
01389     for(ie=0;ie<ExtrapFHC.size();ie++)
01390     {
01391       ExtrapFHC[ie]->SetSinSq2Th13(ss2th13);
01392       ExtrapFHC[ie]->OscillatePrediction();
01393     }
01394     for(ie=0;ie<ExtrapRHC.size();ie++)
01395     {
01396       ExtrapRHC[ie]->SetSinSq2Th13(ss2th13);
01397       ExtrapRHC[ie]->OscillatePrediction();
01398     }
01399     
01400     Bkgd->Reset();
01401     Sig->Reset();
01402     NExp->Reset();
01403     
01404     CombineFHCRHC();
01405     NExp->Add(Bkgd);
01406     NExp->Add(Sig);
01407     
01408     delchi2 = 1e10;
01409     if(FitMethod==0)
01410     {
01411       delchi2 = PoissonChi2(NExp) - minchi2;
01412     }
01413     else if(FitMethod==1)
01414     {
01415       delchi2 = ScaledChi2(Bkgd,Sig) - minchi2;
01416     }
01417     else if(FitMethod==2)
01418     {
01419       delchi2 = StandardChi2(NExp) - minchi2;
01420     }
01421     else if(FitMethod==3)
01422     {
01423       //Likelihood: "Standard" (N syst, N nuisance)
01424       //Calculate the likelihood (x2 for chi)
01425       delchi2 = StandardLikelihood() - minchi2;
01426     }
01427     else if(FitMethod==4)
01428     {
01429       //Likelihood: Bin by Bin Calculation of Systematics
01430       //Calculate the likelihood (x2 for chi)
01431       delchi2 = BinLikelihood() - minchi2;
01432       //cout << "ss2th13 = " << ss2th13 << " Chi2=" << delchi2 << endl;
01433     }
01434     else
01435     {
01436       cout<<"Error in GetSensitivityAt(): Unknown 'FitMethod'."<<endl;
01437     }
01438 
01439     //RBT: How to run the separate options
01440 
01441 //     if (opt==0){
01442 //       //Chi2 standard
01443 //       chi2 = StandardChi2(NExp);
01444 //     } else if (opt==1){
01445 //       //Likelihood: Bin by Bin Calculation of Systematics
01446 //       //Calculate the likelihood (x2 for chi)
01447 //       chi2 = BinLikelihood();
01448 //     } else if (opt==2){
01449 //       //Likelihood: "Standard" (N syst, N nuisance)
01450 //       //Calculate the likelihood (x2 for chi)
01451 //       chi2 = StandardLikelihood();
01452 //     }
01453 
01454     if(delchi2>contourlvl && delchi2prev<contourlvl && TMath::Abs(delchi2prev-delchi2)<0.1)
01455     {
01456       limit = ss2th13 + ((ss2th13-prev)/(delchi2 - delchi2prev))*(contourlvl - delchi2);
01457       sens = limit;
01458       break;
01459     }
01460     delchi2prev = delchi2;
01461     prev = ss2th13;
01462   }
01463   
01464   return sens;
01465 }

void NueFit2D_Joint::ReadGridFiles (  )  [virtual]

Reimplemented from NueFit2D.

Definition at line 240 of file NueFit2D_Joint.cxx.

References ExtrapFHC, ExtrapRHC, Form(), NueFit2D::grid_background, NueFit2D::grid_bin_oscparerr, NueFit2D::grid_bnuecc, NueFit2D::grid_delta, NueFit2D::grid_i_dm2_21, NueFit2D::grid_i_dm2_32, NueFit2D::grid_i_th12, NueFit2D::grid_i_th13, NueFit2D::grid_i_th23, NueFit2D::grid_n_dm2_21, NueFit2D::grid_n_dm2_32, NueFit2D::grid_n_th12, NueFit2D::grid_n_th13, NueFit2D::grid_n_th23, NueFit2D::grid_nc, NueFit2D::grid_nue, NueFit2D::grid_numucc, NueFit2D::grid_nutaucc, NueFit2D::grid_signal, NueFit2D::grid_sinsq2th13, NueFit2D::GridFileName_Inverted, NueFit2D::GridFileName_Normal, NueFit2D::GridNorm, GridScale_Inverted_FHC, GridScale_Inverted_RHC, GridScale_Normal_FHC, GridScale_Normal_RHC, NueFit2D::GridTree_2_Inverted, NueFit2D::GridTree_2_Normal, NueFit2D::GridTree_Inverted, NueFit2D::GridTree_Normal, gSystem(), NueFit2D::IncludeOscParErrs, NueFit2D::nBins, nBinsFHC, NueFit2D::nPts_Inverted, NueFit2D::nPts_Normal, NueFit2D::paramtree_Inverted, and NueFit2D::paramtree_Normal.

Referenced by RunDataGrid(), RunMultiBinFC_MHDeltaFit(), RunMultiBinPseudoExpts(), and RunMultiBinPseudoExpts_MHDeltaFit().

00241 {
00242   double fpfhc;
00243   double fprhc;
00244   unsigned int i;
00245   TTree *temp=0;
00246   unsigned int j;
00247   
00248   for(i=0;i<nBins;i++)
00249   {
00250     grid_bin_oscparerr.push_back(0);
00251   }
00252  
00253 
00254   TFile *fnorm;
00255   GridTree_Normal.clear();
00256   GridTree_2_Normal.clear();
00257   nPts_Normal = 0;
00258   if(!gSystem->AccessPathName(gSystem->ExpandPathName(GridFileName_Normal.c_str())))//if file exists
00259   {
00260     fnorm = new TFile(gSystem->ExpandPathName(GridFileName_Normal.c_str()),"READ");
00261     for(i=0;i<nBins;i++)
00262     {
00263       temp = (TTree*)fnorm->Get(Form("Bin_%i",i));
00264       temp->SetName(Form("Bin_%i_Normal",i));
00265       temp->SetBranchAddress("Background",&grid_background);
00266       temp->SetBranchAddress("Signal",&grid_signal);
00267       temp->SetBranchAddress("Delta",&grid_delta);
00268       temp->SetBranchAddress("Th13Axis",&grid_sinsq2th13);
00269       if(IncludeOscParErrs)
00270       {
00271         temp->SetBranchAddress("DNExp_DOscPars",&grid_bin_oscparerr[i]);
00272       }
00273       GridTree_Normal.push_back(temp);
00274       
00275       GridTree_2_Normal.push_back( vector<TTree*>() );
00276 
00277       if (i<nBinsFHC){
00278         for(j=0;j<ExtrapFHC.size();j++)
00279           {
00280             temp = (TTree*)fnorm->Get(Form("Bin_%i_Run_%i",i,j));
00281             temp->SetName(Form("Bin_%i_Run_%i_Normal",i,j));
00282             temp->SetBranchAddress("NC",&grid_nc);
00283             temp->SetBranchAddress("NuMuCC",&grid_numucc);
00284             temp->SetBranchAddress("BNueCC",&grid_bnuecc);
00285             temp->SetBranchAddress("NuTauCC",&grid_nutaucc);
00286             temp->SetBranchAddress("Signal",&grid_nue);
00287             temp->SetBranchAddress("Delta",&grid_delta);
00288             temp->SetBranchAddress("Th13Axis",&grid_sinsq2th13);
00289             GridTree_2_Normal[i].push_back(temp);
00290           }
00291       } else if (i>=nBinsFHC){
00292         for(j=0;j<ExtrapRHC.size();j++)
00293           {
00294             temp = (TTree*)fnorm->Get(Form("Bin_%i_Run_%i",i,j));
00295             temp->SetName(Form("Bin_%i_Run_%i_Normal",i,j));
00296             temp->SetBranchAddress("NC",&grid_nc);
00297             temp->SetBranchAddress("NuMuCC",&grid_numucc);
00298             temp->SetBranchAddress("BNueCC",&grid_bnuecc);
00299             temp->SetBranchAddress("NuTauCC",&grid_nutaucc);
00300             temp->SetBranchAddress("Signal",&grid_nue);
00301             temp->SetBranchAddress("Delta",&grid_delta);
00302             temp->SetBranchAddress("Th13Axis",&grid_sinsq2th13);
00303             GridTree_2_Normal[i].push_back(temp);
00304           }
00305       }
00306     }
00307     nPts_Normal = GridTree_Normal[0]->GetEntries();
00308     
00309     paramtree_Normal = (TTree*)fnorm->Get("paramtree");
00310     paramtree_Normal->SetName("paramtree_Normal");
00311     paramtree_Normal->SetBranchAddress("farPOTFHC",&fpfhc);
00312     paramtree_Normal->SetBranchAddress("farPOTRHC",&fprhc);
00313     paramtree_Normal->SetBranchAddress("Theta12",&grid_n_th12);
00314     paramtree_Normal->SetBranchAddress("Theta23",&grid_n_th23);
00315     paramtree_Normal->SetBranchAddress("DeltaMSq23",&grid_n_dm2_32);
00316     paramtree_Normal->SetBranchAddress("DeltaMSq12",&grid_n_dm2_21);
00317     paramtree_Normal->SetBranchAddress("Theta13",&grid_n_th13);
00318     paramtree_Normal->GetEntry(0);
00319     
00320     if(GridNorm>0)
00321     {
00322       GridScale_Normal_FHC = GridNorm/fpfhc;
00323       GridScale_Normal_RHC = GridNorm/fprhc;
00324     }
00325   }
00326   else
00327   {
00328     cout<<"Grid file (normal hierarchy) doesn't exist."<<endl;
00329     return;
00330   }
00331 
00332   TFile *finvt;
00333   GridTree_Inverted.clear();
00334   nPts_Inverted = 0;
00335   if(!gSystem->AccessPathName(gSystem->ExpandPathName(GridFileName_Inverted.c_str())))//if file exists
00336   {
00337     finvt = new TFile(gSystem->ExpandPathName(GridFileName_Inverted.c_str()),"READ");
00338     for(i=0;i<nBins;i++)
00339     {
00340       temp = (TTree*)finvt->Get(Form("Bin_%i",i));
00341       temp->SetName(Form("Bin_%i_Inverted",i));
00342       temp->SetBranchAddress("Background",&grid_background);
00343       temp->SetBranchAddress("Signal",&grid_signal);
00344       temp->SetBranchAddress("Delta",&grid_delta);
00345       temp->SetBranchAddress("Th13Axis",&grid_sinsq2th13);
00346       if(IncludeOscParErrs)
00347       {
00348         temp->SetBranchAddress("DNExp_DOscPars",&grid_bin_oscparerr[i]);
00349       }
00350       GridTree_Inverted.push_back(temp);
00351       
00352       GridTree_2_Inverted.push_back( vector<TTree*>() );
00353 
00354       if (i<nBinsFHC){      
00355         for(j=0;j<ExtrapFHC.size();j++)
00356           {
00357             temp = (TTree*)finvt->Get(Form("Bin_%i_Run_%i",i,j));
00358             temp->SetName(Form("Bin_%i_Run_%i_Inverted",i,j));
00359             temp->SetBranchAddress("NC",&grid_nc);
00360             temp->SetBranchAddress("NuMuCC",&grid_numucc);
00361             temp->SetBranchAddress("BNueCC",&grid_bnuecc);
00362             temp->SetBranchAddress("NuTauCC",&grid_nutaucc);
00363             temp->SetBranchAddress("Signal",&grid_nue);
00364             temp->SetBranchAddress("Delta",&grid_delta);
00365             temp->SetBranchAddress("Th13Axis",&grid_sinsq2th13);
00366             GridTree_2_Inverted[i].push_back(temp);
00367           }
00368       } else if (i>=nBinsFHC){
00369         for(j=0;j<ExtrapRHC.size();j++)
00370           {
00371             temp = (TTree*)finvt->Get(Form("Bin_%i_Run_%i",i,j));
00372             temp->SetName(Form("Bin_%i_Run_%i_Inverted",i,j));
00373             temp->SetBranchAddress("NC",&grid_nc);
00374             temp->SetBranchAddress("NuMuCC",&grid_numucc);
00375             temp->SetBranchAddress("BNueCC",&grid_bnuecc);
00376             temp->SetBranchAddress("NuTauCC",&grid_nutaucc);
00377             temp->SetBranchAddress("Signal",&grid_nue);
00378             temp->SetBranchAddress("Delta",&grid_delta);
00379             temp->SetBranchAddress("Th13Axis",&grid_sinsq2th13);
00380             GridTree_2_Inverted[i].push_back(temp);
00381           }
00382       }
00383     }
00384     nPts_Inverted = GridTree_Inverted[0]->GetEntries();
00385     
00386     paramtree_Inverted = (TTree*)finvt->Get("paramtree");
00387     paramtree_Inverted->SetName("paramtree_Inverted");
00388     paramtree_Inverted->SetBranchAddress("farPOTFHC",&fpfhc);
00389     paramtree_Inverted->SetBranchAddress("farPOTRHC",&fprhc);
00390     paramtree_Inverted->SetBranchAddress("Theta12",&grid_i_th12);
00391     paramtree_Inverted->SetBranchAddress("Theta23",&grid_i_th23);
00392     paramtree_Inverted->SetBranchAddress("DeltaMSq23",&grid_i_dm2_32);
00393     paramtree_Inverted->SetBranchAddress("DeltaMSq12",&grid_i_dm2_21);
00394     paramtree_Inverted->SetBranchAddress("Theta13",&grid_i_th13);
00395     paramtree_Inverted->GetEntry(0);
00396     
00397 
00398     if(GridNorm>0)
00399     {
00400       GridScale_Inverted_FHC = GridNorm/fpfhc;
00401       GridScale_Inverted_RHC = GridNorm/fprhc;
00402     }
00403   }
00404   else
00405   {
00406     cout<<"Grid file (inverted hierarchy) doesn't exist."<<endl;
00407     return;
00408   }
00409   
00410   return;
00411 }

void NueFit2D_Joint::RunDataGrid ( string  filenorm,
string  fileinvt 
) [virtual]

Reimplemented from NueFit2D.

Definition at line 484 of file NueFit2D_Joint.cxx.

References NueFit2D::BinLikelihood(), NueFit2D::Bkgd, DefineBinDlnLMinuit(), DefineStdDlnLMinuit(), NueFit2D::ErrCalc, NueFit2D::ErrorMatrix, ExtrapFHC, ExtrapRHC, NueFit2D::FitMethod, NueFit2D::FracErr_Bkgd, NueFit2D::FracErr_Sig, GetMinLikelihood(), NueFit2D::grid_background, NueFit2D::grid_bnuecc, NueFit2D::grid_delta, NueFit2D::grid_nc, NueFit2D::grid_nue, NueFit2D::grid_numucc, NueFit2D::grid_nutaucc, NueFit2D::grid_signal, NueFit2D::grid_sinsq2th13, GridScale_Inverted_FHC, GridScale_Inverted_RHC, GridScale_Normal_FHC, GridScale_Normal_RHC, NueFit2D::GridTree_2_Inverted, NueFit2D::GridTree_2_Normal, NueFit2D::GridTree_Inverted, NueFit2D::GridTree_Normal, gSystem(), NueFit2D::InvErrorMatrix, NueFit2D::nBins, nBinsFHC, nBinsRHC, NueFit2D::NObs, NueFit2D::nPts_Inverted, NueFit2D::nPts_Normal, NueFit2D::PoissonChi2(), ReadGridFiles(), NueFit2D::ScaledChi2(), ErrorCalc::SetGridPred(), ErrorCalc::SetUseGrid(), NueFit2D::Sig, NueFit2D::StandardChi2(), and NueFit2D::StandardLikelihood().

00485 {
00486   if(NObs==0)
00487   {
00488     cout<<"NObs not set.  Quitting..."<<endl;
00489     return;
00490   }
00491   if((ExtrapFHC.size()+ExtrapRHC.size())==0)
00492   {
00493     cout<<"No Extrapolate2D input.  Quitting..."<<endl;
00494     return;
00495   }
00496   for(unsigned int ie=0;ie<ExtrapFHC.size();ie++)
00497   {
00498     ExtrapFHC[ie]->GetPrediction();
00499   }
00500   for(unsigned int ie=0;ie<ExtrapRHC.size();ie++)
00501   {
00502     ExtrapRHC[ie]->GetPrediction();
00503   }
00504   if(FitMethod==1 && (FracErr_Bkgd==0 || FracErr_Sig==0))
00505   {
00506     cout<<"FracErr_Bkgd and FracErr_Sig need to be set for ScaledChi2.  Quitting..."<<endl;
00507     return;
00508   }
00509   if(ErrCalc==0)
00510   {
00511     cout<<"No ErrorCalc object set!  Quitting..."<<endl;
00512     return;
00513   }
00514   if(FitMethod==3) DefineStdDlnLMinuit();
00515   if(FitMethod==4) DefineBinDlnLMinuit();
00516   
00517   nBins = NObs->GetNbinsX();
00518   nBinsFHC=0;
00519   nBinsRHC=0;  
00520   if (ExtrapFHC.size()>0) nBinsFHC = ExtrapFHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
00521   if (ExtrapRHC.size()>0) nBinsRHC = ExtrapRHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
00522   if (nBins!=(nBinsFHC+nBinsRHC)) {
00523     cout << "ERROR: Bin mismatch between NObs and expected distribution.  Aborting..." << endl;
00524     return;
00525   }
00526 
00527   if(ErrorMatrix==0) ErrorMatrix = new TH2D("ErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
00528   if(InvErrorMatrix==0) InvErrorMatrix = new TH2D("InvErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
00529   
00530   ReadGridFiles();
00531   
00532   if(nPts_Normal==0 || nPts_Inverted==0) return;
00533   
00534   TH1D *nexp_bkgd = new TH1D("nexp_bkgd","",nBins,-0.5,nBins-0.5);
00535   TH1D *nexp_signal = new TH1D("nexp_signal","",nBins,-0.5,nBins-0.5);
00536   TH1D *nexp = new TH1D("nexp","",nBins,-0.5,nBins-0.5);
00537   
00538   int i;
00539   unsigned int j,k;
00540   
00541   double chi2data,chi2min;
00542   
00543   vector< vector<double> > nc,numucc,bnuecc,nutaucc,sig;
00544   for(j=0;j<nBins;j++)
00545   {
00546     nc.push_back( vector<double>() );
00547     numucc.push_back( vector<double>() );
00548     bnuecc.push_back( vector<double>() );
00549     nutaucc.push_back( vector<double>() );
00550     sig.push_back( vector<double>() );
00551 
00552   if (j<nBinsFHC){
00553       for(k=0;k<ExtrapFHC.size();k++)
00554         {
00555           nc[j].push_back(0);
00556           numucc[j].push_back(0);
00557           bnuecc[j].push_back(0);
00558           nutaucc[j].push_back(0);
00559           sig[j].push_back(0);
00560         }
00561     } else if (j>=nBinsFHC){
00562       for(k=0;k<ExtrapRHC.size();k++)
00563         {
00564           nc[j].push_back(0);
00565           numucc[j].push_back(0);
00566           bnuecc[j].push_back(0);
00567           nutaucc[j].push_back(0);
00568           sig[j].push_back(0);
00569         }
00570     }
00571   }
00572   
00573   Bkgd = (TH1D*)NObs->Clone("Bkgd");
00574   Bkgd->Reset();
00575   Sig = (TH1D*)NObs->Clone("Sig");
00576   Sig->Reset();
00577   
00578   //normal hierarchy
00579   
00580   ofstream myfile;
00581   myfile.open(gSystem->ExpandPathName(filenorm.c_str()));
00582   
00583   for(i=0;i<nPts_Normal;i++)
00584   {
00585     if(i%100==0) cout<<100.*i/nPts_Normal<<"% complete for normal hierarchy"<<endl;
00586     
00587     nexp_bkgd->Reset();
00588     nexp_signal->Reset();
00589     nexp->Reset();
00590     
00591     for(j=0;j<nBins;j++)
00592     {
00593       GridTree_Normal[j]->GetEntry(i);
00594       if (j<nBinsFHC){
00595         nexp_bkgd->SetBinContent(j+1,grid_background*GridScale_Normal_FHC);
00596         nexp_signal->SetBinContent(j+1,grid_signal*GridScale_Normal_FHC);
00597         for(k=0;k<ExtrapFHC.size();k++)
00598           {
00599             GridTree_2_Normal[j][k]->GetEntry(i);
00600             nc[j][k] = grid_nc*GridScale_Normal_FHC;
00601             numucc[j][k] = grid_numucc*GridScale_Normal_FHC;
00602             bnuecc[j][k] = grid_bnuecc*GridScale_Normal_FHC;
00603             nutaucc[j][k] = grid_nutaucc*GridScale_Normal_FHC;
00604             sig[j][k] = grid_nue*GridScale_Normal_FHC;
00605           }
00606       } else if (j>=nBinsFHC){
00607         nexp_bkgd->SetBinContent(j+1,grid_background*GridScale_Normal_RHC);
00608         nexp_signal->SetBinContent(j+1,grid_signal*GridScale_Normal_RHC);
00609         for(k=0;k<ExtrapRHC.size();k++)
00610           {
00611             GridTree_2_Normal[j][k]->GetEntry(i);
00612             nc[j][k] = grid_nc*GridScale_Normal_RHC;
00613             numucc[j][k] = grid_numucc*GridScale_Normal_RHC;
00614             bnuecc[j][k] = grid_bnuecc*GridScale_Normal_RHC;
00615             nutaucc[j][k] = grid_nutaucc*GridScale_Normal_RHC;
00616             sig[j][k] = grid_nue*GridScale_Normal_RHC;
00617           }
00618       }
00619     }
00620     nexp->Add(nexp_bkgd,nexp_signal,1,1);
00621     ErrCalc->SetGridPred(nBinsFHC+nBinsRHC,nc,numucc,bnuecc,nutaucc,sig);
00622     
00623     chi2min=GetMinLikelihood(grid_delta,true);
00624     
00625     ErrCalc->SetUseGrid(true);//use grid predictions set up above
00626     chi2data = 1e10;
00627     if(FitMethod==0)
00628     {
00629       chi2data = PoissonChi2(nexp) - chi2min;
00630     }
00631     else if(FitMethod==1)
00632     {
00633       chi2data = ScaledChi2(Bkgd,Sig) - chi2min;
00634     }
00635     else if(FitMethod==2)
00636     {
00637       chi2data = StandardChi2(nexp) - chi2min;
00638     }
00639     else if(FitMethod==3)
00640     {
00641       //Likelihood: "Standard" (N syst, N nuisance)
00642       //Calculate the likelihood (x2 for chi)
00643       Bkgd->Reset();
00644       Bkgd->Add(nexp_bkgd);
00645       Sig->Reset();
00646       Sig->Add(nexp_signal);
00647       chi2data = StandardLikelihood() - chi2min;
00648     }
00649     else if(FitMethod==4)
00650     {
00651       //Likelihood: Bin by Bin Calculation of Systematics
00652       //Calculate the likelihood (x2 for chi)
00653       Bkgd->Reset();
00654       Bkgd->Add(nexp_bkgd);
00655       Sig->Reset();
00656       Sig->Add(nexp_signal);
00657       chi2data = BinLikelihood() - chi2min;
00658     }
00659     else
00660     {
00661       cout<<"Error in RunMultiBinFC(): Unknown 'FitMethod'."<<endl;
00662     }
00663     
00664     myfile << grid_sinsq2th13 << " " << grid_delta << " " << nexp_signal->Integral() << " ";
00665     for(j=0;j<nBins;j++)
00666     {
00667       myfile << nexp_signal->GetBinContent(j+1) << " ";
00668     }
00669     myfile << chi2data << endl;
00670   }
00671   
00672   myfile.close();
00673   
00674   //inverted hierarchy
00675   
00676   myfile.open(gSystem->ExpandPathName(fileinvt.c_str()));
00677   
00678   for(i=0;i<nPts_Inverted;i++)
00679   {
00680     if(i%100==0) cout<<100.*i/nPts_Inverted<<"% complete for inverted hierarchy"<<endl;
00681     
00682     nexp_bkgd->Reset();
00683     nexp_signal->Reset();
00684     nexp->Reset();
00685     
00686     for(j=0;j<nBins;j++)
00687     {
00688       GridTree_Inverted[j]->GetEntry(i);
00689      if (j<nBinsFHC){
00690         nexp_bkgd->SetBinContent(j+1,grid_background*GridScale_Inverted_FHC);
00691         nexp_signal->SetBinContent(j+1,grid_signal*GridScale_Inverted_FHC);
00692         for(k=0;k<ExtrapFHC.size();k++)
00693           {
00694             GridTree_2_Inverted[j][k]->GetEntry(i);
00695             nc[j][k] = grid_nc*GridScale_Inverted_FHC;
00696             numucc[j][k] = grid_numucc*GridScale_Inverted_FHC;
00697             bnuecc[j][k] = grid_bnuecc*GridScale_Inverted_FHC;
00698             nutaucc[j][k] = grid_nutaucc*GridScale_Inverted_FHC;
00699             sig[j][k] = grid_nue*GridScale_Inverted_FHC;
00700           }
00701       } else if (j>=nBinsFHC){
00702         nexp_bkgd->SetBinContent(j+1,grid_background*GridScale_Inverted_RHC);
00703         nexp_signal->SetBinContent(j+1,grid_signal*GridScale_Inverted_RHC);
00704         for(k=0;k<ExtrapRHC.size();k++)
00705           {
00706             GridTree_2_Inverted[j][k]->GetEntry(i);
00707             nc[j][k] = grid_nc*GridScale_Inverted_RHC;
00708             numucc[j][k] = grid_numucc*GridScale_Inverted_RHC;
00709             bnuecc[j][k] = grid_bnuecc*GridScale_Inverted_RHC;
00710             nutaucc[j][k] = grid_nutaucc*GridScale_Inverted_RHC;
00711             sig[j][k] = grid_nue*GridScale_Inverted_RHC;
00712           }
00713       }
00714     }
00715     nexp->Add(nexp_bkgd,nexp_signal,1,1);
00716     ErrCalc->SetGridPred(nBins,nc,numucc,bnuecc,nutaucc,sig);
00717     
00718     chi2min=GetMinLikelihood(grid_delta,false);
00719     
00720     ErrCalc->SetUseGrid(true);//use grid predictions set up above
00721     chi2data = 1e10;
00722     if(FitMethod==0)
00723     {
00724       chi2data = PoissonChi2(nexp) - chi2min;
00725     }
00726     else if(FitMethod==1)
00727     {
00728       chi2data = ScaledChi2(Bkgd,Sig) - chi2min;
00729     }
00730     else if(FitMethod==2)
00731     {
00732       chi2data = StandardChi2(nexp) - chi2min;
00733     }
00734     else if(FitMethod==3)
00735     {
00736       //Likelihood: "Standard" (N syst, N nuisance)
00737       //Calculate the likelihood (x2 for chi)
00738       Bkgd->Reset();
00739       Bkgd->Add(nexp_bkgd);
00740       Sig->Reset();
00741       Sig->Add(nexp_signal);
00742       chi2data = StandardLikelihood() - chi2min;
00743     }
00744     else if(FitMethod==4)
00745     {
00746       //Likelihood: Bin by Bin Calculation of Systematics
00747       //Calculate the likelihood (x2 for chi)
00748       Bkgd->Reset();
00749       Bkgd->Add(nexp_bkgd);
00750       Sig->Reset();
00751       Sig->Add(nexp_signal);
00752       chi2data = BinLikelihood() - chi2min;
00753     }
00754     else
00755     {
00756       cout<<"Error in RunMultiBinFC(): Unknown 'FitMethod'."<<endl;
00757     }
00758     
00759     myfile << grid_sinsq2th13 << " " << grid_delta << " " << nexp_signal->Integral() << " ";
00760     for(j=0;j<nBins;j++)
00761     {
00762       myfile << nexp_signal->GetBinContent(j+1) << " ";
00763     }
00764     myfile << chi2data << endl;
00765   }
00766   
00767   myfile.close();
00768   
00769   return;
00770 }

void NueFit2D_Joint::RunDeltaChi2Contour ( int  cl = 0  )  [virtual]

Reimplemented from NueFit2D.

Definition at line 1467 of file NueFit2D_Joint.cxx.

References NueFit2D::BinLikelihood(), NueFit2D::Bkgd, CombineFHCRHC(), DefineBinDlnLMinuit(), DefineStdDlnLMinuit(), NueFit2D::delta, NueFit2D::DeltaHigh, NueFit2D::DeltaLow, NueFit2D::ErrorMatrix, ExtrapFHC, ExtrapRHC, NueFit2D::FitMethod, NueFit2D::FracErr_Bkgd, NueFit2D::FracErr_Sig, gSystem(), NueFit2D::InvErrorMatrix, NueFit2D::nBins, NueFit2D::nDeltaSteps, NueFit2D::NExp, NueFit2D::NObs, NueFit2D::nSinSq2Th13Steps, NueFit2D::outFileName, NueFit2D::PoissonChi2(), NueFit2D::ScaledChi2(), NueFit2D::Sig, NueFit2D::SinSq2Th13High, NueFit2D::SinSq2Th13Low, NueFit2D::StandardChi2(), and NueFit2D::StandardLikelihood().

01468 {
01469   if(NObs==0)
01470   {
01471     cout<<"NObs not set.  Quitting..."<<endl;
01472     return;
01473   }
01474   if(FitMethod==1 && (FracErr_Bkgd==0 || FracErr_Sig==0))
01475   {
01476     cout<<"FracErr_Bkgd and FracErr_Sig need to be set for ScaledChi2.  Quitting..."<<endl;
01477     return;
01478   }
01479   if((ExtrapFHC.size()+ExtrapRHC.size())==0)
01480   {
01481     cout<<"No Extrapolate2D input.  Quitting..."<<endl;
01482     return;
01483   }
01484   
01485   if(FitMethod==3) DefineStdDlnLMinuit();
01486   if(FitMethod==4) DefineBinDlnLMinuit();
01487   
01488   Bkgd = (TH1D*)NObs->Clone("Bkgd");
01489   Bkgd->Reset();
01490   Sig = (TH1D*)NObs->Clone("Sig");
01491   Sig->Reset();
01492   NExp = (TH1D*)NObs->Clone("NExp");
01493   NExp->Reset();
01494   
01495   unsigned int ie;
01496   for(ie=0;ie<ExtrapFHC.size();ie++) ExtrapFHC[ie]->GetPrediction();
01497   for(ie=0;ie<ExtrapRHC.size();ie++) ExtrapRHC[ie]->GetPrediction();
01498   
01499   if(ErrorMatrix==0) ErrorMatrix = new TH2D("ErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
01500   if(InvErrorMatrix==0) InvErrorMatrix = new TH2D("InvErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
01501   
01502   Int_t i,j;
01503   
01504   Double_t ss2th13 = 0;
01505   Double_t delta = 0;
01506   Double_t Deltaincrement = 0;
01507   if(nDeltaSteps>0) Deltaincrement = (DeltaHigh - DeltaLow)/(nDeltaSteps);
01508   Double_t Th13increment = 0;
01509   if(nSinSq2Th13Steps>0) Th13increment = (SinSq2Th13High - SinSq2Th13Low)/(nSinSq2Th13Steps);
01510   
01511   double chi2[3],val[3];
01512   TGraph *g3;
01513   TF1* fit;
01514   double best;
01515   double minchi2;
01516   
01517   double limit;
01518   double delchi2;
01519   double sprev,delchi2prev;
01520   double contourlvl = 0;
01521   if(cl==0)//90% CL
01522   {
01523     contourlvl = 2.71;
01524   }
01525   else if(cl==1)//68.3% cL
01526   {
01527     contourlvl = 1.0;
01528   }
01529   else
01530   {
01531     cout<<"Error in RunDeltaChi2Contour(): Input value should be 0 or 1 for 90% or 68.3%.  Quitting..."<<endl;
01532     return;
01533   }
01534   
01535   cout<<"Seeking ";
01536   if(cl==0) cout<<"90% ";
01537   else cout<<"68% ";
01538   cout<<" CL upper limit"<<endl;
01539   
01540   vector<double> deltapts;
01541   vector<double> bestfit_norm;
01542   vector<double> bestfit_invt;
01543   vector<double> limit_norm;
01544   vector<double> limit_invt;
01545   
01546   for(j=0;j<nDeltaSteps+1;j++)
01547   {
01548     delta = j*Deltaincrement*TMath::Pi() + DeltaLow;
01549     for(ie=0;ie<ExtrapFHC.size();ie++) ExtrapFHC[ie]->SetDeltaCP(delta);
01550     for(ie=0;ie<ExtrapRHC.size();ie++) ExtrapRHC[ie]->SetDeltaCP(delta);
01551 
01552     deltapts.push_back(delta/TMath::Pi());
01553     
01554     best=-1.;
01555     minchi2=100;
01556     for(i=0;i<3;i++)
01557     {
01558       chi2[i]=-1.;
01559       val[i]=-1.;
01560     }
01561     for(i=0;i<nSinSq2Th13Steps+1;i++)
01562     {
01563       chi2[0] = chi2[1];
01564       chi2[1] = chi2[2];
01565       val[0] = val[1];
01566       val[1] = val[2];
01567       ss2th13 = i*Th13increment + SinSq2Th13Low;
01568 
01569       for(ie=0;ie<ExtrapFHC.size();ie++)
01570       {
01571         ExtrapFHC[ie]->SetSinSq2Th13(ss2th13);
01572         ExtrapFHC[ie]->OscillatePrediction();
01573       }
01574       for(ie=0;ie<ExtrapRHC.size();ie++)
01575       {
01576         ExtrapRHC[ie]->SetSinSq2Th13(ss2th13);
01577         ExtrapRHC[ie]->OscillatePrediction();
01578       }
01579 
01580       Bkgd->Reset();
01581       Sig->Reset();
01582       NExp->Reset();
01583 
01584       CombineFHCRHC();
01585       NExp->Add(Bkgd);
01586       NExp->Add(Sig);
01587       
01588       chi2[2] = 1e10;
01589       if(FitMethod==0)
01590       {
01591         chi2[2] = PoissonChi2(NExp);
01592       }
01593       else if(FitMethod==1)
01594       {
01595         chi2[2] = ScaledChi2(Bkgd,Sig);
01596       }
01597       else if(FitMethod==2)
01598       {
01599         chi2[2] = StandardChi2(NExp);
01600       }
01601       else if(FitMethod==3)
01602       {
01603         //Likelihood: "Standard" (N syst, N nuisance)
01604         //Calculate the likelihood (x2 for chi)
01605         chi2[2] = StandardLikelihood();
01606       }
01607       else if(FitMethod==4)
01608       {
01609         //Likelihood: Bin by Bin Calculation of Systematics
01610         //Calculate the likelihood (x2 for chi)
01611         chi2[2] = BinLikelihood();
01612       }
01613       else
01614       {
01615         cout<<"Error in RunDeltaChi2Contour(): Unknown 'FitMethod'."<<endl;
01616       }
01617       
01618       val[2] = ss2th13;
01619       
01620       if(i<2) continue;
01621       
01622       if(i<3 && chi2[2]>chi2[1] && chi2[1]>chi2[0])//first three points are increasing, first point is minimum.
01623       {
01624         best = val[0];
01625         minchi2 = chi2[0];
01626         cout<<"minimum at 1st point: "<<minchi2<<" at "<<best<<endl;
01627         break;
01628       }
01629       
01630       if(chi2[2]>chi2[1] && chi2[0]>chi2[1])//found minimum
01631       {
01632         g3 = new TGraph(3, val, chi2);
01633         fit = new TF1("pol2", "pol2");
01634         g3->Fit(fit, "Q");//fit to second order polynominal
01635         if(fit->GetParameter(2) > 0)//if the x^2 term is nonzero
01636         {
01637           best = -fit->GetParameter(1)/(2*fit->GetParameter(2));//the location of the minimum is -p1/(2*p2)
01638           minchi2 = fit->GetParameter(0) + fit->GetParameter(1)*best + fit->GetParameter(2)*best*best;
01639           cout<<"minimum with fit: "<<minchi2<<" at "<<best<<endl;
01640         }
01641         else//if the x^2 term is zero, then just use the minimum you got by scanning
01642         {
01643           best = val[1];
01644           minchi2 = chi2[1];
01645           cout<<"minimum with scan: "<<minchi2<<" at "<<best<<endl;
01646         }
01647         break;
01648       }
01649     }
01650     bestfit_norm.push_back(best);
01651     
01652     limit = -1.;
01653     delchi2prev = 1000;
01654     sprev = 0;
01655     for(i=0;i<nSinSq2Th13Steps+1;i++)
01656     {
01657       ss2th13 = i*Th13increment + SinSq2Th13Low;
01658 
01659       for(ie=0;ie<ExtrapFHC.size();ie++)
01660       {
01661         ExtrapFHC[ie]->SetSinSq2Th13(ss2th13);
01662         ExtrapFHC[ie]->OscillatePrediction();
01663       }
01664       for(ie=0;ie<ExtrapRHC.size();ie++)
01665       {
01666         ExtrapRHC[ie]->SetSinSq2Th13(ss2th13);
01667         ExtrapRHC[ie]->OscillatePrediction();
01668       }
01669 
01670       Bkgd->Reset();
01671       Sig->Reset();
01672       NExp->Reset();
01673 
01674       CombineFHCRHC();
01675       NExp->Add(Bkgd);
01676       NExp->Add(Sig);
01677       
01678       delchi2 = 1e10;
01679       if(FitMethod==0)
01680       {
01681         delchi2 = PoissonChi2(NExp) - minchi2;
01682       }
01683       else if(FitMethod==1)
01684       {
01685         delchi2 = ScaledChi2(Bkgd,Sig) - minchi2;
01686       }
01687       else if(FitMethod==2)
01688       {
01689         delchi2 = StandardChi2(NExp) - minchi2;
01690       }
01691       else if(FitMethod==3)
01692       {
01693         //Likelihood: "Standard" (N syst, N nuisance)
01694         //Calculate the likelihood (x2 for chi)
01695         delchi2 = StandardLikelihood() - minchi2;
01696       }
01697       else if(FitMethod==4)
01698       {
01699         //Likelihood: Bin by Bin Calculation of Systematics
01700         //Calculate the likelihood (x2 for chi)
01701         delchi2 = BinLikelihood() - minchi2;
01702       }
01703       else
01704       {
01705         cout<<"Error in RunDeltaChi2Contour(): Unknown 'FitMethod'."<<endl;
01706       }
01707       
01708       if(i==1) continue;
01709       
01710       if(delchi2>contourlvl && delchi2prev<contourlvl)
01711       {
01712         limit = ss2th13 + ((ss2th13-sprev)/(delchi2 - delchi2prev))*(contourlvl - delchi2);
01713         cout<<delchi2prev<<", "<<sprev<<", "<<delchi2<<", "<<ss2th13<<", "<<limit<<endl;
01714         break;
01715       }
01716       delchi2prev = delchi2;
01717       sprev = ss2th13;
01718     }
01719     limit_norm.push_back(limit);
01720   }
01721   
01722   for(ie=0;ie<ExtrapFHC.size();ie++) ExtrapFHC[ie]->InvertMassHierarchy();
01723   for(ie=0;ie<ExtrapRHC.size();ie++) ExtrapRHC[ie]->InvertMassHierarchy();
01724 
01725   
01726   for(j=0;j<nDeltaSteps+1;j++)
01727   {
01728     delta = j*Deltaincrement*TMath::Pi() + DeltaLow;
01729     for(ie=0;ie<ExtrapFHC.size();ie++) ExtrapFHC[ie]->SetDeltaCP(delta);
01730     for(ie=0;ie<ExtrapRHC.size();ie++) ExtrapRHC[ie]->SetDeltaCP(delta);
01731     
01732     best=-1.;
01733     minchi2=100;
01734     for(i=0;i<3;i++)
01735     {
01736       chi2[i]=-1.;
01737       val[i]=-1.;
01738     }
01739     for(i=0;i<nSinSq2Th13Steps+1;i++)
01740     {
01741       chi2[0] = chi2[1];
01742       chi2[1] = chi2[2];
01743       val[0] = val[1];
01744       val[1] = val[2];
01745       ss2th13 = i*Th13increment + SinSq2Th13Low;
01746 
01747       for(ie=0;ie<ExtrapFHC.size();ie++)
01748       {
01749         ExtrapFHC[ie]->SetSinSq2Th13(ss2th13);
01750         ExtrapFHC[ie]->OscillatePrediction();
01751       }
01752 
01753       for(ie=0;ie<ExtrapRHC.size();ie++)
01754       {
01755         ExtrapRHC[ie]->SetSinSq2Th13(ss2th13);
01756         ExtrapRHC[ie]->OscillatePrediction();
01757       }
01758 
01759       Bkgd->Reset();
01760       Sig->Reset();
01761       NExp->Reset();
01762 
01763       CombineFHCRHC();
01764       NExp->Add(Bkgd);
01765       NExp->Add(Sig);
01766       
01767       chi2[2] = 1e10;
01768       if(FitMethod==0)
01769       {
01770         chi2[2] = PoissonChi2(NExp);
01771       }
01772       else if(FitMethod==1)
01773       {
01774         chi2[2] = ScaledChi2(Bkgd,Sig);
01775       }
01776       else if(FitMethod==2)
01777       {
01778         chi2[2] = StandardChi2(NExp);
01779       }
01780       else if(FitMethod==3)
01781       {
01782         //Likelihood: "Standard" (N syst, N nuisance)
01783         //Calculate the likelihood (x2 for chi)
01784         chi2[2] = StandardLikelihood();
01785       }
01786       else if(FitMethod==4)
01787       {
01788         //Likelihood: Bin by Bin Calculation of Systematics
01789         //Calculate the likelihood (x2 for chi)
01790         chi2[2] = BinLikelihood();
01791       }
01792       else
01793       {
01794         cout<<"Error in RunDeltaChi2Contour(): Unknown 'FitMethod'."<<endl;
01795       }
01796       
01797       val[2] = ss2th13;
01798       
01799       if(i<2) continue;
01800       
01801       if(i<3 && chi2[2]>chi2[1] && chi2[1]>chi2[0])//first three points are increasing, first point is minimum.
01802       {
01803         best = val[0];
01804         minchi2 = chi2[0];
01805         cout<<"minimum at 1st point: "<<minchi2<<" at "<<best<<endl;
01806         break;
01807       }
01808       
01809       if(chi2[2]>chi2[1] && chi2[0]>chi2[1])//found minimum
01810       {
01811         g3 = new TGraph(3, val, chi2);
01812         fit = new TF1("pol2", "pol2");
01813         g3->Fit(fit, "Q");//fit to second order polynominal
01814         if(fit->GetParameter(2) > 0)//if the x^2 term is nonzero
01815         {
01816           best = -fit->GetParameter(1)/(2*fit->GetParameter(2));//the location of the minimum is -p1/(2*p2)
01817           minchi2 = fit->GetParameter(0) + fit->GetParameter(1)*best + fit->GetParameter(2)*best*best;
01818           cout<<"minimum with fit: "<<minchi2<<" at "<<best<<endl;
01819         }
01820         else//if the x^2 term is zero, then just use the minimum you got by scanning
01821         {
01822           best = val[1];
01823           minchi2 = chi2[1];
01824           cout<<"minimum with scan: "<<minchi2<<" at "<<best<<endl;
01825         }
01826         break;
01827       }
01828     }
01829     bestfit_invt.push_back(best);
01830     
01831     limit = -1.;
01832     delchi2prev = 1000;
01833     sprev = 0;
01834     for(i=0;i<nSinSq2Th13Steps+1;i++)
01835     {
01836       ss2th13 = i*Th13increment + SinSq2Th13Low;
01837 
01838       for(ie=0;ie<ExtrapFHC.size();ie++)
01839       {
01840         ExtrapFHC[ie]->SetSinSq2Th13(ss2th13);
01841         ExtrapFHC[ie]->OscillatePrediction();
01842       }
01843       for(ie=0;ie<ExtrapRHC.size();ie++)
01844       {
01845         ExtrapRHC[ie]->SetSinSq2Th13(ss2th13);
01846         ExtrapRHC[ie]->OscillatePrediction();
01847       }
01848 
01849       Bkgd->Reset();
01850       Sig->Reset();
01851       NExp->Reset();
01852 
01853       CombineFHCRHC();
01854       NExp->Add(Bkgd);
01855       NExp->Add(Sig);
01856       
01857       delchi2 = 1e10;
01858       if(FitMethod==0)
01859       {
01860         delchi2 = PoissonChi2(NExp) - minchi2;
01861       }
01862       else if(FitMethod==1)
01863       {
01864         delchi2 = ScaledChi2(Bkgd,Sig) - minchi2;
01865       }
01866       else if(FitMethod==2)
01867       {
01868         delchi2 = StandardChi2(NExp) - minchi2;
01869       }
01870       else if(FitMethod==3)
01871       {
01872         //Likelihood: "Standard" (N syst, N nuisance)
01873         //Calculate the likelihood (x2 for chi)
01874         delchi2 = StandardLikelihood() - minchi2;
01875       }
01876       else if(FitMethod==4)
01877       {
01878         //Likelihood: Bin by Bin Calculation of Systematics
01879         //Calculate the likelihood (x2 for chi)
01880         delchi2 = BinLikelihood() - minchi2;
01881       }
01882       else
01883       {
01884         cout<<"Error in RunDeltaChi2Contour(): Unknown 'FitMethod'."<<endl;
01885       }
01886       
01887       if(i==1) continue;
01888       
01889       if(delchi2>contourlvl && delchi2prev<contourlvl)
01890       {
01891         limit = ss2th13 + ((ss2th13-sprev)/(delchi2 - delchi2prev))*(contourlvl - delchi2);
01892         cout<<delchi2prev<<", "<<sprev<<", "<<delchi2<<", "<<ss2th13<<", "<<limit<<endl;
01893         break;
01894       }
01895       delchi2prev = delchi2;
01896       sprev = ss2th13;
01897     }
01898     limit_invt.push_back(limit);
01899   }
01900   
01901   double *s_best_n = new double[nDeltaSteps+1];
01902   double *s_best_i = new double[nDeltaSteps+1];
01903   double *s_limit_n = new double[nDeltaSteps+1];
01904   double *s_limit_i = new double[nDeltaSteps+1];
01905   double *d = new double[nDeltaSteps+1];
01906   for(i=0;i<nDeltaSteps+1;i++)
01907   {
01908     d[i] = deltapts.at(i);
01909     s_best_n[i] = bestfit_norm.at(i);
01910     s_best_i[i] = bestfit_invt.at(i);
01911     s_limit_n[i] = limit_norm.at(i);
01912     s_limit_i[i] = limit_invt.at(i);
01913   }
01914   
01915   TGraph *gn = new TGraph(nDeltaSteps+1,s_best_n,d);
01916   gn->SetMarkerStyle(20);
01917   gn->SetTitle("");
01918   gn->GetYaxis()->SetTitle("#delta");
01919   gn->GetXaxis()->SetTitle("sin^{2}2#theta_{13}");
01920   gn->GetXaxis()->SetLimits(0,0.6);
01921   gn->SetLineWidth(4);
01922   gn->SetMaximum(2);
01923   gn->SetName("BestFit_Normal");
01924   
01925   TGraph *gi = new TGraph(nDeltaSteps+1,s_best_i,d);
01926   gi->SetMarkerStyle(20);
01927   gi->SetTitle("");
01928   gi->GetYaxis()->SetTitle("#delta");
01929   gi->GetXaxis()->SetTitle("sin^{2}2#theta_{13}");
01930   gi->GetXaxis()->SetLimits(0,0.6);
01931   gi->SetLineWidth(4);
01932   gi->SetLineStyle(2);
01933   gi->SetMaximum(2);
01934   gi->SetName("BestFit_Inverted");
01935   
01936   TGraph *gn_limit = new TGraph(nDeltaSteps+1,s_limit_n,d);
01937   gn_limit->SetMarkerStyle(20);
01938   gn_limit->SetTitle("");
01939   gn_limit->GetYaxis()->SetTitle("#delta");
01940   gn_limit->GetXaxis()->SetTitle("sin^{2}2#theta_{13}");
01941   gn_limit->GetXaxis()->SetLimits(0,0.6);
01942   gn_limit->SetLineWidth(4);
01943   gn_limit->SetLineColor(kBlue);
01944   gn_limit->SetMarkerColor(kBlue);
01945   gn_limit->SetMaximum(2);
01946   gn_limit->SetName("Limit_Normal");
01947   
01948   TGraph *gi_limit = new TGraph(nDeltaSteps+1,s_limit_i,d);
01949   gi_limit->SetMarkerStyle(20);
01950   gi_limit->SetTitle("");
01951   gi_limit->GetYaxis()->SetTitle("#delta");
01952   gi_limit->GetXaxis()->SetTitle("sin^{2}2#theta_{13}");
01953   gi_limit->GetXaxis()->SetLimits(0,0.6);
01954   gi_limit->SetLineWidth(4);
01955   gi_limit->SetLineColor(kRed);
01956   gi_limit->SetMarkerColor(kRed);
01957   gi_limit->SetMaximum(2);
01958   gi_limit->SetName("Limit_Inverted");
01959   
01960   if(cl==0) cout<<"90% ";
01961   if(cl==1) cout<<"68% ";
01962   cout<<"confidence level limit = "<<limit_norm.at(0)<<", "<<limit_invt.at(0)<<endl;
01963   
01964   TFile *fout = new TFile(gSystem->ExpandPathName(outFileName.c_str()),"RECREATE");
01965   gn->Write();
01966   gi->Write();
01967   gn_limit->Write();
01968   gi_limit->Write();
01969   fout->Write();
01970   fout->Close();
01971   
01972   delete [] s_best_n;
01973   delete [] s_best_i;
01974   delete [] s_limit_n;
01975   delete [] s_limit_i;
01976   delete [] d;
01977   
01978   return;
01979 }

void NueFit2D_Joint::RunMultiBinFC_MHDeltaFit (  )  [virtual]

Reimplemented from NueFit2D.

Definition at line 3186 of file NueFit2D_Joint.cxx.

References NueFit2D::BinLikelihood(), NueFit2D::Bkgd, DefineBinDlnLMinuit(), DefineStdDlnLMinuit(), NueFit2D::ErrCalc, NueFit2D::ErrorMatrix, ExtrapFHC, ExtrapRHC, NueFit2D::FitMethod, Form(), NueFit2D::FracErr_Bkgd, NueFit2D::FracErr_Sig, GetMinLikelihood_Delta(), NueFit2D::grid_background, NueFit2D::grid_bnuecc, NueFit2D::grid_delta, NueFit2D::grid_nc, NueFit2D::grid_nue, NueFit2D::grid_numucc, NueFit2D::grid_nutaucc, NueFit2D::grid_signal, GridScale_Inverted_FHC, GridScale_Inverted_RHC, GridScale_Normal_FHC, GridScale_Normal_RHC, NueFit2D::GridTree_2_Inverted, NueFit2D::GridTree_2_Normal, NueFit2D::GridTree_Inverted, NueFit2D::GridTree_Normal, gSystem(), NueFit2D::InvErrorMatrix, NueFit2D::nBins, nBinsFHC, nBinsRHC, NueFit2D::NObs, NueFit2D::nPts_Inverted, NueFit2D::nPts_Normal, NueFit2D::outFileName, NueFit2D::PoissonChi2(), NueFit2D::PseudoExpFile, ReadGridFiles(), NueFit2D::ScaledChi2(), ErrorCalc::SetGridPred(), NueFit2D::SetupChi2Hists(), ErrorCalc::SetUseGrid(), NueFit2D::Sig, NueFit2D::StandardChi2(), and NueFit2D::StandardLikelihood().

03187 {
03188   if(NObs==0)
03189   {
03190     cout<<"NObs not set.  Quitting..."<<endl;
03191     return;
03192   }
03193   if((ExtrapFHC.size()+ExtrapRHC.size())==0)
03194   {
03195     cout<<"No Extrapolate2D input.  Quitting..."<<endl;
03196     return;
03197   }
03198   for(unsigned int ie=0;ie<ExtrapFHC.size();ie++)
03199   {
03200     ExtrapFHC[ie]->GetPrediction();
03201   }
03202   for(unsigned int ie=0;ie<ExtrapRHC.size();ie++)
03203   {
03204     ExtrapRHC[ie]->GetPrediction();
03205   }
03206   if(FitMethod==1 && (FracErr_Bkgd==0 || FracErr_Sig==0))
03207   {
03208     cout<<"FracErr_Bkgd and FracErr_Sig need to be set for ScaledChi2.  Quitting..."<<endl;
03209     return;
03210   }
03211   if(ErrCalc==0)
03212   {
03213     cout<<"No ErrorCalc object set!  Quitting..."<<endl;
03214     return;
03215   }
03216   if(FitMethod==3) DefineStdDlnLMinuit();
03217   if(FitMethod==4) DefineBinDlnLMinuit();
03218   
03219   nBinsFHC=0;
03220   nBinsRHC=0;  
03221   if (ExtrapFHC.size()>0) nBinsFHC = ExtrapFHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
03222   if (ExtrapRHC.size()>0) nBinsRHC = ExtrapRHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
03223   nBins = nBinsFHC + nBinsRHC;
03224   
03225   if(ErrorMatrix==0) ErrorMatrix = new TH2D("ErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
03226   if(InvErrorMatrix==0) InvErrorMatrix = new TH2D("InvErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
03227   
03228   ReadGridFiles();
03229   SetupChi2Hists();
03230   
03231   if(nPts_Normal==0 || nPts_Inverted==0) return;
03232   
03233   TH1D *nexp_bkgd = new TH1D("nexp_bkgd","",nBins,-0.5,nBins-0.5);
03234   TH1D *nexp_signal = new TH1D("nexp_signal","",nBins,-0.5,nBins-0.5);
03235   TH1D *nexp = new TH1D("nexp","",nBins,-0.5,nBins-0.5);
03236   
03237   int i;
03238   unsigned int j,k;
03239   
03240   if(gSystem->AccessPathName(gSystem->ExpandPathName(PseudoExpFile.c_str())))
03241   {
03242     cout<<"Pseudo-experiment file doesn't exist."<<endl;
03243     return;
03244   }
03245   
03246   TFile *f = new TFile(gSystem->ExpandPathName(PseudoExpFile.c_str()),"READ");
03247   
03248   double *chi2data_NH = new double[nPts_Normal];
03249   double *chi2data_IH = new double[nPts_Inverted];
03250   double *chi2min_NH = new double[nPts_Normal];
03251   double *chi2min_IH = new double[nPts_Inverted];
03252   
03253   int chi2databin;
03254   TH1D *chi2hist = (TH1D*)f->Get("Chi2Hist_MHFit_Normal_0");
03255   double chi2binwidth = chi2hist->GetBinWidth(1);
03256   double chi2start = chi2hist->GetXaxis()->GetBinLowEdge(1);
03257   double frac;
03258   delete chi2hist;
03259   
03260   double *deltapts = new double[nPts_Normal];
03261   
03262   vector< vector<double> > nc,numucc,bnuecc,nutaucc,sig;
03263   for(j=0;j<nBins;j++)
03264   {
03265     nc.push_back( vector<double>() );
03266     numucc.push_back( vector<double>() );
03267     bnuecc.push_back( vector<double>() );
03268     nutaucc.push_back( vector<double>() );
03269     sig.push_back( vector<double>() );
03270     
03271     if (j<nBinsFHC){
03272       for(k=0;k<ExtrapFHC.size();k++)
03273       {
03274         nc[j].push_back(0);
03275         numucc[j].push_back(0);
03276         bnuecc[j].push_back(0);
03277         nutaucc[j].push_back(0);
03278         sig[j].push_back(0);
03279       }
03280     } else if (j>=nBinsFHC){
03281       for(k=0;k<ExtrapRHC.size();k++)
03282       {
03283         nc[j].push_back(0);
03284         numucc[j].push_back(0);
03285         bnuecc[j].push_back(0);
03286         nutaucc[j].push_back(0);
03287         sig[j].push_back(0);
03288       }
03289     }
03290   }
03291   
03292   Bkgd = (TH1D*)NObs->Clone("Bkgd");
03293   Bkgd->Reset();
03294   Sig = (TH1D*)NObs->Clone("Sig");
03295   Sig->Reset();
03296   
03297   //normal hierarchy
03298   
03299   for(i=0;i<nPts_Normal;i++)
03300   {
03301     if(i%100==0) cout<<100.*i/nPts_Normal<<"% complete for normal hierarchy"<<endl;
03302     
03303     nexp_bkgd->Reset();
03304     nexp_signal->Reset();
03305     nexp->Reset();
03306     
03307     for(j=0;j<nBins;j++)
03308     {
03309       GridTree_Normal[j]->GetEntry(i);
03310       if (j<nBinsFHC){
03311         nexp_bkgd->SetBinContent(j+1,grid_background*GridScale_Normal_FHC);
03312         nexp_signal->SetBinContent(j+1,grid_signal*GridScale_Normal_FHC);
03313         for(k=0;k<ExtrapFHC.size();k++)
03314         {
03315           GridTree_2_Normal[j][k]->GetEntry(i);
03316           nc[j][k] = grid_nc*GridScale_Normal_FHC;
03317           numucc[j][k] = grid_numucc*GridScale_Normal_FHC;
03318           bnuecc[j][k] = grid_bnuecc*GridScale_Normal_FHC;
03319           nutaucc[j][k] = grid_nutaucc*GridScale_Normal_FHC;
03320           sig[j][k] = grid_nue*GridScale_Normal_FHC;
03321         }
03322       } else if (j>=nBinsFHC){
03323         nexp_bkgd->SetBinContent(j+1,grid_background*GridScale_Normal_RHC);
03324         nexp_signal->SetBinContent(j+1,grid_signal*GridScale_Normal_RHC);
03325         for(k=0;k<ExtrapRHC.size();k++)
03326         {
03327           GridTree_2_Normal[j][k]->GetEntry(i);
03328           nc[j][k] = grid_nc*GridScale_Normal_RHC;
03329           numucc[j][k] = grid_numucc*GridScale_Normal_RHC;
03330           bnuecc[j][k] = grid_bnuecc*GridScale_Normal_RHC;
03331           nutaucc[j][k] = grid_nutaucc*GridScale_Normal_RHC;
03332           sig[j][k] = grid_nue*GridScale_Normal_RHC;
03333         }
03334       }
03335     }
03336     nexp->Add(nexp_bkgd,nexp_signal,1,1);
03337     ErrCalc->SetGridPred(nBins,nc,numucc,bnuecc,nutaucc,sig);
03338     
03339     chi2min_NH[i]=GetMinLikelihood_Delta(true);
03340     
03341     ErrCalc->SetUseGrid(true);//use grid predictions set up above
03342     chi2data_NH[i] = 1e10;
03343     if(FitMethod==0)
03344     {
03345       chi2data_NH[i] = PoissonChi2(nexp);
03346     }
03347     else if(FitMethod==1)
03348     {
03349       chi2data_NH[i] = ScaledChi2(Bkgd,Sig);
03350     }
03351     else if(FitMethod==2)
03352     {
03353       chi2data_NH[i] = StandardChi2(nexp);
03354     }
03355     else if(FitMethod==3)
03356     {
03357       //Likelihood: "Standard" (N syst, N nuisance)
03358       //Calculate the likelihood (x2 for chi)
03359       Bkgd->Reset();
03360       Bkgd->Add(nexp_bkgd);
03361       Sig->Reset();
03362       Sig->Add(nexp_signal);
03363       chi2data_NH[i] = StandardLikelihood();
03364     }
03365     else if(FitMethod==4)
03366     {
03367       //Likelihood: Bin by Bin Calculation of Systematics
03368       //Calculate the likelihood (x2 for chi)
03369       Bkgd->Reset();
03370       Bkgd->Add(nexp_bkgd);
03371       Sig->Reset();
03372       Sig->Add(nexp_signal);
03373       chi2data_NH[i] = BinLikelihood();
03374     }
03375     else
03376     {
03377       cout<<"Error in RunMultiBinFC_MHDeltaFit(): Unknown 'FitMethod'."<<endl;
03378     }
03379     
03380     deltapts[i] = grid_delta/TMath::Pi();
03381   }
03382   
03383   //inverted hierarchy
03384   
03385   for(i=0;i<nPts_Inverted;i++)
03386   {
03387     if(i%100==0) cout<<100.*i/nPts_Inverted<<"% complete for inverted hierarchy"<<endl;
03388     
03389     nexp_bkgd->Reset();
03390     nexp_signal->Reset();
03391     nexp->Reset();
03392     
03393     for(j=0;j<nBins;j++)
03394     {
03395       GridTree_Inverted[j]->GetEntry(i);
03396       if (j<nBinsFHC){
03397         nexp_bkgd->SetBinContent(j+1,grid_background*GridScale_Inverted_FHC);
03398         nexp_signal->SetBinContent(j+1,grid_signal*GridScale_Inverted_FHC);
03399         for(k=0;k<ExtrapFHC.size();k++)
03400         {
03401           GridTree_2_Inverted[j][k]->GetEntry(i);
03402           nc[j][k] = grid_nc*GridScale_Inverted_FHC;
03403           numucc[j][k] = grid_numucc*GridScale_Inverted_FHC;
03404           bnuecc[j][k] = grid_bnuecc*GridScale_Inverted_FHC;
03405           nutaucc[j][k] = grid_nutaucc*GridScale_Inverted_FHC;
03406           sig[j][k] = grid_nue*GridScale_Inverted_FHC;
03407         }
03408       } else if (j>=nBinsFHC){
03409         nexp_bkgd->SetBinContent(j+1,grid_background*GridScale_Inverted_RHC);
03410         nexp_signal->SetBinContent(j+1,grid_signal*GridScale_Inverted_RHC);
03411         for(k=0;k<ExtrapRHC.size();k++)
03412         {
03413           GridTree_2_Inverted[j][k]->GetEntry(i);
03414           nc[j][k] = grid_nc*GridScale_Inverted_RHC;
03415           numucc[j][k] = grid_numucc*GridScale_Inverted_RHC;
03416           bnuecc[j][k] = grid_bnuecc*GridScale_Inverted_RHC;
03417           nutaucc[j][k] = grid_nutaucc*GridScale_Inverted_RHC;
03418           sig[j][k] = grid_nue*GridScale_Inverted_RHC;
03419         }
03420       }
03421     }
03422     nexp->Add(nexp_bkgd,nexp_signal,1,1);
03423     ErrCalc->SetGridPred(nBins,nc,numucc,bnuecc,nutaucc,sig);
03424     
03425     chi2min_IH[i]=GetMinLikelihood_Delta(false);
03426     
03427     ErrCalc->SetUseGrid(true);//use grid predictions set up above
03428     chi2data_IH[i] = 1e10;
03429     if(FitMethod==0)
03430     {
03431       chi2data_IH[i] = PoissonChi2(nexp);
03432     }
03433     else if(FitMethod==1)
03434     {
03435       chi2data_IH[i] = ScaledChi2(Bkgd,Sig);
03436     }
03437     else if(FitMethod==2)
03438     {
03439       chi2data_IH[i] = StandardChi2(nexp);
03440     }
03441     else if(FitMethod==3)
03442     {
03443       //Likelihood: "Standard" (N syst, N nuisance)
03444       //Calculate the likelihood (x2 for chi)
03445       Bkgd->Reset();
03446       Bkgd->Add(nexp_bkgd);
03447       Sig->Reset();
03448       Sig->Add(nexp_signal);
03449       chi2data_IH[i] = StandardLikelihood();
03450     }
03451     else if(FitMethod==4)
03452     {
03453       //Likelihood: Bin by Bin Calculation of Systematics
03454       //Calculate the likelihood (x2 for chi)
03455       Bkgd->Reset();
03456       Bkgd->Add(nexp_bkgd);
03457       Sig->Reset();
03458       Sig->Add(nexp_signal);
03459       chi2data_IH[i] = BinLikelihood();
03460     }
03461     else
03462     {
03463       cout<<"Error in RunMultiBinFC_MHDeltaFit(): Unknown 'FitMethod'."<<endl;
03464     }
03465   }
03466   
03467   if(nPts_Normal!=nPts_Inverted) cout<<"different number of points"<<endl;
03468   
03469   double *deltafit_nh = new double[nPts_Normal];
03470   double *deltafit_ih = new double[nPts_Inverted];
03471   double *mhfit_nh = new double[nPts_Normal];
03472   double *mhfit_ih = new double[nPts_Inverted];
03473   double delchi2,delchi2_norm,delchi2_invt,delchi2_bestmh;
03474   
03475   //delta fit
03476   for(i=0;i<nPts_Normal;i++)
03477   {
03478     //normal
03479     delchi2 = chi2data_NH[i] - chi2min_NH[i];
03480     chi2databin = int((delchi2-chi2start)/chi2binwidth)+1;
03481     
03482     TH1D *chi2hist = (TH1D*)f->Get(Form("Chi2Hist_DeltaFit_Normal_%i",i));
03483     
03484     if(chi2hist->Integral()<1)
03485     {
03486       cout<<"Warning, chi2hist is empty."<<endl;
03487       frac=0;
03488     }
03489     else
03490     {
03491       frac = chi2hist->Integral(1,chi2databin-1)/chi2hist->Integral();
03492     }
03493     if(chi2databin==1) frac=0;
03494     
03495     deltafit_nh[i] = frac;
03496     
03497     //inverted
03498     delchi2 = chi2data_IH[i] - chi2min_IH[i];
03499     chi2databin = int((delchi2-chi2start)/chi2binwidth)+1;
03500     
03501     chi2hist = (TH1D*)f->Get(Form("Chi2Hist_DeltaFit_Inverted_%i",i));
03502     
03503     if(chi2hist->Integral()<1)
03504     {
03505       cout<<"Warning, chi2hist is empty."<<endl;
03506       frac=0;
03507     }
03508     else
03509     {
03510       frac = chi2hist->Integral(1,chi2databin-1)/chi2hist->Integral();
03511     }
03512     if(chi2databin==1) frac=0;
03513     
03514     deltafit_ih[i] = frac;
03515     
03516     delete chi2hist;
03517   }
03518   
03519   //mh fit
03520   for(i=0;i<nPts_Normal;i++)
03521   {
03522     delchi2_norm = chi2data_NH[i] - chi2min_NH[i];
03523     delchi2_invt = chi2data_IH[i] - chi2min_IH[i];
03524     delchi2_bestmh = delchi2_invt;
03525     if(delchi2_norm<delchi2_bestmh) delchi2_bestmh = delchi2_norm;
03526     cout<<"NH: chi2min = "<<chi2min_NH[i]<<", IH: chi2min = "<<chi2min_IH[i]<<endl;
03527     
03528     delchi2 = delchi2_norm - delchi2_bestmh;
03529     chi2databin = int((delchi2-chi2start)/chi2binwidth)+1;
03530     
03531     TH1D *chi2hist = (TH1D*)f->Get(Form("Chi2Hist_MHFit_Normal_%i",i));
03532     
03533     if(chi2hist->Integral()<1)
03534     {
03535       cout<<"Warning, chi2hist is empty."<<endl;
03536       frac=0;
03537     }
03538     else
03539     {
03540       frac = chi2hist->Integral(1,chi2databin-1)/chi2hist->Integral();
03541     }
03542     if(chi2databin==1) frac=0;
03543     
03544     mhfit_nh[i] = frac;
03545     
03546     delchi2 = delchi2_invt - delchi2_bestmh;
03547     chi2databin = int((delchi2-chi2start)/chi2binwidth)+1;
03548     
03549     chi2hist = (TH1D*)f->Get(Form("Chi2Hist_MHFit_Inverted_%i",i));
03550     
03551     if(chi2hist->Integral()<1)
03552     {
03553       cout<<"Warning, chi2hist is empty."<<endl;
03554       frac=0;
03555     }
03556     else
03557     {
03558       frac = chi2hist->Integral(1,chi2databin-1)/chi2hist->Integral();
03559     }
03560     if(chi2databin==1) frac=0;
03561     
03562     mhfit_ih[i] = frac;
03563     
03564     delete chi2hist;
03565   }
03566 
03567   f->Close();
03568   
03569   TGraph *g_delta_NH = new TGraph(nPts_Normal,deltapts,deltafit_nh);
03570   g_delta_NH->SetName("delta_NH");
03571   
03572   TGraph *g_delta_IH = new TGraph(nPts_Normal,deltapts,deltafit_ih);
03573   g_delta_IH->SetName("delta_IH");
03574   
03575   TGraph *g_MH_NH = new TGraph(nPts_Normal,deltapts,mhfit_nh);
03576   g_MH_NH->SetName("mh_NH");
03577   
03578   TGraph *g_MH_IH = new TGraph(nPts_Normal,deltapts,mhfit_ih);
03579   g_MH_IH->SetName("mh_IH");
03580   
03581   TFile *fout = new TFile(gSystem->ExpandPathName(outFileName.c_str()),"RECREATE");
03582   g_delta_NH->Write();
03583   g_delta_IH->Write();
03584   g_MH_NH->Write();
03585   g_MH_IH->Write();
03586   fout->Close();
03587   
03588   delete [] chi2data_NH;
03589   delete [] chi2data_IH;
03590   delete [] chi2min_NH;
03591   delete [] chi2min_IH;
03592   delete [] deltafit_nh;
03593   delete [] deltafit_ih;
03594   delete [] mhfit_nh;
03595   delete [] mhfit_ih;
03596   
03597   return;
03598 }

void NueFit2D_Joint::RunMultiBinPseudoExpts ( bool  Print = true  )  [virtual]

Reimplemented from NueFit2D.

Definition at line 772 of file NueFit2D_Joint.cxx.

References NueFit2D::BinLikelihood(), NueFit2D::Bkgd, ErrorCalc::CalculateSystErrorMatrix(), ErrorCalc::CovMatrix, DefineBinDlnLMinuit(), DefineStdDlnLMinuit(), NueFit2D::ErrCalc, NueFit2D::ErrorMatrix, ExtrapFHC, ExtrapRHC, NueFit2D::FitMethod, Form(), NueFit2D::FracErr_Bkgd, NueFit2D::FracErr_Sig, NueFit2D::GenerateOneCorrelatedExp(), GetMinLikelihood(), NueFit2D::grid_background, NueFit2D::grid_bin_oscparerr, NueFit2D::grid_bnuecc, NueFit2D::grid_delta, NueFit2D::grid_nc, NueFit2D::grid_nue, NueFit2D::grid_numucc, NueFit2D::grid_nutaucc, NueFit2D::grid_signal, NueFit2D::grid_sinsq2th13, GridScale_Inverted_FHC, GridScale_Inverted_RHC, GridScale_Normal_FHC, GridScale_Normal_RHC, NueFit2D::GridTree_2_Inverted, NueFit2D::GridTree_2_Normal, NueFit2D::GridTree_Inverted, NueFit2D::GridTree_Normal, gSystem(), NueFit2D::IncludeOscParErrs, NueFit2D::InvErrorMatrix, NueFit2D::nBins, nBinsFHC, nBinsRHC, NueFit2D::NObs, NueFit2D::nPts_Inverted, NueFit2D::nPts_Normal, NueFit2D::NumExpts, NueFit2D::outFileName, NueFit2D::PoissonChi2(), ReadGridFiles(), NueFit2D::ScaledChi2(), ErrorCalc::SetGridPred(), ErrorCalc::SetUseGrid(), NueFit2D::Sig, NueFit2D::StandardChi2(), and NueFit2D::StandardLikelihood().

00773 {
00774   if(FitMethod==1 && (FracErr_Bkgd==0 || FracErr_Sig==0))
00775   {
00776     cout<<"FracErr_Bkgd and FracErr_Sig need to be set for ScaledChi2.  Quitting..."<<endl;
00777     return;
00778   }
00779   if((ExtrapFHC.size()+ExtrapRHC.size())==0)
00780   {
00781     cout<<"No Extrapolate2D input.  Quitting..."<<endl;
00782     return;
00783   }
00784   for(unsigned int ie=0;ie<ExtrapFHC.size();ie++)
00785   {
00786     ExtrapFHC[ie]->GetPrediction();
00787   }
00788   for(unsigned int ie=0;ie<ExtrapRHC.size();ie++)
00789   {
00790     ExtrapRHC[ie]->GetPrediction();
00791   }
00792   if(ErrCalc==0)
00793   {
00794     cout<<"Need to set ErrorCalc object!  Quitting..."<<endl;
00795     return;
00796   }
00797   if(FitMethod==3) DefineStdDlnLMinuit();
00798   if(FitMethod==4) DefineBinDlnLMinuit();
00799 
00800   nBinsFHC=0;
00801   nBinsRHC=0;  
00802   if (ExtrapFHC.size()>0) nBinsFHC = ExtrapFHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
00803   if (ExtrapRHC.size()>0) nBinsRHC = ExtrapRHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
00804   nBins = nBinsFHC + nBinsRHC;
00805 
00806   if(ErrorMatrix==0) ErrorMatrix = new TH2D("ErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
00807   if(InvErrorMatrix==0) InvErrorMatrix = new TH2D("InvErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
00808   
00809   TH2D *Error4Expts = new TH2D("Error4Expts","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
00810   
00811   ReadGridFiles();
00812   
00813   if(nPts_Normal==0 || nPts_Inverted==0) return;
00814   
00815   gRandom->SetSeed(0);
00816   
00817   int i,u;
00818   unsigned int j,k;
00819   TH1D *nexp_bkgd = new TH1D("nexp_bkgd","",nBins,-0.5,nBins-0.5);
00820   TH1D *nexp_signal = new TH1D("nexp_signal","",nBins,-0.5,nBins-0.5);
00821   TH1D *nexp = new TH1D("nexp","",nBins,-0.5,nBins-0.5);
00822   double delchi2,chi2min;
00823   TH1D *chi2hist = new TH1D("chi2hist","",110000,-10,100);
00824   double ele;
00825   int noff;
00826   
00827   vector< vector<double> > nc,numucc,bnuecc,nutaucc,sig;
00828   for(j=0;j<nBins;j++)
00829   {
00830     nc.push_back( vector<double>() );
00831     numucc.push_back( vector<double>() );
00832     bnuecc.push_back( vector<double>() );
00833     nutaucc.push_back( vector<double>() );
00834     sig.push_back( vector<double>() );
00835     if (j<nBinsFHC){
00836       for(k=0;k<ExtrapFHC.size();k++)
00837         {
00838           nc[j].push_back(0);
00839           numucc[j].push_back(0);
00840           bnuecc[j].push_back(0);
00841           nutaucc[j].push_back(0);
00842           sig[j].push_back(0);
00843         }
00844     } else if (j>=nBinsFHC){
00845       for(k=0;k<ExtrapRHC.size();k++)
00846         {
00847           nc[j].push_back(0);
00848           numucc[j].push_back(0);
00849           bnuecc[j].push_back(0);
00850           nutaucc[j].push_back(0);
00851           sig[j].push_back(0);
00852         }
00853     }
00854   }
00855   
00856   Bkgd = (TH1D*)nexp->Clone("Bkgd");
00857   Bkgd->Reset();
00858   Sig = (TH1D*)nexp->Clone("Sig");
00859   Sig->Reset();
00860   
00861   ofstream myfile;
00862   string file,ofile;
00863   
00864   //normal hierarchy
00865   
00866   TFile *f = new TFile(gSystem->ExpandPathName(outFileName.c_str()),"RECREATE");
00867   
00868   if(Print)
00869   {
00870     ofile = gSystem->ExpandPathName(outFileName.c_str());
00871     file = ofile.substr(0,ofile.length()-5) + "_Normal.dat";
00872     myfile.open(gSystem->ExpandPathName(file.c_str()));
00873   }
00874   
00875   for(i=0;i<nPts_Normal;i++)
00876   {
00877     cout<<"point "<<(i+1)<<"/"<<nPts_Normal<<" (normal hierarchy)"<<endl;
00878     
00879     nexp_bkgd->Reset();
00880     nexp_signal->Reset();
00881     nexp->Reset();
00882     
00883     for(j=0;j<nBins;j++)
00884     {
00885       GridTree_Normal[j]->GetEntry(i);
00886       if (j<nBinsFHC){
00887         nexp_bkgd->SetBinContent(j+1,grid_background*GridScale_Normal_FHC);
00888         nexp_signal->SetBinContent(j+1,grid_signal*GridScale_Normal_FHC);
00889         for(k=0;k<ExtrapFHC.size();k++)
00890           {
00891             GridTree_2_Normal[j][k]->GetEntry(i);
00892             nc[j][k] = grid_nc*GridScale_Normal_FHC;
00893             numucc[j][k] = grid_numucc*GridScale_Normal_FHC;
00894             bnuecc[j][k] = grid_bnuecc*GridScale_Normal_FHC;
00895             nutaucc[j][k] = grid_nutaucc*GridScale_Normal_FHC;
00896             sig[j][k] = grid_nue*GridScale_Normal_FHC;
00897           }
00898       } else if (j>=nBinsFHC){
00899         nexp_bkgd->SetBinContent(j+1,grid_background*GridScale_Normal_RHC);
00900         nexp_signal->SetBinContent(j+1,grid_signal*GridScale_Normal_RHC);
00901         for(k=0;k<ExtrapRHC.size();k++)
00902           {
00903             GridTree_2_Normal[j][k]->GetEntry(i);
00904             nc[j][k] = grid_nc*GridScale_Normal_RHC;
00905             numucc[j][k] = grid_numucc*GridScale_Normal_RHC;
00906             bnuecc[j][k] = grid_bnuecc*GridScale_Normal_RHC;
00907             nutaucc[j][k] = grid_nutaucc*GridScale_Normal_RHC;
00908             sig[j][k] = grid_nue*GridScale_Normal_RHC;
00909           }
00910       }
00911     }
00912 
00913 
00914     nexp->Add(nexp_bkgd,nexp_signal,1,1);
00915     ErrCalc->SetGridPred(nBinsFHC+nBinsRHC,nc,numucc,bnuecc,nutaucc,sig);
00916     
00917     Error4Expts->Reset();
00918     ErrCalc->SetUseGrid(true);
00919     ErrCalc->CalculateSystErrorMatrix();
00920     Error4Expts->Add(ErrCalc->CovMatrix);
00921     //ErrCalc->CalculateHOOError();
00922     //Error4Expts->Add(ErrCalc->CovMatrix_Decomp);
00923 
00924     if(IncludeOscParErrs)
00925     {
00926       noff=0;
00927       for(j=0;j<nBins;j++)
00928       {
00929         ele=Error4Expts->GetBinContent(j+1,j+1);
00930         ele+=(grid_bin_oscparerr[j]*grid_bin_oscparerr[j]*nexp->GetBinContent(j+1)*nexp->GetBinContent(j+1));
00931         Error4Expts->SetBinContent(j+1,j+1,ele);
00932         
00933         for(k=0;k<nBins;k++)
00934         {
00935           if(k>j)
00936           {
00937             ele=Error4Expts->GetBinContent(j+1,k+1);
00938             ele+=(grid_bin_oscparerr[j]*grid_bin_oscparerr[k]*nexp->GetBinContent(j+1)*nexp->GetBinContent(k+1));
00939             Error4Expts->SetBinContent(j+1,k+1,ele);
00940             
00941             ele=Error4Expts->GetBinContent(k+1,j+1);
00942             ele+=(grid_bin_oscparerr[j]*grid_bin_oscparerr[k]*nexp->GetBinContent(j+1)*nexp->GetBinContent(k+1));
00943             Error4Expts->SetBinContent(k+1,j+1,ele);
00944             
00945             noff++;
00946           }
00947         }
00948       }
00949     }
00950     
00951     chi2hist->Reset();
00952     chi2hist->SetName(Form("Chi2Hist_Normal_%i",i));
00953     
00954     for(u=0;u<NumExpts;u++)
00955     {
00956       cout<<"expt "<<(u+1)<<"/"<<NumExpts<<endl;
00957       
00958       GenerateOneCorrelatedExp(nexp,Error4Expts);
00959       if(Print)
00960       {
00961         myfile << grid_sinsq2th13 << " " << grid_delta << " ";
00962         for(j=0;j<nBins;j++)
00963         {
00964           myfile << NObs->GetBinContent(j+1) << " ";
00965         }
00966         myfile << endl;
00967       }
00968       
00969       chi2min=GetMinLikelihood(grid_delta,true);
00970       
00971       ErrCalc->SetUseGrid(true);//will use the grid predictions set above
00972       delchi2 = 1e10;
00973       if(FitMethod==0)
00974       {
00975         delchi2 = PoissonChi2(nexp) - chi2min;
00976       }
00977       else if(FitMethod==1)
00978       {
00979         delchi2 = ScaledChi2(nexp_bkgd,nexp_signal) - chi2min;
00980       }
00981       else if(FitMethod==2)
00982       {
00983         delchi2 = StandardChi2(nexp) - chi2min;
00984       }
00985       else if(FitMethod==3)
00986       {
00987         //Likelihood: "Standard" (N syst, N nuisance)
00988         //Calculate the likelihood (x2 for chi)
00989         Bkgd->Reset();
00990         Bkgd->Add(nexp_bkgd);
00991         Sig->Reset();
00992         Sig->Add(nexp_signal);
00993         delchi2 = StandardLikelihood() - chi2min;
00994       }
00995       else if(FitMethod==4)
00996       {
00997         //Likelihood: Bin by Bin Calculation of Systematics
00998         //Calculate the likelihood (x2 for chi)
00999         Bkgd->Reset();
01000         Bkgd->Add(nexp_bkgd);
01001         Sig->Reset();
01002         Sig->Add(nexp_signal);
01003 
01004         delchi2 = BinLikelihood() - chi2min;
01005       }
01006       else
01007       {
01008         cout<<"Error in RunMultiBinPseudoExpts(): Unknown 'FitMethod'."<<endl;
01009       }
01010       chi2hist->Fill(delchi2);
01011       
01012     }
01013     f->cd();
01014     chi2hist->Write();
01015     f->Close();
01016     
01017     f = new TFile(gSystem->ExpandPathName(outFileName.c_str()),"UPDATE");
01018   }
01019   
01020   if(Print) myfile.close();
01021   
01022   for(j=0;j<nBins;j++)
01023   {
01024     GridTree_Normal[j]->Write();
01025     
01026     if (j<nBinsFHC){
01027       for(k=0;k<ExtrapFHC.size();k++)
01028         {
01029           GridTree_2_Normal[j][k]->Write();
01030         }
01031     } else if (j>=nBinsFHC){
01032       for(k=0;k<ExtrapRHC.size();k++)
01033         {
01034           GridTree_2_Normal[j][k]->Write();
01035         }
01036     }
01037   }
01038   
01039   f->Close();
01040   
01041   //inverted hierarchy
01042   
01043   f = new TFile(gSystem->ExpandPathName(outFileName.c_str()),"UPDATE");
01044   
01045   if(Print)
01046   {
01047     ofile = gSystem->ExpandPathName(outFileName.c_str());
01048     file = ofile.substr(0,ofile.length()-5) + "_Inverted.dat";
01049     myfile.open(gSystem->ExpandPathName(file.c_str()));
01050   }
01051   
01052   for(i=0;i<nPts_Inverted;i++)
01053   {
01054     cout<<"point "<<(i+1)<<"/"<<nPts_Inverted<<" (inverted hierarchy)"<<endl;
01055     
01056     nexp_bkgd->Reset();
01057     nexp_signal->Reset();
01058     nexp->Reset();
01059     
01060     for(j=0;j<nBins;j++)
01061     {
01062       GridTree_Inverted[j]->GetEntry(i);
01063       if (j<nBinsFHC){
01064         nexp_bkgd->SetBinContent(j+1,grid_background*GridScale_Inverted_FHC);
01065         nexp_signal->SetBinContent(j+1,grid_signal*GridScale_Inverted_FHC);
01066         for(k=0;k<ExtrapFHC.size();k++)
01067           {
01068             GridTree_2_Inverted[j][k]->GetEntry(i);
01069             nc[j][k] = grid_nc*GridScale_Inverted_FHC;
01070             numucc[j][k] = grid_numucc*GridScale_Inverted_FHC;
01071             bnuecc[j][k] = grid_bnuecc*GridScale_Inverted_FHC;
01072             nutaucc[j][k] = grid_nutaucc*GridScale_Inverted_FHC;
01073             sig[j][k] = grid_nue*GridScale_Inverted_FHC;
01074           }
01075       } else if (j>=nBinsFHC){
01076         nexp_bkgd->SetBinContent(j+1,grid_background*GridScale_Inverted_RHC);
01077         nexp_signal->SetBinContent(j+1,grid_signal*GridScale_Inverted_RHC);
01078         for(k=0;k<ExtrapRHC.size();k++)
01079           {
01080             GridTree_2_Inverted[j][k]->GetEntry(i);
01081             nc[j][k] = grid_nc*GridScale_Inverted_RHC;
01082             numucc[j][k] = grid_numucc*GridScale_Inverted_RHC;
01083             bnuecc[j][k] = grid_bnuecc*GridScale_Inverted_RHC;
01084             nutaucc[j][k] = grid_nutaucc*GridScale_Inverted_RHC;
01085             sig[j][k] = grid_nue*GridScale_Inverted_RHC;
01086           }
01087       }
01088     }
01089 
01090     nexp->Add(nexp_bkgd,nexp_signal,1,1);
01091     ErrCalc->SetGridPred(nBinsFHC+nBinsRHC,nc,numucc,bnuecc,nutaucc,sig);
01092     
01093     Error4Expts->Reset();
01094     ErrCalc->SetUseGrid(true);
01095     ErrCalc->CalculateSystErrorMatrix();
01096     Error4Expts->Add(ErrCalc->CovMatrix);
01097     //ErrCalc->CalculateHOOError();
01098     //Error4Expts->Add(ErrCalc->CovMatrix_Decomp);
01099     if(IncludeOscParErrs)
01100     {
01101       noff=0;
01102       for(j=0;j<nBins;j++)
01103       {
01104         ele=Error4Expts->GetBinContent(j+1,j+1);
01105         ele+=(grid_bin_oscparerr[j]*grid_bin_oscparerr[j]*nexp->GetBinContent(j+1)*nexp->GetBinContent(j+1));
01106         Error4Expts->SetBinContent(j+1,j+1,ele);
01107         
01108         for(k=0;k<nBins;k++)
01109         {
01110           if(k>j)
01111           {
01112             ele=Error4Expts->GetBinContent(j+1,k+1);
01113             ele+=(grid_bin_oscparerr[j]*grid_bin_oscparerr[k]*nexp->GetBinContent(j+1)*nexp->GetBinContent(k+1));
01114             Error4Expts->SetBinContent(j+1,k+1,ele);
01115             
01116             ele=Error4Expts->GetBinContent(k+1,j+1);
01117             ele+=(grid_bin_oscparerr[j]*grid_bin_oscparerr[k]*nexp->GetBinContent(j+1)*nexp->GetBinContent(k+1));
01118             Error4Expts->SetBinContent(k+1,j+1,ele);
01119             
01120             noff++;
01121           }
01122         }
01123       }
01124     }
01125     
01126     chi2hist->Reset();
01127     chi2hist->SetName(Form("Chi2Hist_Inverted_%i",i));
01128     
01129     for(u=0;u<NumExpts;u++)
01130     {
01131       cout<<"expt "<<(u+1)<<"/"<<NumExpts<<endl;
01132       
01133       GenerateOneCorrelatedExp(nexp,Error4Expts);
01134       if(Print)
01135       {
01136         myfile << grid_sinsq2th13 << " " << grid_delta << " ";
01137         for(j=0;j<nBins;j++)
01138         {
01139           myfile << NObs->GetBinContent(j+1) << " ";
01140         }
01141         myfile << endl;
01142       }
01143       
01144       chi2min=GetMinLikelihood(grid_delta,false);
01145       
01146       ErrCalc->SetUseGrid(true);//will use the grid predictions set above
01147       delchi2 = 1e10;
01148       if(FitMethod==0)
01149       {
01150         delchi2 = PoissonChi2(nexp) - chi2min;
01151       }
01152       else if(FitMethod==1)
01153       {
01154         delchi2 = ScaledChi2(nexp_bkgd,nexp_signal) - chi2min;
01155       }
01156       else if(FitMethod==2)
01157       {
01158         delchi2 = StandardChi2(nexp) - chi2min;
01159       }
01160       else if(FitMethod==3)
01161       {
01162         //Likelihood: "Standard" (N syst, N nuisance)
01163         //Calculate the likelihood (x2 for chi)
01164         Bkgd->Reset();
01165         Bkgd->Add(nexp_bkgd);
01166         Sig->Reset();
01167         Sig->Add(nexp_signal);
01168         delchi2 = StandardLikelihood() - chi2min;
01169       }
01170       else if(FitMethod==4)
01171       {
01172         //Likelihood: Bin by Bin Calculation of Systematics
01173         //Calculate the likelihood (x2 for chi)
01174         Bkgd->Reset();
01175         Bkgd->Add(nexp_bkgd);
01176         Sig->Reset();
01177         Sig->Add(nexp_signal);
01178         delchi2 = BinLikelihood() - chi2min;
01179       }
01180       else
01181       {
01182         cout<<"Error in RunMultiBinPseudoExpts(): Unknown 'FitMethod'."<<endl;
01183       }
01184       chi2hist->Fill(delchi2);
01185     }
01186     f->cd();
01187     chi2hist->Write();
01188     f->Close();
01189     f = new TFile(gSystem->ExpandPathName(outFileName.c_str()),"UPDATE");
01190   }
01191   
01192   if(Print) myfile.close();
01193   
01194   for(j=0;j<nBins;j++)
01195   {
01196     GridTree_Inverted[j]->Write();
01197     if (j<nBinsFHC){
01198       for(k=0;k<ExtrapFHC.size();k++)
01199         {
01200           GridTree_2_Inverted[j][k]->Write();
01201         }
01202     } else if (j>=nBinsFHC){
01203       for(k=0;k<ExtrapRHC.size();k++)
01204         {
01205           GridTree_2_Inverted[j][k]->Write();
01206         }
01207     }
01208   }
01209   
01210   f->Close();
01211   
01212   return;
01213 }

void NueFit2D_Joint::RunMultiBinPseudoExpts_MHDeltaFit ( bool  Print = true  )  [virtual]

Reimplemented from NueFit2D.

Definition at line 2535 of file NueFit2D_Joint.cxx.

References NueFit2D::BinLikelihood(), NueFit2D::Bkgd, ErrorCalc::CalculateSystErrorMatrix(), ErrorCalc::CovMatrix, DefineBinDlnLMinuit(), DefineStdDlnLMinuit(), NueFit2D::ErrCalc, NueFit2D::ErrorMatrix, ExtrapFHC, ExtrapRHC, NueFit2D::FitMethod, Form(), NueFit2D::FracErr_Bkgd, NueFit2D::FracErr_Sig, NueFit2D::GenerateOneCorrelatedExp(), GetMinLikelihood_Delta(), NueFit2D::grid_background, NueFit2D::grid_bin_oscparerr, NueFit2D::grid_bnuecc, NueFit2D::grid_delta, NueFit2D::grid_nc, NueFit2D::grid_nue, NueFit2D::grid_numucc, NueFit2D::grid_nutaucc, NueFit2D::grid_signal, NueFit2D::grid_sinsq2th13, GridScale_Inverted_FHC, GridScale_Inverted_RHC, GridScale_Normal_FHC, GridScale_Normal_RHC, NueFit2D::GridTree_2_Inverted, NueFit2D::GridTree_2_Normal, NueFit2D::GridTree_Inverted, NueFit2D::GridTree_Normal, gSystem(), NueFit2D::IncludeOscParErrs, NueFit2D::InvErrorMatrix, NueFit2D::nBins, nBinsFHC, nBinsRHC, NueFit2D::NObs, NueFit2D::nPts_Inverted, NueFit2D::nPts_Normal, NueFit2D::NumExpts, NueFit2D::outFileName, NueFit2D::PoissonChi2(), ReadGridFiles(), NueFit2D::ScaledChi2(), ErrorCalc::SetGridPred(), ErrorCalc::SetUseGrid(), NueFit2D::Sig, NueFit2D::StandardChi2(), and NueFit2D::StandardLikelihood().

02536 {
02537   if(FitMethod==1 && (FracErr_Bkgd==0 || FracErr_Sig==0))
02538   {
02539     cout<<"FracErr_Bkgd and FracErr_Sig need to be set for ScaledChi2.  Quitting..."<<endl;
02540     return;
02541   }
02542   if((ExtrapFHC.size()+ExtrapRHC.size())==0)
02543   {
02544     cout<<"No Extrapolate2D input.  Quitting..."<<endl;
02545     return;
02546   }
02547   for(unsigned int ie=0;ie<ExtrapFHC.size();ie++)
02548   {
02549     ExtrapFHC[ie]->GetPrediction();
02550   }
02551   for(unsigned int ie=0;ie<ExtrapRHC.size();ie++)
02552   {
02553     ExtrapRHC[ie]->GetPrediction();
02554   }
02555   if(ErrCalc==0)
02556   {
02557     cout<<"Need to set ErrorCalc object!  Quitting..."<<endl;
02558     return;
02559   }
02560   if(FitMethod==3) DefineStdDlnLMinuit();
02561   if(FitMethod==4) DefineBinDlnLMinuit();
02562   
02563   nBinsFHC=0;
02564   nBinsRHC=0;  
02565   if (ExtrapFHC.size()>0) nBinsFHC = ExtrapFHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
02566   if (ExtrapRHC.size()>0) nBinsRHC = ExtrapRHC[0]->Pred_TotalBkgd_VsBinNumber->GetNbinsX();
02567   nBins = nBinsFHC + nBinsRHC;
02568   
02569   if(ErrorMatrix==0) ErrorMatrix = new TH2D("ErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
02570   if(InvErrorMatrix==0) InvErrorMatrix = new TH2D("InvErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
02571   
02572   TH2D *Error4Expts = new TH2D("Error4Expts","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
02573   
02574   ReadGridFiles();
02575   
02576   if(nPts_Normal==0 || nPts_Inverted==0) return;
02577   
02578   gRandom->SetSeed(0);
02579   
02580   int i,u;
02581   unsigned int j,k;
02582   TH1D *nexp_bkgd_n = new TH1D("nexp_bkgd_n","",nBins,-0.5,nBins-0.5);
02583   TH1D *nexp_signal_n = new TH1D("nexp_signal_n","",nBins,-0.5,nBins-0.5);
02584   TH1D *nexp_n = new TH1D("nexp_n","",nBins,-0.5,nBins-0.5);
02585   
02586   TH1D *nexp_bkgd_i = new TH1D("nexp_bkgd_i","",nBins,-0.5,nBins-0.5);
02587   TH1D *nexp_signal_i = new TH1D("nexp_signal_i","",nBins,-0.5,nBins-0.5);
02588   TH1D *nexp_i = new TH1D("nexp_i","",nBins,-0.5,nBins-0.5);
02589   
02590   double chi2_norm,chi2_invt,chi2min_norm,chi2min_invt;
02591   double delchi2_norm, delchi2_invt,delchi2_bestmh;
02592   
02593   TH1D *chi2hist_deltafit = new TH1D("chi2hist_deltafit","",110000,-10,100);
02594   TH1D *chi2hist_mhfit = new TH1D("chi2hist_mhfit","",110000,-10,100);
02595   double ele;
02596   int noff;
02597   
02598   vector< vector<double> > nc_n,numucc_n,bnuecc_n,nutaucc_n,sig_n;
02599   vector< vector<double> > nc_i,numucc_i,bnuecc_i,nutaucc_i,sig_i;
02600   for(j=0;j<nBins;j++)
02601   {
02602     nc_n.push_back( vector<double>() );
02603     numucc_n.push_back( vector<double>() );
02604     bnuecc_n.push_back( vector<double>() );
02605     nutaucc_n.push_back( vector<double>() );
02606     sig_n.push_back( vector<double>() );
02607     
02608     nc_i.push_back( vector<double>() );
02609     numucc_i.push_back( vector<double>() );
02610     bnuecc_i.push_back( vector<double>() );
02611     nutaucc_i.push_back( vector<double>() );
02612     sig_i.push_back( vector<double>() );
02613     
02614     if (j<nBinsFHC){
02615       for(k=0;k<ExtrapFHC.size();k++)
02616       {
02617         nc_n[j].push_back(0);
02618         numucc_n[j].push_back(0);
02619         bnuecc_n[j].push_back(0);
02620         nutaucc_n[j].push_back(0);
02621         sig_n[j].push_back(0);
02622         
02623         nc_i[j].push_back(0);
02624         numucc_i[j].push_back(0);
02625         bnuecc_i[j].push_back(0);
02626         nutaucc_i[j].push_back(0);
02627         sig_i[j].push_back(0);
02628       }
02629     } else if (j>=nBinsFHC){
02630       for(k=0;k<ExtrapRHC.size();k++)
02631       {
02632         nc_n[j].push_back(0);
02633         numucc_n[j].push_back(0);
02634         bnuecc_n[j].push_back(0);
02635         nutaucc_n[j].push_back(0);
02636         sig_n[j].push_back(0);
02637         
02638         nc_i[j].push_back(0);
02639         numucc_i[j].push_back(0);
02640         bnuecc_i[j].push_back(0);
02641         nutaucc_i[j].push_back(0);
02642         sig_i[j].push_back(0);
02643       }
02644     }
02645   }
02646   
02647   Bkgd = (TH1D*)nexp_n->Clone("Bkgd");
02648   Bkgd->Reset();
02649   Sig = (TH1D*)nexp_n->Clone("Sig");
02650   Sig->Reset();
02651   
02652   ofstream myfile;
02653   string file,ofile;
02654   
02655   //normal hierarchy
02656   
02657   TFile *f = new TFile(gSystem->ExpandPathName(outFileName.c_str()),"RECREATE");
02658   
02659   if(Print)
02660   {
02661     ofile = gSystem->ExpandPathName(outFileName.c_str());
02662     file = ofile.substr(0,ofile.length()-5) + "_Normal.dat";
02663     myfile.open(gSystem->ExpandPathName(file.c_str()));
02664   }
02665   
02666   for(i=0;i<nPts_Normal;i++)
02667   {
02668     cout<<"point "<<(i+1)<<"/"<<nPts_Normal<<" (normal hierarchy)"<<endl;
02669     
02670     nexp_bkgd_n->Reset();
02671     nexp_signal_n->Reset();
02672     nexp_n->Reset();
02673     
02674     nexp_bkgd_i->Reset();
02675     nexp_signal_i->Reset();
02676     nexp_i->Reset();
02677     
02678     for(j=0;j<nBins;j++)
02679     {
02680       GridTree_Normal[j]->GetEntry(i);
02681       if (j<nBinsFHC){
02682         nexp_bkgd_n->SetBinContent(j+1,grid_background*GridScale_Normal_FHC);
02683         nexp_signal_n->SetBinContent(j+1,grid_signal*GridScale_Normal_FHC);
02684         for(k=0;k<ExtrapFHC.size();k++)
02685         {
02686           GridTree_2_Normal[j][k]->GetEntry(i);
02687           nc_n[j][k] = grid_nc*GridScale_Normal_FHC;
02688           numucc_n[j][k] = grid_numucc*GridScale_Normal_FHC;
02689           bnuecc_n[j][k] = grid_bnuecc*GridScale_Normal_FHC;
02690           nutaucc_n[j][k] = grid_nutaucc*GridScale_Normal_FHC;
02691           sig_n[j][k] = grid_nue*GridScale_Normal_FHC;
02692         }
02693       } else if (j>=nBinsFHC){
02694         nexp_bkgd_n->SetBinContent(j+1,grid_background*GridScale_Normal_RHC);
02695         nexp_signal_n->SetBinContent(j+1,grid_signal*GridScale_Normal_RHC);
02696         for(k=0;k<ExtrapRHC.size();k++)
02697         {
02698           GridTree_2_Normal[j][k]->GetEntry(i);
02699           nc_n[j][k] = grid_nc*GridScale_Normal_RHC;
02700           numucc_n[j][k] = grid_numucc*GridScale_Normal_RHC;
02701           bnuecc_n[j][k] = grid_bnuecc*GridScale_Normal_RHC;
02702           nutaucc_n[j][k] = grid_nutaucc*GridScale_Normal_RHC;
02703           sig_n[j][k] = grid_nue*GridScale_Normal_RHC;
02704         }
02705       }
02706     }
02707     for(j=0;j<nBins;j++)
02708     {
02709       GridTree_Inverted[j]->GetEntry(i);
02710       if (j<nBinsFHC){
02711         nexp_bkgd_i->SetBinContent(j+1,grid_background*GridScale_Inverted_FHC);
02712         nexp_signal_i->SetBinContent(j+1,grid_signal*GridScale_Inverted_FHC);
02713         for(k=0;k<ExtrapFHC.size();k++)
02714         {
02715           GridTree_2_Inverted[j][k]->GetEntry(i);
02716           nc_i[j][k] = grid_nc*GridScale_Inverted_FHC;
02717           numucc_i[j][k] = grid_numucc*GridScale_Inverted_FHC;
02718           bnuecc_i[j][k] = grid_bnuecc*GridScale_Inverted_FHC;
02719           nutaucc_i[j][k] = grid_nutaucc*GridScale_Inverted_FHC;
02720           sig_i[j][k] = grid_nue*GridScale_Inverted_FHC;
02721         }
02722       } else if (j>=nBinsFHC){
02723         nexp_bkgd_i->SetBinContent(j+1,grid_background*GridScale_Inverted_RHC);
02724         nexp_signal_i->SetBinContent(j+1,grid_signal*GridScale_Inverted_RHC);
02725         for(k=0;k<ExtrapRHC.size();k++)
02726         {
02727           GridTree_2_Inverted[j][k]->GetEntry(i);
02728           nc_i[j][k] = grid_nc*GridScale_Inverted_RHC;
02729           numucc_i[j][k] = grid_numucc*GridScale_Inverted_RHC;
02730           bnuecc_i[j][k] = grid_bnuecc*GridScale_Inverted_RHC;
02731           nutaucc_i[j][k] = grid_nutaucc*GridScale_Inverted_RHC;
02732           sig_i[j][k] = grid_nue*GridScale_Inverted_RHC;
02733         }
02734       }
02735     }
02736     nexp_n->Add(nexp_bkgd_n,nexp_signal_n,1,1);
02737     nexp_i->Add(nexp_bkgd_i,nexp_signal_i,1,1);
02738     
02739     //want to generate pseudos for normal hierarchy
02740     ErrCalc->SetGridPred(nBins,nc_n,numucc_n,bnuecc_n,nutaucc_n,sig_n);
02741     
02742     Error4Expts->Reset();
02743     ErrCalc->SetUseGrid(true);
02744     ErrCalc->CalculateSystErrorMatrix();
02745     Error4Expts->Add(ErrCalc->CovMatrix);
02746     //ErrCalc->CalculateHOOError();
02747     //Error4Expts->Add(ErrCalc->CovMatrix_Decomp);
02748     
02749     if(IncludeOscParErrs)
02750     {
02751       noff=0;
02752       for(j=0;j<nBins;j++)
02753       {
02754         GridTree_Normal[j]->GetEntry(i);
02755         ele=Error4Expts->GetBinContent(j+1,j+1);
02756         ele+=(grid_bin_oscparerr[j]*grid_bin_oscparerr[j]*nexp_n->GetBinContent(j+1)*nexp_n->GetBinContent(j+1));
02757         Error4Expts->SetBinContent(j+1,j+1,ele);
02758         
02759         for(k=0;k<nBins;k++)
02760         {
02761           if(k>j)
02762           {
02763             ele=Error4Expts->GetBinContent(j+1,k+1);
02764             ele+=(grid_bin_oscparerr[j]*grid_bin_oscparerr[k]*nexp_n->GetBinContent(j+1)*nexp_n->GetBinContent(k+1));
02765             Error4Expts->SetBinContent(j+1,k+1,ele);
02766             
02767             ele=Error4Expts->GetBinContent(k+1,j+1);
02768             ele+=(grid_bin_oscparerr[j]*grid_bin_oscparerr[k]*nexp_n->GetBinContent(j+1)*nexp_n->GetBinContent(k+1));
02769             Error4Expts->SetBinContent(k+1,j+1,ele);
02770             
02771             noff++;
02772           }
02773         }
02774       }
02775     }
02776     
02777     chi2hist_deltafit->Reset();
02778     chi2hist_deltafit->SetName(Form("Chi2Hist_DeltaFit_Normal_%i",i));
02779     
02780     chi2hist_mhfit->Reset();
02781     chi2hist_mhfit->SetName(Form("Chi2Hist_MHFit_Normal_%i",i));
02782     
02783     for(u=0;u<NumExpts;u++)
02784     {
02785       cout<<"expt "<<(u+1)<<"/"<<NumExpts<<endl;
02786       
02787       GenerateOneCorrelatedExp(nexp_n,Error4Expts);
02788       if(Print)
02789       {
02790         myfile << grid_sinsq2th13 << " " << grid_delta << " ";
02791         for(j=0;j<nBins;j++)
02792         {
02793           myfile << NObs->GetBinContent(j+1) << " ";
02794         }
02795       }
02796       
02797       chi2min_norm=GetMinLikelihood_Delta(true);
02798       chi2min_invt=GetMinLikelihood_Delta(false);
02799       
02800       ErrCalc->SetGridPred(nBins,nc_n,numucc_n,bnuecc_n,nutaucc_n,sig_n);
02801       ErrCalc->SetUseGrid(true);
02802       chi2_norm = 1e10;
02803       if(FitMethod==0)
02804       {
02805         chi2_norm = PoissonChi2(nexp_n);
02806       }
02807       else if(FitMethod==1)
02808       {
02809         chi2_norm = ScaledChi2(nexp_bkgd_n,nexp_signal_n);
02810       }
02811       else if(FitMethod==2)
02812       {
02813         chi2_norm = StandardChi2(nexp_n);
02814       }
02815       else if(FitMethod==3)
02816       {
02817         //Likelihood: "Standard" (N syst, N nuisance)
02818         //Calculate the likelihood (x2 for chi)
02819         Bkgd->Reset();
02820         Bkgd->Add(nexp_bkgd_n);
02821         Sig->Reset();
02822         Sig->Add(nexp_signal_n);
02823         chi2_norm = StandardLikelihood();
02824       }
02825       else if(FitMethod==4)
02826       {
02827         //Likelihood: Bin by Bin Calculation of Systematics
02828         //Calculate the likelihood (x2 for chi)
02829         Bkgd->Reset();
02830         Bkgd->Add(nexp_bkgd_n);
02831         Sig->Reset();
02832         Sig->Add(nexp_signal_n);
02833         chi2_norm = BinLikelihood();
02834       }
02835       else
02836       {
02837         cout<<"Error in RunMultiBinPseudoExpts_MHDeltaFit(): Unknown 'FitMethod'."<<endl;
02838       }
02839       
02840       ErrCalc->SetGridPred(nBins,nc_i,numucc_i,bnuecc_i,nutaucc_i,sig_i);
02841       ErrCalc->SetUseGrid(true);
02842       chi2_invt = 1e10;
02843       if(FitMethod==0)
02844       {
02845         chi2_invt = PoissonChi2(nexp_i);
02846       }
02847       else if(FitMethod==1)
02848       {
02849         chi2_invt = ScaledChi2(nexp_bkgd_i,nexp_signal_i);
02850       }
02851       else if(FitMethod==2)
02852       {
02853         chi2_invt = StandardChi2(nexp_i);
02854       }
02855       else if(FitMethod==3)
02856       {
02857         //Likelihood: "Standard" (N syst, N nuisance)
02858         //Calculate the likelihood (x2 for chi)
02859         Bkgd->Reset();
02860         Bkgd->Add(nexp_bkgd_i);
02861         Sig->Reset();
02862         Sig->Add(nexp_signal_i);
02863         chi2_invt = StandardLikelihood();
02864       }
02865       else if(FitMethod==4)
02866       {
02867         //Likelihood: Bin by Bin Calculation of Systematics
02868         //Calculate the likelihood (x2 for chi)
02869         Bkgd->Reset();
02870         Bkgd->Add(nexp_bkgd_i);
02871         Sig->Reset();
02872         Sig->Add(nexp_signal_i);
02873         chi2_invt = BinLikelihood();
02874       }
02875       else
02876       {
02877         cout<<"Error in RunMultiBinPseudoExpts_MHDeltaFit(): Unknown 'FitMethod'."<<endl;
02878       }
02879       
02880       if(Print)
02881       {
02882         myfile << chi2_norm << " " << chi2_invt << " " << chi2min_norm << " " << chi2min_invt << endl;
02883       }
02884       
02885       delchi2_norm = chi2_norm - chi2min_norm;
02886       delchi2_invt = chi2_invt - chi2min_invt;
02887       delchi2_bestmh = delchi2_norm;
02888       if(delchi2_invt<delchi2_bestmh) delchi2_bestmh = delchi2_invt;
02889       
02890       chi2hist_deltafit->Fill(delchi2_norm);
02891       chi2hist_mhfit->Fill(delchi2_norm - delchi2_bestmh);
02892     }
02893     f->cd();
02894     chi2hist_deltafit->Write();
02895     chi2hist_mhfit->Write();
02896     f->Close();
02897     
02898     f = new TFile(gSystem->ExpandPathName(outFileName.c_str()),"UPDATE");
02899   }
02900   
02901   if(Print) myfile.close();
02902   
02903   for(j=0;j<nBins;j++)
02904   {
02905     GridTree_Normal[j]->Write();
02906     
02907     if (j<nBinsFHC){
02908       for(k=0;k<ExtrapFHC.size();k++)
02909       {
02910         GridTree_2_Normal[j][k]->Write();
02911       }
02912     } else if (j>=nBinsFHC){
02913       for(k=0;k<ExtrapRHC.size();k++)
02914       {
02915         GridTree_2_Normal[j][k]->Write();
02916       }
02917     }
02918   }
02919   
02920   f->Close();
02921   
02922   //inverted hierarchy
02923   
02924   f = new TFile(gSystem->ExpandPathName(outFileName.c_str()),"UPDATE");
02925   
02926   if(Print)
02927   {
02928     ofile = gSystem->ExpandPathName(outFileName.c_str());
02929     file = ofile.substr(0,ofile.length()-5) + "_Inverted.dat";
02930     myfile.open(gSystem->ExpandPathName(file.c_str()));
02931   }
02932   
02933   for(i=0;i<nPts_Inverted;i++)
02934   {
02935     cout<<"point "<<(i+1)<<"/"<<nPts_Inverted<<" (inverted hierarchy)"<<endl;
02936     
02937     nexp_bkgd_n->Reset();
02938     nexp_signal_n->Reset();
02939     nexp_n->Reset();
02940     
02941     nexp_bkgd_i->Reset();
02942     nexp_signal_i->Reset();
02943     nexp_i->Reset();
02944     
02945     for(j=0;j<nBins;j++)
02946     {
02947       GridTree_Normal[j]->GetEntry(i);
02948       if (j<nBinsFHC){
02949         nexp_bkgd_n->SetBinContent(j+1,grid_background*GridScale_Normal_FHC);
02950         nexp_signal_n->SetBinContent(j+1,grid_signal*GridScale_Normal_FHC);
02951         for(k=0;k<ExtrapFHC.size();k++)
02952         {
02953           GridTree_2_Normal[j][k]->GetEntry(i);
02954           nc_n[j][k] = grid_nc*GridScale_Normal_FHC;
02955           numucc_n[j][k] = grid_numucc*GridScale_Normal_FHC;
02956           bnuecc_n[j][k] = grid_bnuecc*GridScale_Normal_FHC;
02957           nutaucc_n[j][k] = grid_nutaucc*GridScale_Normal_FHC;
02958           sig_n[j][k] = grid_nue*GridScale_Normal_FHC;
02959         }
02960       } else if (j>=nBinsFHC){
02961         nexp_bkgd_n->SetBinContent(j+1,grid_background*GridScale_Normal_RHC);
02962         nexp_signal_n->SetBinContent(j+1,grid_signal*GridScale_Normal_RHC);
02963         for(k=0;k<ExtrapRHC.size();k++)
02964         {
02965           GridTree_2_Normal[j][k]->GetEntry(i);
02966           nc_n[j][k] = grid_nc*GridScale_Normal_RHC;
02967           numucc_n[j][k] = grid_numucc*GridScale_Normal_RHC;
02968           bnuecc_n[j][k] = grid_bnuecc*GridScale_Normal_RHC;
02969           nutaucc_n[j][k] = grid_nutaucc*GridScale_Normal_RHC;
02970           sig_n[j][k] = grid_nue*GridScale_Normal_RHC;
02971         }
02972       }
02973     }
02974     for(j=0;j<nBins;j++)
02975     {
02976       GridTree_Inverted[j]->GetEntry(i);
02977       if (j<nBinsFHC){
02978         nexp_bkgd_i->SetBinContent(j+1,grid_background*GridScale_Inverted_FHC);
02979         nexp_signal_i->SetBinContent(j+1,grid_signal*GridScale_Inverted_FHC);
02980         for(k=0;k<ExtrapFHC.size();k++)
02981         {
02982           GridTree_2_Inverted[j][k]->GetEntry(i);
02983           nc_i[j][k] = grid_nc*GridScale_Inverted_FHC;
02984           numucc_i[j][k] = grid_numucc*GridScale_Inverted_FHC;
02985           bnuecc_i[j][k] = grid_bnuecc*GridScale_Inverted_FHC;
02986           nutaucc_i[j][k] = grid_nutaucc*GridScale_Inverted_FHC;
02987           sig_i[j][k] = grid_nue*GridScale_Inverted_FHC;
02988         }
02989       } else if (j>=nBinsFHC){
02990         nexp_bkgd_i->SetBinContent(j+1,grid_background*GridScale_Inverted_RHC);
02991         nexp_signal_i->SetBinContent(j+1,grid_signal*GridScale_Inverted_RHC);
02992         for(k=0;k<ExtrapRHC.size();k++)
02993         {
02994           GridTree_2_Inverted[j][k]->GetEntry(i);
02995           nc_i[j][k] = grid_nc*GridScale_Inverted_RHC;
02996           numucc_i[j][k] = grid_numucc*GridScale_Inverted_RHC;
02997           bnuecc_i[j][k] = grid_bnuecc*GridScale_Inverted_RHC;
02998           nutaucc_i[j][k] = grid_nutaucc*GridScale_Inverted_RHC;
02999           sig_i[j][k] = grid_nue*GridScale_Inverted_RHC;
03000         }
03001       }
03002     }
03003     nexp_n->Add(nexp_bkgd_n,nexp_signal_n,1,1);
03004     nexp_i->Add(nexp_bkgd_i,nexp_signal_i,1,1);
03005     
03006     //want to generate pseudos for inverted hierarchy
03007     ErrCalc->SetGridPred(nBins,nc_i,numucc_i,bnuecc_i,nutaucc_i,sig_i);
03008     Error4Expts->Reset();
03009     ErrCalc->SetUseGrid(true);
03010     ErrCalc->CalculateSystErrorMatrix();
03011     Error4Expts->Add(ErrCalc->CovMatrix);
03012     //ErrCalc->CalculateHOOError();
03013     //Error4Expts->Add(ErrCalc->CovMatrix_Decomp);
03014     if(IncludeOscParErrs)
03015     {
03016       noff=0;
03017       for(j=0;j<nBins;j++)
03018       {
03019         ele=Error4Expts->GetBinContent(j+1,j+1);
03020         ele+=(grid_bin_oscparerr[j]*grid_bin_oscparerr[j]*nexp_i->GetBinContent(j+1)*nexp_i->GetBinContent(j+1));
03021         Error4Expts->SetBinContent(j+1,j+1,ele);
03022         
03023         for(k=0;k<nBins;k++)
03024         {
03025           if(k>j)
03026           {
03027             ele=Error4Expts->GetBinContent(j+1,k+1);
03028             ele+=(grid_bin_oscparerr[j]*grid_bin_oscparerr[k]*nexp_i->GetBinContent(j+1)*nexp_i->GetBinContent(k+1));
03029             Error4Expts->SetBinContent(j+1,k+1,ele);
03030             
03031             ele=Error4Expts->GetBinContent(k+1,j+1);
03032             ele+=(grid_bin_oscparerr[j]*grid_bin_oscparerr[k]*nexp_i->GetBinContent(j+1)*nexp_i->GetBinContent(k+1));
03033             Error4Expts->SetBinContent(k+1,j+1,ele);
03034             
03035             noff++;
03036           }
03037         }
03038       }
03039     }
03040     
03041     chi2hist_deltafit->Reset();
03042     chi2hist_deltafit->SetName(Form("Chi2Hist_DeltaFit_Inverted_%i",i));
03043     
03044     chi2hist_mhfit->Reset();
03045     chi2hist_mhfit->SetName(Form("Chi2Hist_MHFit_Inverted_%i",i));
03046     
03047     for(u=0;u<NumExpts;u++)
03048     {
03049       cout<<"expt "<<(u+1)<<"/"<<NumExpts<<endl;
03050       
03051       GenerateOneCorrelatedExp(nexp_i,Error4Expts);
03052       if(Print)
03053       {
03054         myfile << grid_sinsq2th13 << " " << grid_delta << " ";
03055         for(j=0;j<nBins;j++)
03056         {
03057           myfile << NObs->GetBinContent(j+1) << " ";
03058         }
03059       }
03060       
03061       chi2min_norm=GetMinLikelihood_Delta(true);
03062       chi2min_invt=GetMinLikelihood_Delta(false);
03063       
03064       ErrCalc->SetGridPred(nBins,nc_n,numucc_n,bnuecc_n,nutaucc_n,sig_n);
03065       ErrCalc->SetUseGrid(true);
03066       chi2_norm = 1e10;
03067       if(FitMethod==0)
03068       {
03069         chi2_norm = PoissonChi2(nexp_n);
03070       }
03071       else if(FitMethod==1)
03072       {
03073         chi2_norm = ScaledChi2(nexp_bkgd_n,nexp_signal_n);
03074       }
03075       else if(FitMethod==2)
03076       {
03077         chi2_norm = StandardChi2(nexp_n);
03078       }
03079       else if(FitMethod==3)
03080       {
03081         //Likelihood: "Standard" (N syst, N nuisance)
03082         //Calculate the likelihood (x2 for chi)
03083         Bkgd->Reset();
03084         Bkgd->Add(nexp_bkgd_n);
03085         Sig->Reset();
03086         Sig->Add(nexp_signal_n);
03087         chi2_norm = StandardLikelihood();
03088       }
03089       else if(FitMethod==4)
03090       {
03091         //Likelihood: Bin by Bin Calculation of Systematics
03092         //Calculate the likelihood (x2 for chi)
03093         Bkgd->Reset();
03094         Bkgd->Add(nexp_bkgd_n);
03095         Sig->Reset();
03096         Sig->Add(nexp_signal_n);
03097         chi2_norm = BinLikelihood();
03098       }
03099       else
03100       {
03101         cout<<"Error in RunMultiBinPseudoExpts_MHDeltaFit(): Unknown 'FitMethod'."<<endl;
03102       }
03103       
03104       ErrCalc->SetGridPred(nBins,nc_i,numucc_i,bnuecc_i,nutaucc_i,sig_i);
03105       ErrCalc->SetUseGrid(true);
03106       chi2_invt = 1e10;
03107       if(FitMethod==0)
03108       {
03109         chi2_invt = PoissonChi2(nexp_i);
03110       }
03111       else if(FitMethod==1)
03112       {
03113         chi2_invt = ScaledChi2(nexp_bkgd_i,nexp_signal_i);
03114       }
03115       else if(FitMethod==2)
03116       {
03117         chi2_invt = StandardChi2(nexp_i);
03118       }
03119       else if(FitMethod==3)
03120       {
03121         //Likelihood: "Standard" (N syst, N nuisance)
03122         //Calculate the likelihood (x2 for chi)
03123         Bkgd->Reset();
03124         Bkgd->Add(nexp_bkgd_i);
03125         Sig->Reset();
03126         Sig->Add(nexp_signal_i);
03127         chi2_invt = StandardLikelihood();
03128       }
03129       else if(FitMethod==4)
03130       {
03131         //Likelihood: Bin by Bin Calculation of Systematics
03132         //Calculate the likelihood (x2 for chi)
03133         Bkgd->Reset();
03134         Bkgd->Add(nexp_bkgd_i);
03135         Sig->Reset();
03136         Sig->Add(nexp_signal_i);
03137         chi2_invt = BinLikelihood();
03138       }
03139       else
03140       {
03141         cout<<"Error in RunMultiBinPseudoExpts_MHDeltaFit(): Unknown 'FitMethod'."<<endl;
03142       }
03143       
03144       if(Print)
03145       {
03146         myfile << chi2_norm << " " << chi2_invt << " " << chi2min_norm << " " << chi2min_invt << endl;
03147       }
03148       
03149       delchi2_norm = chi2_norm - chi2min_norm;
03150       delchi2_invt = chi2_invt - chi2min_invt;
03151       delchi2_bestmh = delchi2_invt;
03152       if(delchi2_norm<delchi2_bestmh) delchi2_bestmh = delchi2_norm;
03153       
03154       chi2hist_deltafit->Fill(delchi2_invt);
03155       chi2hist_mhfit->Fill(delchi2_invt - delchi2_bestmh);
03156     }
03157     f->cd();
03158     chi2hist_deltafit->Write();
03159     chi2hist_mhfit->Write();
03160     f->Close();
03161     f = new TFile(gSystem->ExpandPathName(outFileName.c_str()),"UPDATE");
03162   }
03163   
03164   if(Print) myfile.close();
03165   
03166   for(j=0;j<nBins;j++)
03167   {
03168     GridTree_Inverted[j]->Write();
03169     if (j<nBinsFHC){
03170       for(k=0;k<ExtrapFHC.size();k++)
03171       {
03172         GridTree_2_Inverted[j][k]->Write();
03173       }
03174     } else if (j>=nBinsFHC){
03175       for(k=0;k<ExtrapRHC.size();k++)
03176       {
03177         GridTree_2_Inverted[j][k]->Write();
03178       }
03179     }
03180   }
03181   
03182   f->Close();
03183   
03184   return;
03185 }

void NueFit2D_Joint::RunNSIDeltaChi2Contour ( int  cl = 0  ) 

Definition at line 4660 of file NueFit2D_Joint.cxx.

References NueFit2D::BinLikelihood(), NueFit2D::Bkgd, CombineFHCRHC(), DefineBinDlnLMinuit(), DefineStdDlnLMinuit(), NueFit2D::delta, NueFit2D::DeltaHigh, NueFit2D::DeltaLow, NueFit2D::epsetauHigh, NueFit2D::epsetauLow, NueFit2D::ErrorMatrix, ExtrapFHC, ExtrapRHC, NueFit2D::FitMethod, NueFit2D::FracErr_Bkgd, NueFit2D::FracErr_Sig, gSystem(), NueFit2D::InvErrorMatrix, NueFit2D::nBins, NueFit2D::nDeltaSteps, NueFit2D::nepsetauSteps, NueFit2D::NExp, NueFit2D::NObs, NueFit2D::outFileName, PhaseBroom, NueFit2D::PoissonChi2(), NueFit2D::ScaledChi2(), NueFit2D::Sig, NueFit2D::StandardChi2(), and NueFit2D::StandardLikelihood().

04661 {
04662   if(NObs==0)
04663   {
04664     cout<<"NObs not set.  Quitting..."<<endl;
04665     return;
04666   }
04667   if(FitMethod==1 && (FracErr_Bkgd==0 || FracErr_Sig==0))
04668   {
04669     cout<<"FracErr_Bkgd and FracErr_Sig need to be set for ScaledChi2.  Quitting..."<<endl;
04670     return;
04671   }
04672   if((ExtrapFHC.size()+ExtrapRHC.size())==0)
04673   {
04674     cout<<"No Extrapolate2D input.  Quitting..."<<endl;
04675     return;
04676   }
04677   
04678   if(FitMethod==3) DefineStdDlnLMinuit();
04679   if(FitMethod==4) DefineBinDlnLMinuit();
04680   
04681   Bkgd = (TH1D*)NObs->Clone("Bkgd");
04682   Bkgd->Reset();
04683   Sig = (TH1D*)NObs->Clone("Sig");
04684   Sig->Reset();
04685   NExp = (TH1D*)NObs->Clone("NExp");
04686   NExp->Reset();
04687   
04688   unsigned int ie;
04689   for(ie=0;ie<ExtrapFHC.size();ie++) ExtrapFHC[ie]->GetPrediction();
04690   for(ie=0;ie<ExtrapRHC.size();ie++) ExtrapRHC[ie]->GetPrediction();
04691   
04692   if(ErrorMatrix==0) ErrorMatrix = new TH2D("ErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
04693   if(InvErrorMatrix==0) InvErrorMatrix = new TH2D("InvErrorMatrix","",nBins,-0.5,nBins-0.5,nBins,-0.5,nBins-0.5);
04694   
04695   Int_t i,j,k,q,s;
04696   
04697   Double_t epsetau;
04698   Double_t delta = 0;
04699   Double_t Deltaincrement = 0;
04700   if(nDeltaSteps>0) Deltaincrement = (DeltaHigh - DeltaLow)/(nDeltaSteps);
04701   Double_t epsetauincrement = 0;
04702   if(nepsetauSteps>0) epsetauincrement = (epsetauHigh - epsetauLow)/(nepsetauSteps);
04703   
04704   double chi2[3],val[3];
04705   double ch[61],va[61];
04706   TGraph *g3;
04707   TF1* fit;
04708   double best;
04709   double vbt;
04710   double minchi2;
04711   double mc2;
04712   
04713   double limit;
04714   double delchi2;
04715   double sprev,delchi2prev;
04716   double contourlvl = 0;
04717   if(cl==0)//90% CL
04718   {
04719     contourlvl = 2.71;
04720   }
04721   else if(cl==1)//68.3% cL
04722   {
04723     contourlvl = 1.0;
04724   }
04725   else if(cl==2)//95% cL
04726   {
04727     contourlvl = 3.84;
04728   }
04729   else if(cl==3)//99% cL
04730   {
04731     contourlvl = 6.63;
04732   }
04733   else
04734   {
04735     cout<<"Error in RunNSIDeltaChi2Contour(): Input value should be 0,1,2,or 3 for 90%, 68.3%, 95%, or 99%.  Quitting..."<<endl;
04736     return;
04737   }
04738   
04739   cout<<"Seeking ";
04740   if(cl==0) cout<<"90% ";
04741   else if(cl==1)cout<<"68% ";
04742   else if(cl==2)cout<<"95% ";
04743   else if(cl==3)cout<<"99% ";
04744   else cout<< "??? ";
04745   cout<<" CL limits"<<endl;
04746   cout<<"PhaseBroom set to "<<PhaseBroom<<endl;
04747   
04748   vector<double> deltapts;
04749   vector<double> bestfit_norm;
04750   vector<double> bestfit_invt;
04751   vector<double> limit_norm;
04752   vector<double> limit_invt;
04753   vector<double> limit_normlow;
04754   vector<double> limit_invtlow;
04755   
04756   for(j=0;j<nDeltaSteps+1;j++)
04757   {
04758     delta = j*Deltaincrement*TMath::Pi() + DeltaLow;
04759     for(ie=0;ie<ExtrapFHC.size();ie++){
04760       ExtrapFHC[ie]->SetDeltaCP((1.0-PhaseBroom)*delta);
04761       ExtrapFHC[ie]->SetDelta_etau(PhaseBroom*delta);
04762     }
04763     for(ie=0;ie<ExtrapRHC.size();ie++){
04764       ExtrapRHC[ie]->SetDeltaCP((1.0-PhaseBroom)*delta);
04765       ExtrapRHC[ie]->SetDelta_etau(PhaseBroom*delta);
04766     }
04767     deltapts.push_back(delta/TMath::Pi());
04768     
04769     best=-1.;
04770     minchi2=100;
04771     mc2=1000;
04772     for(i=0;i<3;i++)
04773     {
04774       chi2[i]=-1.;
04775       val[i]=-1.;
04776     }
04777     for(k=0;k<=60;k++){
04778       ch[k]=-1.;
04779       va[k]=-1.;
04780     }
04781     for(q=0;q<=60;q++){
04782       epsetau = epsetauLow+q*0.1;
04783       for(ie=0;ie<ExtrapFHC.size();ie++)
04784       {
04785         ExtrapFHC[ie]->SetEps_etau(epsetau);
04786         ExtrapFHC[ie]->OscillatePrediction();
04787       }
04788       for(ie=0;ie<ExtrapRHC.size();ie++)
04789       {
04790         ExtrapRHC[ie]->SetEps_etau(epsetau);
04791         ExtrapRHC[ie]->OscillatePrediction();
04792       }
04793 
04794       Bkgd->Reset();
04795       Sig->Reset();
04796       NExp->Reset();
04797 
04798       CombineFHCRHC();
04799       NExp->Add(Bkgd);
04800       NExp->Add(Sig);
04801       
04802       if(FitMethod==0)
04803       {
04804         ch[q] = PoissonChi2(NExp);
04805       }
04806       else if(FitMethod==1)
04807       {
04808         ch[q] = ScaledChi2(Bkgd,Sig);
04809       }
04810       else if(FitMethod==2)
04811       {
04812         ch[q] = StandardChi2(NExp);
04813       }
04814       else if(FitMethod==3)
04815       {
04816         //Likelihood: "Standard" (N syst, N nuisance)
04817         //Calculate the likelihood (x2 for chi)
04818         ch[q] = StandardLikelihood();
04819       }
04820       else if(FitMethod==4)
04821       {
04822         //Likelihood: Bin by Bin Calculation of Systematics
04823         //Calculate the likelihood (x2 for chi)
04824         ch[q] = BinLikelihood();
04825       }
04826       else
04827       {
04828         cout<<"Error in RunNSIDeltaChi2Contour(): Unknown 'FitMethod'."<<endl;
04829       }
04830      
04831       va[q]=epsetau;
04832     }//coarse sweep
04833 
04834     for(s=0;s<61;s++){
04835       if(ch[s]<mc2){
04836         mc2=ch[s];
04837         vbt=va[s];
04838       }
04839     }//find junk
04840     cout << "Found CBF for Delta: "<<delta<< " at " <<vbt<< " with chi2 of "<<mc2<< endl;
04841     for(i=0;i<=nepsetauSteps/5;i++)
04842     {
04843       chi2[0] = chi2[1];
04844       chi2[1] = chi2[2];
04845       val[0] = val[1];
04846       val[1] = val[2];
04847       epsetau = vbt + (-(nepsetauSteps/10)+i)*epsetauincrement;
04848  
04849       for(ie=0;ie<ExtrapFHC.size();ie++)
04850       {
04851         ExtrapFHC[ie]->SetEps_etau(epsetau);
04852         ExtrapFHC[ie]->OscillatePrediction();
04853       }
04854       for(ie=0;ie<ExtrapRHC.size();ie++)
04855       {
04856         ExtrapRHC[ie]->SetEps_etau(epsetau);
04857         ExtrapRHC[ie]->OscillatePrediction();
04858       }
04859 
04860       Bkgd->Reset();
04861       Sig->Reset();
04862       NExp->Reset();
04863 
04864       CombineFHCRHC();
04865       NExp->Add(Bkgd);
04866       NExp->Add(Sig);
04867       
04868       chi2[2] = 1e10;
04869       if(FitMethod==0)
04870       {
04871         chi2[2] = PoissonChi2(NExp);
04872       }
04873       else if(FitMethod==1)
04874       {
04875         chi2[2] = ScaledChi2(Bkgd,Sig);
04876       }
04877       else if(FitMethod==2)
04878       {
04879         chi2[2] = StandardChi2(NExp);
04880       }
04881       else if(FitMethod==3)
04882       {
04883         //Likelihood: "Standard" (N syst, N nuisance)
04884         //Calculate the likelihood (x2 for chi)
04885         chi2[2] = StandardLikelihood();
04886       }
04887       else if(FitMethod==4)
04888       {
04889         //Likelihood: Bin by Bin Calculation of Systematics
04890         //Calculate the likelihood (x2 for chi)
04891         chi2[2] = BinLikelihood();
04892       }
04893       else
04894       {
04895         cout<<"Error in RunNSIDeltaChi2Contour(): Unknown 'FitMethod'."<<endl;
04896       }
04897       
04898       val[2] = epsetau;
04899       
04900       if(i<2) continue;
04901       
04902       if(i<3 && chi2[2]>chi2[1] && chi2[1]>chi2[0])//first three points are increasing, first point is minimum.
04903       {
04904         best = val[0];
04905         minchi2 = chi2[0];
04906         cout<<"minimum at 1st point: "<<minchi2<<" at "<<best<<endl;
04907         break;
04908       }
04909       
04910       if(chi2[2]>chi2[1] && chi2[0]>chi2[1])//found minimum
04911       {
04912         g3 = new TGraph(3, val, chi2);
04913         fit = new TF1("pol2", "pol2");
04914         g3->Fit(fit, "Q");//fit to second order polynominal
04915         if(fit->GetParameter(2) > 0)//if the x^2 term is nonzero
04916         {
04917           best = -fit->GetParameter(1)/(2*fit->GetParameter(2));//the location of the minimum is -p1/(2*p2)
04918           minchi2 = fit->GetParameter(0) + fit->GetParameter(1)*best + fit->GetParameter(2)*best*best;
04919           cout<<"minimum with fit: "<<minchi2<<" at "<<best<<endl;
04920         }
04921         else//if the x^2 term is zero, then just use the minimum you got by scanning
04922         {
04923           best = val[1];
04924           minchi2 = chi2[1];
04925           cout<<"minimum with scan: "<<minchi2<<" at "<<best<<endl;
04926         }
04927         break;
04928       }
04929     }
04930     bestfit_norm.push_back(best);
04931     
04932     limit = -1.;
04933     delchi2prev = 1000;
04934     sprev = best;
04935     for(i=0;i<(epsetauHigh-best)/epsetauincrement;i++)
04936     {
04937       //epsetau = i*epsetauincrement + epsetauLow;
04938       epsetau = i*epsetauincrement + best;
04939       for(ie=0;ie<ExtrapFHC.size();ie++)
04940       {
04941         ExtrapFHC[ie]->SetEps_etau(epsetau);
04942         ExtrapFHC[ie]->OscillatePrediction();
04943       }
04944       for(ie=0;ie<ExtrapRHC.size();ie++)
04945       {
04946         ExtrapRHC[ie]->SetEps_etau(epsetau);
04947         ExtrapRHC[ie]->OscillatePrediction();
04948       }
04949 
04950       Bkgd->Reset();
04951       Sig->Reset();
04952       NExp->Reset();
04953 
04954       CombineFHCRHC();
04955       NExp->Add(Bkgd);
04956       NExp->Add(Sig);
04957       
04958       delchi2 = 1e10;
04959       if(FitMethod==0)
04960       {
04961         delchi2 = PoissonChi2(NExp) - minchi2;
04962       }
04963       else if(FitMethod==1)
04964       {
04965         delchi2 = ScaledChi2(Bkgd,Sig) - minchi2;
04966       }
04967       else if(FitMethod==2)
04968       {
04969         delchi2 = StandardChi2(NExp) - minchi2;
04970       }
04971       else if(FitMethod==3)
04972       {
04973         //Likelihood: "Standard" (N syst, N nuisance)
04974         //Calculate the likelihood (x2 for chi)
04975         delchi2 = StandardLikelihood() - minchi2;
04976       }
04977       else if(FitMethod==4)
04978       {
04979         //Likelihood: Bin by Bin Calculation of Systematics
04980         //Calculate the likelihood (x2 for chi)
04981         delchi2 = BinLikelihood() - minchi2;
04982       }
04983       else
04984       {
04985         cout<<"Error in RunNSIDeltaChi2Contour(): Unknown 'FitMethod'."<<endl;
04986       }
04987       
04988       if(i==1) continue;
04989       
04990       if(delchi2>contourlvl && delchi2prev<contourlvl)
04991       {
04992         limit = epsetau + ((epsetau-sprev)/(delchi2 - delchi2prev))*(contourlvl - delchi2);
04993         cout<<delchi2prev<<", "<<sprev<<", "<<delchi2<<", "<<epsetau<<", "<<limit<<endl;
04994         break;
04995       }
04996       delchi2prev = delchi2;
04997       sprev = epsetau;
04998     }
04999     limit_norm.push_back(limit);
05000 
05001     limit = -1.;
05002     delchi2prev = 1000;
05003     sprev = best;
05004     for(i=0;i<(best-epsetauLow)/epsetauincrement;i++)
05005     {
05006       //epsetau = i*epsetauincrement + epsetauLow;
05007       epsetau = -i*epsetauincrement + best;
05008       for(ie=0;ie<ExtrapFHC.size();ie++)
05009       {
05010         ExtrapFHC[ie]->SetEps_etau(epsetau);
05011         ExtrapFHC[ie]->OscillatePrediction();
05012       }
05013       for(ie=0;ie<ExtrapRHC.size();ie++)
05014       {
05015         ExtrapRHC[ie]->SetEps_etau(epsetau);
05016         ExtrapRHC[ie]->OscillatePrediction();
05017       }
05018 
05019       Bkgd->Reset();
05020       Sig->Reset();
05021       NExp->Reset();
05022 
05023       CombineFHCRHC();
05024       NExp->Add(Bkgd);
05025       NExp->Add(Sig);
05026       
05027       delchi2 = 1e10;
05028       if(FitMethod==0)
05029       {
05030         delchi2 = PoissonChi2(NExp) - minchi2;
05031       }
05032       else if(FitMethod==1)
05033       {
05034         delchi2 = ScaledChi2(Bkgd,Sig) - minchi2;
05035       }
05036       else if(FitMethod==2)
05037       {
05038         delchi2 = StandardChi2(NExp) - minchi2;
05039       }
05040       else if(FitMethod==3)
05041       {
05042         //Likelihood: "Standard" (N syst, N nuisance)
05043         //Calculate the likelihood (x2 for chi)
05044         delchi2 = StandardLikelihood() - minchi2;
05045       }
05046       else if(FitMethod==4)
05047       {
05048         //Likelihood: Bin by Bin Calculation of Systematics
05049         //Calculate the likelihood (x2 for chi)
05050         delchi2 = BinLikelihood() - minchi2;
05051       }
05052       else
05053       {
05054         cout<<"Error in RunNSIDeltaChi2Contour(): Unknown 'FitMethod'."<<endl;
05055       }
05056       
05057       if(i==1) continue;
05058       
05059       if(delchi2>contourlvl && delchi2prev<contourlvl)
05060       {
05061         limit = epsetau + ((epsetau-sprev)/(delchi2 - delchi2prev))*(contourlvl - delchi2);
05062         cout<<delchi2prev<<", "<<sprev<<", "<<delchi2<<", "<<epsetau<<", "<<limit<<endl;
05063         break;
05064       }
05065       delchi2prev = delchi2;
05066       sprev = epsetau;
05067     }
05068     limit_normlow.push_back(limit);
05069   }
05070   
05071   for(ie=0;ie<ExtrapFHC.size();ie++) ExtrapFHC[ie]->InvertMassHierarchy();
05072   for(ie=0;ie<ExtrapRHC.size();ie++) ExtrapRHC[ie]->InvertMassHierarchy();
05073 
05074   
05075   for(j=0;j<nDeltaSteps+1;j++)
05076   {
05077     delta = j*Deltaincrement*TMath::Pi() + DeltaLow;
05078     for(ie=0;ie<ExtrapFHC.size();ie++){
05079       ExtrapFHC[ie]->SetDeltaCP((1.0-PhaseBroom)*delta);
05080       ExtrapFHC[ie]->SetDelta_etau(PhaseBroom*delta);
05081     }
05082     for(ie=0;ie<ExtrapRHC.size();ie++){
05083       ExtrapRHC[ie]->SetDeltaCP((1.0-PhaseBroom)*delta);
05084       ExtrapRHC[ie]->SetDelta_etau(PhaseBroom*delta);
05085     }
05086 
05087     best=-1.;
05088     minchi2=100;
05089     mc2=1000;
05090     for(i=0;i<3;i++)
05091     {
05092       chi2[i]=-1.;
05093       val[i]=-1.;
05094     }
05095     for(k=0;k<=60;k++){
05096       ch[k]=-1.;
05097       va[k]=-1.;
05098     }
05099     for(q=0;q<=60;q++){
05100       epsetau = epsetauLow+q*0.1;
05101       for(ie=0;ie<ExtrapFHC.size();ie++)
05102       {
05103         ExtrapFHC[ie]->SetEps_etau(epsetau);
05104         ExtrapFHC[ie]->OscillatePrediction();
05105       }
05106       for(ie=0;ie<ExtrapRHC.size();ie++)
05107       {
05108         ExtrapRHC[ie]->SetEps_etau(epsetau);
05109         ExtrapRHC[ie]->OscillatePrediction();
05110       }
05111 
05112       Bkgd->Reset();
05113       Sig->Reset();
05114       NExp->Reset();
05115 
05116       CombineFHCRHC();
05117       NExp->Add(Bkgd);
05118       NExp->Add(Sig);
05119       
05120       if(FitMethod==0)
05121       {
05122         ch[q] = PoissonChi2(NExp);
05123       }
05124       else if(FitMethod==1)
05125       {
05126         ch[q] = ScaledChi2(Bkgd,Sig);
05127       }
05128       else if(FitMethod==2)
05129       {
05130         ch[q] = StandardChi2(NExp);
05131       }
05132       else if(FitMethod==3)
05133       {
05134         //Likelihood: "Standard" (N syst, N nuisance)
05135         //Calculate the likelihood (x2 for chi)
05136         ch[q] = StandardLikelihood();
05137       }
05138       else if(FitMethod==4)
05139       {
05140         //Likelihood: Bin by Bin Calculation of Systematics
05141         //Calculate the likelihood (x2 for chi)
05142         ch[q] = BinLikelihood();
05143       }
05144       else
05145       {
05146         cout<<"Error in RunNSIDeltaChi2Contour(): Unknown 'FitMethod'."<<endl;
05147       }
05148      
05149       va[q]=epsetau;
05150     }//coarse sweep
05151 
05152     for(s=0;s<61;s++){
05153       if(ch[s]<mc2){
05154         mc2=ch[s];
05155         vbt=va[s];
05156       }
05157     }//find junk
05158     cout << "Found CBF for Delta: "<<delta<< " at " <<vbt<< " with chi2 of "<<mc2<< endl;
05159     for(i=0;i<=nepsetauSteps/5;i++)
05160     {
05161       chi2[0] = chi2[1];
05162       chi2[1] = chi2[2];
05163       val[0] = val[1];
05164       val[1] = val[2];
05165       epsetau = vbt + (-(nepsetauSteps/10)+i)*epsetauincrement;
05166 
05167       for(ie=0;ie<ExtrapFHC.size();ie++)
05168       {
05169         ExtrapFHC[ie]->SetEps_etau(epsetau);
05170         ExtrapFHC[ie]->OscillatePrediction();
05171       }
05172 
05173       for(ie=0;ie<ExtrapRHC.size();ie++)
05174       {
05175         ExtrapRHC[ie]->SetEps_etau(epsetau);
05176         ExtrapRHC[ie]->OscillatePrediction();
05177       }
05178 
05179       Bkgd->Reset();
05180       Sig->Reset();
05181       NExp->Reset();
05182 
05183       CombineFHCRHC();
05184       NExp->Add(Bkgd);
05185       NExp->Add(Sig);
05186       
05187       chi2[2] = 1e10;
05188       if(FitMethod==0)
05189       {
05190         chi2[2] = PoissonChi2(NExp);
05191       }
05192       else if(FitMethod==1)
05193       {
05194         chi2[2] = ScaledChi2(Bkgd,Sig);
05195       }
05196       else if(FitMethod==2)
05197       {
05198         chi2[2] = StandardChi2(NExp);
05199       }
05200       else if(FitMethod==3)
05201       {
05202         //Likelihood: "Standard" (N syst, N nuisance)
05203         //Calculate the likelihood (x2 for chi)
05204         chi2[2] = StandardLikelihood();
05205       }
05206       else if(FitMethod==4)
05207       {
05208         //Likelihood: Bin by Bin Calculation of Systematics
05209         //Calculate the likelihood (x2 for chi)
05210         chi2[2] = BinLikelihood();
05211       }
05212       else
05213       {
05214         cout<<"Error in RunDeltaChi2Contour(): Unknown 'FitMethod'."<<endl;
05215       }
05216       
05217       val[2] = epsetau;
05218       
05219       if(i<2) continue;
05220       
05221       if(i<3 && chi2[2]>chi2[1] && chi2[1]>chi2[0])//first three points are increasing, first point is minimum.
05222       {
05223         best = val[0];
05224         minchi2 = chi2[0];
05225         cout<<"minimum at 1st point: "<<minchi2<<" at "<<best<<endl;
05226         break;
05227       }
05228       
05229       if(chi2[2]>chi2[1] && chi2[0]>chi2[1])//found minimum
05230       {
05231         g3 = new TGraph(3, val, chi2);
05232         fit = new TF1("pol2", "pol2");
05233         g3->Fit(fit, "Q");//fit to second order polynominal
05234         if(fit->GetParameter(2) > 0)//if the x^2 term is nonzero
05235         {
05236           best = -fit->GetParameter(1)/(2*fit->GetParameter(2));//the location of the minimum is -p1/(2*p2)
05237           minchi2 = fit->GetParameter(0) + fit->GetParameter(1)*best + fit->GetParameter(2)*best*best;
05238           cout<<"minimum with fit: "<<minchi2<<" at "<<best<<endl;
05239         }
05240         else//if the x^2 term is zero, then just use the minimum you got by scanning
05241         {
05242           best = val[1];
05243           minchi2 = chi2[1];
05244           cout<<"minimum with scan: "<<minchi2<<" at "<<best<<endl;
05245         }
05246         break;
05247       }
05248     }
05249     bestfit_invt.push_back(best);
05250     
05251     limit = -1.;
05252     delchi2prev = 1000;
05253     sprev = best;
05254     for(i=0;i<(epsetauHigh-best)/epsetauincrement;i++)
05255     {
05256       //epsetau = i*epsetauincrement + epsetauLow;
05257       epsetau = i*epsetauincrement + best;
05258       for(ie=0;ie<ExtrapFHC.size();ie++)
05259       {
05260         ExtrapFHC[ie]->SetEps_etau(epsetau);
05261         ExtrapFHC[ie]->OscillatePrediction();
05262       }
05263       for(ie=0;ie<ExtrapRHC.size();ie++)
05264       {
05265         ExtrapRHC[ie]->SetEps_etau(epsetau);
05266         ExtrapRHC[ie]->OscillatePrediction();
05267       }
05268 
05269       Bkgd->Reset();
05270       Sig->Reset();
05271       NExp->Reset();
05272 
05273       CombineFHCRHC();
05274       NExp->Add(Bkgd);
05275       NExp->Add(Sig);
05276       
05277       delchi2 = 1e10;
05278       if(FitMethod==0)
05279       {
05280         delchi2 = PoissonChi2(NExp) - minchi2;
05281       }
05282       else if(FitMethod==1)
05283       {
05284         delchi2 = ScaledChi2(Bkgd,Sig) - minchi2;
05285       }
05286       else if(FitMethod==2)
05287       {
05288         delchi2 = StandardChi2(NExp) - minchi2;
05289       }
05290       else if(FitMethod==3)
05291       {
05292         //Likelihood: "Standard" (N syst, N nuisance)
05293         //Calculate the likelihood (x2 for chi)
05294         delchi2 = StandardLikelihood() - minchi2;
05295       }
05296       else if(FitMethod==4)
05297       {
05298         //Likelihood: Bin by Bin Calculation of Systematics
05299         //Calculate the likelihood (x2 for chi)
05300         delchi2 = BinLikelihood() - minchi2;
05301       }
05302       else
05303       {
05304         cout<<"Error in RunNSIDeltaChi2Contour(): Unknown 'FitMethod'."<<endl;
05305       }
05306       
05307       if(i==1) continue;
05308       
05309       if(delchi2>contourlvl && delchi2prev<contourlvl)
05310       {
05311         limit = epsetau + ((epsetau-sprev)/(delchi2 - delchi2prev))*(contourlvl - delchi2);
05312         cout<<delchi2prev<<", "<<sprev<<", "<<delchi2<<", "<<epsetau<<", "<<limit<<endl;
05313         break;
05314       }
05315       delchi2prev = delchi2;
05316       sprev = epsetau;
05317     }
05318     limit_invt.push_back(limit);
05319 
05320 
05321     limit = -1.;
05322     delchi2prev = 1000;
05323     sprev = best;
05324     for(i=0;i<(best-epsetauLow)/epsetauincrement;i++)
05325     {
05326       //epsetau = i*epsetauincrement + epsetauLow;
05327       epsetau = -i*epsetauincrement + best;
05328       for(ie=0;ie<ExtrapFHC.size();ie++)
05329       {
05330         ExtrapFHC[ie]->SetEps_etau(epsetau);
05331         ExtrapFHC[ie]->OscillatePrediction();
05332       }
05333       for(ie=0;ie<ExtrapRHC.size();ie++)
05334       {
05335         ExtrapRHC[ie]->SetEps_etau(epsetau);
05336         ExtrapRHC[ie]->OscillatePrediction();
05337       }
05338 
05339       Bkgd->Reset();
05340       Sig->Reset();
05341       NExp->Reset();
05342 
05343       CombineFHCRHC();
05344       NExp->Add(Bkgd);
05345       NExp->Add(Sig);
05346       
05347       delchi2 = 1e10;
05348       if(FitMethod==0)
05349       {
05350         delchi2 = PoissonChi2(NExp) - minchi2;
05351       }
05352       else if(FitMethod==1)
05353       {
05354         delchi2 = ScaledChi2(Bkgd,Sig) - minchi2;
05355       }
05356       else if(FitMethod==2)
05357       {
05358         delchi2 = StandardChi2(NExp) - minchi2;
05359       }
05360       else if(FitMethod==3)
05361       {
05362         //Likelihood: "Standard" (N syst, N nuisance)
05363         //Calculate the likelihood (x2 for chi)
05364         delchi2 = StandardLikelihood() - minchi2;
05365       }
05366       else if(FitMethod==4)
05367       {
05368         //Likelihood: Bin by Bin Calculation of Systematics
05369         //Calculate the likelihood (x2 for chi)
05370         delchi2 = BinLikelihood() - minchi2;
05371       }
05372       else
05373       {
05374         cout<<"Error in RunNSIDeltaChi2Contour(): Unknown 'FitMethod'."<<endl;
05375       }
05376       
05377       if(i==1) continue;
05378       
05379       if(delchi2>contourlvl && delchi2prev<contourlvl)
05380       {
05381         limit = epsetau + ((epsetau-sprev)/(delchi2 - delchi2prev))*(contourlvl - delchi2);
05382         cout<<delchi2prev<<", "<<sprev<<", "<<delchi2<<", "<<epsetau<<", "<<limit<<endl;
05383         break;
05384       }
05385       delchi2prev = delchi2;
05386       sprev = epsetau;
05387     }
05388     limit_invtlow.push_back(limit);
05389   }
05390   
05391   double *s_best_n = new double[nDeltaSteps+1];
05392   double *s_best_i = new double[nDeltaSteps+1];
05393   double *s_limit_n = new double[nDeltaSteps+1];
05394   double *s_limit_i = new double[nDeltaSteps+1];
05395   double *s_limit_nlow = new double[nDeltaSteps+1];
05396   double *s_limit_ilow = new double[nDeltaSteps+1];
05397   double *d = new double[nDeltaSteps+1];
05398   for(i=0;i<nDeltaSteps+1;i++)
05399   {
05400     d[i] = deltapts.at(i);
05401     s_best_n[i] = bestfit_norm.at(i);
05402     s_best_i[i] = bestfit_invt.at(i);
05403     s_limit_n[i] = limit_norm.at(i);
05404     s_limit_i[i] = limit_invt.at(i);
05405     
05406     s_limit_nlow[i] = limit_normlow.at(i);
05407     s_limit_ilow[i] = limit_invtlow.at(i);
05408   }
05409   
05410   TGraph *gn = new TGraph(nDeltaSteps+1,s_best_n,d);
05411   gn->SetMarkerStyle(20);
05412   gn->SetTitle("");
05413   gn->GetYaxis()->SetTitle("#delta");
05414   gn->GetXaxis()->SetTitle("#epsilon_{e#tau}");
05415   gn->GetXaxis()->SetLimits(-3.0,3.0);
05416   gn->SetLineWidth(4);
05417   gn->SetMaximum(2);
05418   gn->SetName("BestFit_Normal");
05419   
05420   TGraph *gi = new TGraph(nDeltaSteps+1,s_best_i,d);
05421   gi->SetMarkerStyle(20);
05422   gi->SetTitle("");
05423   gi->GetYaxis()->SetTitle("#delta");
05424   gi->GetXaxis()->SetTitle("#epsilon_{e#tau}");
05425   gi->GetXaxis()->SetLimits(-3.0,3.0);
05426   gi->SetLineWidth(4);
05427   gi->SetLineStyle(2);
05428   gi->SetMaximum(2);
05429   gi->SetName("BestFit_Inverted");
05430   
05431   TGraph *gn_limit = new TGraph(nDeltaSteps+1,s_limit_n,d);
05432   gn_limit->SetMarkerStyle(20);
05433   gn_limit->SetTitle("");
05434   gn_limit->GetYaxis()->SetTitle("#delta");
05435   gn_limit->GetXaxis()->SetTitle("#epsilon_{e#tau}");
05436   gn_limit->GetXaxis()->SetLimits(-3.0,3.0);
05437   gn_limit->SetLineWidth(4);
05438   gn_limit->SetLineColor(kBlue);
05439   gn_limit->SetMarkerColor(kBlue);
05440   gn_limit->SetMaximum(2);
05441   gn_limit->SetName("Limit_Normal");
05442 
05443  TGraph *gn_limitlow = new TGraph(nDeltaSteps+1,s_limit_nlow,d);
05444   gn_limitlow->SetMarkerStyle(20);
05445   gn_limitlow->SetTitle("");
05446   gn_limitlow->GetYaxis()->SetTitle("#delta");
05447   gn_limitlow->GetXaxis()->SetTitle("#epsilon_{e#tau}");
05448   gn_limitlow->GetXaxis()->SetLimits(-3.0,3.0);
05449   gn_limitlow->SetLineWidth(4);
05450   gn_limitlow->SetLineColor(kBlue);
05451   gn_limitlow->SetMarkerColor(kBlue);
05452   gn_limitlow->SetMaximum(2);
05453   gn_limitlow->SetName("Limit_NormalLow");
05454   
05455   TGraph *gi_limit = new TGraph(nDeltaSteps+1,s_limit_i,d);
05456   gi_limit->SetMarkerStyle(20);
05457   gi_limit->SetTitle("");
05458   gi_limit->GetYaxis()->SetTitle("#delta");
05459   gi_limit->GetXaxis()->SetTitle("#epsilon_{e#tau}");
05460   gi_limit->GetXaxis()->SetLimits(-3.0,3.0);
05461   gi_limit->SetLineWidth(4);
05462   gi_limit->SetLineColor(kRed);
05463   gi_limit->SetMarkerColor(kRed);
05464   gi_limit->SetMaximum(2);
05465   gi_limit->SetName("Limit_Inverted");
05466  
05467  TGraph *gi_limitlow = new TGraph(nDeltaSteps+1,s_limit_ilow,d);
05468   gi_limitlow->SetMarkerStyle(20);
05469   gi_limitlow->SetTitle("");
05470   gi_limitlow->GetYaxis()->SetTitle("#delta");
05471   gi_limitlow->GetXaxis()->SetTitle("#epsilon_{e#tau}");
05472   gi_limitlow->GetXaxis()->SetLimits(-3.0,3.0);
05473   gi_limitlow->SetLineWidth(4);
05474   gi_limitlow->SetLineColor(kRed);
05475   gi_limitlow->SetMarkerColor(kRed);
05476   gi_limitlow->SetMaximum(2);
05477   gi_limitlow->SetName("Limit_InvertedLow");
05478  
05479   if(cl==0) cout<<"90% ";
05480   if(cl==1) cout<<"68% ";
05481   if(cl==2) cout<<"95% ";
05482   if(cl==3) cout<<"99% ";
05483   cout<<"confidence level limit = "<<limit_norm.at(0)<<", "<<limit_invt.at(0)<<endl;
05484   
05485   TFile *fout = new TFile(gSystem->ExpandPathName(outFileName.c_str()),"RECREATE");
05486   gn->Write();
05487   gi->Write();
05488   gn_limit->Write();
05489   gi_limit->Write();
05490   gn_limitlow->Write();
05491   gi_limitlow->Write();
05492   fout->Write();
05493   fout->Close();
05494   
05495   delete [] s_best_n;
05496   delete [] s_best_i;
05497   delete [] s_limit_n;
05498   delete [] s_limit_i;
05499   delete [] s_limit_nlow;
05500   delete [] s_limit_ilow;
05501   delete [] d;
05502   
05503   return;
05504 }

void NueFit2D_Joint::ScalePredictionFHC ( double  scale = 1.0  )  [inline]

Definition at line 124 of file NueFit2D_Joint.h.

References potScaleFHC.

00124 { potScaleFHC=scale; };

void NueFit2D_Joint::ScalePredictionRHC ( double  scale = 1.0  )  [inline]

Definition at line 125 of file NueFit2D_Joint.h.

References potScaleRHC.

00125 { potScaleRHC=scale; };

void NueFit2D_Joint::SetPhaseBroom ( double  pb = 0.0  )  [inline]

Definition at line 155 of file NueFit2D_Joint.h.

References PhaseBroom.

00155                                      { 
00156       if(pb>=0.0 && pb <=1.0){ 
00157        PhaseBroom = pb; //0 cycles delta_CP, 1 cylces delta_etau, inbetween cycles by fraction
00158       }
00159       else{
00160        PhaseBroom = 0.0;
00161       }
00162     }


Member Data Documentation

unsigned int NueFit2D_Joint::nBinsFHC [private]
unsigned int NueFit2D_Joint::nBinsRHC [private]
int NueFit2D_Joint::NHIH [private]

Definition at line 170 of file NueFit2D_Joint.h.

Referenced by DoIH(), DoNSIFitForever(), and GetNSIFFX2().

double NueFit2D_Joint::PhaseBroom [private]

Definition at line 169 of file NueFit2D_Joint.h.

Referenced by RunNSIDeltaChi2Contour(), and SetPhaseBroom().

double NueFit2D_Joint::potScaleFHC [private]
double NueFit2D_Joint::potScaleRHC [private]

Definition at line 177 of file NueFit2D_Joint.h.

Referenced by CalculateDlnLMatrix(), and NueFit2D_Joint().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1