NuMatrixFitter Class Reference

#include <NuMatrixFitter.h>

List of all members.

Public Types

enum  NuModel_t { kUndefined = 0, kOscillation = 1, kDecay = 2, kDecoherence = 3 }

Public Member Functions

 NuMatrixFitter ()
 NuMatrixFitter (const Double_t FluxPoT, const string ndData, const string fdData, const string helperFile, const string xsecFile, const NuXMLConfig *xmlConfig=0)
 NuMatrixFitter (const Double_t fluxPoT, const std::vector< std::string > ndDataFiles, const std::vector< std::string > fdDataFiles, const std::vector< std::string > helperFiles, const std::string xsecFile, const NuXMLConfig *xmlConfig=0)
virtual ~NuMatrixFitter ()
virtual Double_t CalculateChi2 (const TH1D *fdPred, const TH1D *fdData) const
virtual Double_t CalculateLikelihood (const TH1D *fdPred, const TH1D *fdData) const
virtual Double_t CalculateLikelihoodNormPenalty (const TH1D *fdPred, const TH1D *fdData, const Double_t normalisation) const
virtual Double_t CalculateLikelihood (const TH1D *fdPred, const TH1D *fdData, const Double_t nubarEcut) const
Double_t NormalisationPenaltyTerm (const Double_t normalisation) const
TH1D RebinForFit (const TH1D *fineHist)
virtual void CCAnalysis ()
virtual void CPTAnalysis ()
virtual void AppearanceAnalysis ()
virtual NuMatrixOutputDoCPTFit (NuMatrixInput *mmInput, TH1D *numubarfdData, const Double_t dm2, const Double_t sn2)
virtual NuMatrixOutputDoNoChargeCutFit (NuMatrixInput &mmInput, TH1D *numubarfdData, const NuModel_t model=kOscillation)
virtual NuMatrixOutputDoCCFitChargeCut (const NuMatrixInput &mmInput, const TH1D *numufdData, const TH1D *numubarfdData, const Double_t nubarEcut=-1.0, const NuModel_t model=kOscillation)
virtual NuMatrixOutputDoPRLCCFit (const NuMatrixInput &mmInput, const TH1D *fdData, const NuModel_t model=kOscillation)
virtual NuMatrixOutputDoTransitionFit (const char *inputFilename, const TH1D *numubarfdData, const Double_t dm2, const Double_t sn2) const
virtual void MultiRunCPTFit (const std::vector< std::string > inputFilenames, const std::vector< std::string > fdDataFilenames, const Double_t dm2, const Double_t sn2, const std::vector< std::string > outputFilenames, const std::string combinedOutputFilename)
virtual void MultiRunPRLCCFit (const vector< string > inputFilenames, const vector< string > fdDataFilenames, const std::vector< std::string > outputFilenames, const std::string combinedOutputFilename, const NuModel_t model=kOscillation)
virtual const std::pair
< std::vector< NuMatrixOutput * >
, NuMatrixOutput * > 
DoMultiRunCPTFit (const vector< NuMatrixInput * > mmInput, const vector< TH1D > numubarfdData, const Double_t dm2, const Double_t sn2) const
virtual const std::pair
< std::vector< NuMatrixOutput * >
, NuMatrixOutput * > 
DoMultiRunPRLCCFit (const vector< NuMatrixInput * > mmInput, const vector< TH1D > numufdData, const NuModel_t model=kOscillation)
virtual const NuMatrixOutputDoMultiRunTransitionFit (const vector< NuMatrixInput > mmInput, const vector< TH1D > numubarfdData, const Double_t dm2, const Double_t sn2) const
virtual void CPTFit (const std::string inputFilename, const string fdDataFilename, const Double_t dm2, const Double_t sn2, const std::string outputFilename)
virtual void CCFitChargeCut (const std::string inputFilename, const string fdDataFilename, const std::string outputFilename, const Double_t nubarEcut=-1.0, const NuModel_t model=kOscillation)
virtual void NoChargeCutFit (const std::string inputFilename, const string fdDataFilename, const std::string outputFilename, const NuModel_t model=kOscillation)
void PRLCCFit (const string inputFilename, const string fdDataFilename, const string outputFilename, const NuModel_t model=kOscillation)
void TransitionFit (const string inputFilename, const string fdDataFilename, const string outputFilename, const Int_t experiments)
void WriteNoChargeCutHistosForFitter (std::string outputFileName="")
void WriteHistosForFitter (TString outputFileName="$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutputForFitter.root")
virtual void WritePRLCCHistosForFitter (std::string outputFileName)
virtual void WriteMultiRunHistosForFitter (const vector< TString > outputFileNames)
virtual void NuMuBarAppearanceAnalysis ()
void SetXMLConfig (const NuXMLConfig *xmlConfig)

Private Member Functions

void FillGridSearchParams (const NuXMLConfig *xmlConfig)
void SetGridParamsDefault ()

Private Attributes

Double_t ffluxPoT
std::string fndDataFilename
std::string ffdDataFilename
std::string fhelperFilename
std::string fxsecFilename
std::string foutputFilename
std::vector< std::string > fvndDataFilenames
std::vector< std::string > fvfdDataFilenames
std::vector< std::string > fvhelperFilenames
const NuXMLConfigfxmlConfig
Float_t fdm2NuLow
Float_t fdm2NuHigh
Float_t fdm2NuGranularity
Float_t fsn2NuLow
Float_t fsn2NuHigh
Float_t fsn2NuGranularity
Float_t fdm2BarLow
Float_t fdm2BarHigh
Float_t fdm2BarGranularity
Float_t fsn2BarLow
Float_t fsn2BarHigh
Float_t fsn2BarGranularity
Float_t ftransitionProbLow
Float_t ftransitionProbHigh
Float_t ftransitionProbGranularity
Float_t fmu2NuLow
Float_t fmu2NuHigh
Float_t fmu2NuGranularity

Detailed Description

Definition at line 217 of file NuMatrixFitter.h.


Member Enumeration Documentation

Enumerator:
kUndefined 
kOscillation 
kDecay 
kDecoherence 

Definition at line 221 of file NuMatrixFitter.h.

00221               {
00222     kUndefined = 0,
00223       kOscillation = 1,
00224       kDecay = 2,
00225       kDecoherence = 3
00226   } NuModel_t;


Constructor & Destructor Documentation

NuMatrixFitter::NuMatrixFitter (  ) 

Definition at line 318 of file NuMatrixFitter.cxx.

00319 {
00320   fndDataFilename = "";
00321   ffdDataFilename = "";
00322   fhelperFilename = "";
00323   fxsecFilename = "";
00324 
00325   SetGridParamsDefault();
00326   fxmlConfig=0;
00327 }

NuMatrixFitter::NuMatrixFitter ( const Double_t  FluxPoT,
const string  ndData,
const string  fdData,
const string  helperFile,
const string  xsecFile,
const NuXMLConfig xmlConfig = 0 
)

Definition at line 330 of file NuMatrixFitter.cxx.

References ffdDataFilename, ffluxPoT, fhelperFilename, FillGridSearchParams(), fndDataFilename, foutputFilename, fxmlConfig, fxsecFilename, and SetGridParamsDefault().

00336 {
00337   //This constructor is the setup for extrapolating single runs
00338   ffluxPoT = FluxPoT;
00339   fndDataFilename = ndData;
00340   ffdDataFilename = fdData;
00341   fhelperFilename = helperFile;
00342   fxsecFilename = xsecFile;
00343   foutputFilename =
00344     "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutput.root";
00345   fxmlConfig=xmlConfig;
00346   if (fxmlConfig) {
00347     FillGridSearchParams(fxmlConfig);
00348   } else {
00349     SetGridParamsDefault();
00350   }
00351 }

NuMatrixFitter::NuMatrixFitter ( const Double_t  fluxPoT,
const std::vector< std::string >  ndDataFiles,
const std::vector< std::string >  fdDataFiles,
const std::vector< std::string >  helperFiles,
const std::string  xsecFile,
const NuXMLConfig xmlConfig = 0 
) [explicit]
NuMatrixFitter::~NuMatrixFitter (  )  [virtual]

Definition at line 381 of file NuMatrixFitter.cxx.

00382 {
00383 }


Member Function Documentation

void NuMatrixFitter::AppearanceAnalysis (  )  [virtual]

Definition at line 2187 of file NuMatrixFitter.cxx.

References ffdDataFilename, ffluxPoT, fhelperFilename, fndDataFilename, fxsecFilename, NuMatrixMethod::GetFDOtherCCTrueUnoscBackground(), NuMatrixMethod::GetFDPotentialAppearanceFidFlux(), NuMMExtrapolation::kCCNuMuBarOscBackgroundCombined, NuMMExtrapolation::kCCNuMuOscBackgroundCombined, Munits::kilotonne, Msg::kInfo, NuBinningScheme::kUnknown, MSG, NuMatrixMethod::PredictFDFluxUnosc(), NuMatrixMethod::PredictFDSpectrumAppearance(), NuMatrixMethod::SetExtrapolationScheme(), NuMatrixMethod::SetFDAppearedFidFlux(), NuMatrixMethod::SetFDCCTrueUnoscBackground(), and Munits::tonne.

02188 {
02189   //Get the FD data & PoT histograms
02190   //Data:
02191   TFile fdDataFile(ffdDataFilename.c_str(),"READ");
02192   TH1D* numufdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
02193   numufdData->SetDirectory(0);
02194   numufdData->Sumw2();
02195   TH1D* numubarfdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
02196   numubarfdData->SetDirectory(0);
02197   numubarfdData->Sumw2();
02198   //PoT:
02199   TH1F* hfdPoTTorTGT = (TH1F*) fdDataFile.Get("hPottortgt");
02200   hfdPoTTorTGT->SetDirectory(0);
02201   TH1F* hfdSpillsPerFile = (TH1F*) fdDataFile.Get("hSpillsPerFile");
02202   hfdSpillsPerFile->SetDirectory(0);
02203   TH1F* hfdRun = (TH1F*) fdDataFile.Get("hRun");
02204   hfdRun->SetDirectory(0);
02205   fdDataFile.Close();
02206 
02207   //Get the ND PoT histograms
02208   TFile ndDataFile(fndDataFilename.c_str(),"READ");
02209   TH1F* hndPoTTorTGT = (TH1F*) ndDataFile.Get("hPottortgt");
02210   hndPoTTorTGT->SetDirectory(0);
02211   ndDataFile.Close();
02212   
02213   //Calculate the total PoT
02214   Double_t ndPoT = hndPoTTorTGT->GetMean()*hndPoTTorTGT->GetEntries()*1e12;
02215   Double_t fdPoT = hfdRun->GetEntries()*hfdPoTTorTGT->GetMean()*1e8*1e12;
02216   //   Double_t fdPoT = 6.5e20;
02217   MSG("NuMatrixFitter",Msg::kInfo)  << "ndPoT: " << ndPoT
02218        << "; fdPoT: " << fdPoT
02219        << endl;
02220   
02221   //Hard coded, unchecked fiducial masses. Should make configurable.
02222   Double_t ndFidMass = 45.128*Munits::tonne;
02223   Double_t fdFidMass = 5.487*Munits::kilotonne*14.0*TMath::Pi()*448/(52.92*486);
02224   
02225   NuMatrixMethod mmNuMu(NuBinningScheme::kUnknown,
02226                         ndPoT,
02227                         fdPoT,
02228                         ffluxPoT,
02229                         ndFidMass,
02230                         fdFidMass,
02231                         fhelperFilename,
02232                         fxsecFilename,
02233                         fndDataFilename,
02234                         "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/AppearanceOutputNuMu.root");
02235   mmNuMu.SetExtrapolationScheme
02236     (NuMMExtrapolation::kCCNuMuOscBackgroundCombined);
02237   
02238   NuMatrixMethod mmNuMuBar(NuBinningScheme::kUnknown,
02239                            ndPoT,
02240                            fdPoT,
02241                            ffluxPoT,
02242                            ndFidMass,
02243                            fdFidMass,
02244                            fhelperFilename,
02245                            fxsecFilename,
02246                            fndDataFilename,
02247                            "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/AppearanceOutputNuMuBar.root");
02248   mmNuMuBar.SetExtrapolationScheme
02249     (NuMMExtrapolation::kCCNuMuBarOscBackgroundCombined);
02250   
02251   mmNuMu.PredictFDFluxUnosc();
02252   mmNuMuBar.PredictFDFluxUnosc();
02253   
02254   mmNuMu.SetFDCCTrueUnoscBackground(mmNuMuBar.GetFDOtherCCTrueUnoscBackground());
02255   mmNuMuBar.SetFDCCTrueUnoscBackground(mmNuMu.GetFDOtherCCTrueUnoscBackground());
02256   Double_t dm2 = 2.5e-3;
02257   Double_t sn2 = 1.0;
02258   Double_t dm2bar = 0.0;//2.5e-3;
02259   Double_t sn2bar = 0.0;//1.0;
02260   mmNuMu.SetFDAppearedFidFlux(mmNuMuBar.GetFDPotentialAppearanceFidFlux());
02261   mmNuMuBar.SetFDAppearedFidFlux(mmNuMu.GetFDPotentialAppearanceFidFlux());
02262   
02263   
02264   //  const TH1D* numufdPred =
02265   mmNuMu.PredictFDSpectrumAppearance(dm2,
02266                                      sn2,
02267                                      dm2bar,
02268                                      sn2bar,
02269                                      0.2);
02270   //  const TH1D* numubarfdPred =
02271   mmNuMuBar.PredictFDSpectrumAppearance(dm2bar,
02272                                         sn2bar,
02273                                         dm2,
02274                                         sn2,
02275                                         0.2);
02276   
02277   
02278   TFile *sfile = new TFile("$SRT_PRIVATE_CONTEXT/NtupleUtils/results/AppearanceOutputNuMu.root","RECREATE");
02279   numufdData->Write();
02280   sfile->Close();
02281   if (sfile){delete sfile; sfile = 0;}
02282   
02283   TFile *sfilebar = new TFile("$SRT_PRIVATE_CONTEXT/NtupleUtils/results/AppearanceOutputNuMuBar.root","RECREATE");
02284   numubarfdData->Write();
02285   sfilebar->Close();
02286   if (sfilebar){delete sfilebar; sfilebar = 0;}
02287 }

Double_t NuMatrixFitter::CalculateChi2 ( const TH1D *  fdPred,
const TH1D *  fdData 
) const [virtual]

Definition at line 449 of file NuMatrixFitter.cxx.

00451 {
00452   //Return Chi2.
00453 
00454   Double_t chi2 = 0;
00455 
00456   for (Int_t i=1; i<=fdPred->GetNbinsX(); ++i){
00457     //Bizarre limits because root histograms are silly
00458     Double_t mnu = fdPred->GetBinContent(i);
00459     Double_t dnu = fdData->GetBinContent(i);
00460 //     MSG("NuMatrixFitter",Msg::kInfo)  << "Bin " << i
00461 //       << ": MC " << mnu
00462 //       << "; data " << dnu
00463 //       << endl;
00464     if (mnu<0.0001){mnu = 0.0001;}
00465     chi2 += (dnu-mnu)*(dnu-mnu)/mnu;
00466   }
00467   return chi2;
00468 }

Double_t NuMatrixFitter::CalculateLikelihood ( const TH1D *  fdPred,
const TH1D *  fdData,
const Double_t  nubarEcut 
) const [virtual]

Definition at line 528 of file NuMatrixFitter.cxx.

00531 {
00532   //Aim to minimise -2lnL. That's what I'm returning.
00533 
00534   Double_t like = 0;
00535   Int_t i = 1;
00536   if (i>0.0){
00537     i = fdPred->GetXaxis()->FindBin(Ecut);
00538   }
00539   for (; i<=fdPred->GetNbinsX(); ++i){
00540     //Bizarre limits because root histograms are silly
00541     Double_t mnu = fdPred->GetBinContent(i);
00542     Double_t dnu = fdData->GetBinContent(i);
00543     if (mnu<0.0001){mnu = 0.0001;}
00544     if (dnu){
00545       like += 2*(mnu - dnu + dnu*log(dnu/mnu));
00546     }
00547     else{like += 2*(mnu - dnu);}
00548   }
00549   return like;
00550 }

Double_t NuMatrixFitter::CalculateLikelihood ( const TH1D *  fdPred,
const TH1D *  fdData 
) const [virtual]

Definition at line 471 of file NuMatrixFitter.cxx.

Referenced by CCAnalysis(), CPTAnalysis(), DoCPTFit(), DoMultiRunCPTFit(), DoMultiRunPRLCCFit(), DoNoChargeCutFit(), DoPRLCCFit(), and DoTransitionFit().

00473 {
00474   //Aim to minimise -2lnL. That's what I'm returning.
00475 
00476   Double_t like = 0;
00477 
00478   for (Int_t i=1; i<=fdPred->GetNbinsX(); ++i){
00479     //Bizarre limits because root histograms are silly
00480     Double_t mnu = fdPred->GetBinContent(i);
00481     Double_t dnu = fdData->GetBinContent(i);
00482     if (mnu<0.0001){mnu = 0.0001;}
00483     if (dnu){
00484       like += 2*(mnu - dnu + dnu*log(dnu/mnu));
00485     }
00486     else{like += 2*(mnu - dnu);}
00487   }
00488 //   MSG("NuMatrixFitter",Msg::kInfo)  << "Like = " << like << endl;
00489   return like;
00490 }

Double_t NuMatrixFitter::CalculateLikelihoodNormPenalty ( const TH1D *  fdPred,
const TH1D *  fdData,
const Double_t  normalisation 
) const [virtual]

Definition at line 493 of file NuMatrixFitter.cxx.

00496 {
00497   //Aim to minimise -2lnL. That's what I'm returning.
00498 
00499   Double_t normOneSigma = 0.04;
00500   Double_t like = 0;
00501 
00502   for (Int_t i=1; i<=fdPred->GetNbinsX(); ++i){
00503     //Bizarre limits because root histograms are silly
00504     Double_t mnu = fdPred->GetBinContent(i);
00505     Double_t dnu = fdData->GetBinContent(i);
00506     if (mnu<0.0001){mnu = 0.0001;}
00507     if (dnu){
00508       like += 2*(mnu - dnu + dnu*log(dnu/mnu)) +
00509         ((normalisation - 1.0)*(normalisation-1.0))/
00510         (normOneSigma*normOneSigma);
00511     }
00512     else{like += 2*(mnu - dnu);}
00513   }
00514 //   MSG("NuMatrixFitter",Msg::kInfo)  << "Like = " << like << endl;
00515   return like;
00516 }

void NuMatrixFitter::CCAnalysis (  )  [virtual]

Definition at line 553 of file NuMatrixFitter.cxx.

References CalculateLikelihood(), fdm2NuGranularity, fdm2NuHigh, fdm2NuLow, ffdDataFilename, ffluxPoT, fhelperFilename, fndDataFilename, foutputFilename, fsn2NuGranularity, fsn2NuHigh, fsn2NuLow, fxsecFilename, NuMMExtrapolation::kCCNuMuNoOscBackground, Munits::kilotonne, Msg::kInfo, NuBinningScheme::kUnknown, MSG, NuMatrixMethod::PredictFDFluxUnosc(), NuMatrixMethod::PredictFDSpectrumBackgroundNoOsc(), NuMatrixMethod::SetExtrapolationScheme(), and Munits::tonne.

00554 {
00555   //Get the FD data & PoT histograms
00556   TFile fdDataFile(ffdDataFilename.c_str(),"READ");
00557   TH1D* fdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
00558   //  TH1D* fdData = (TH1D*) fdDataFile.Get("RecoEnergyAll_FD");
00559   fdData->SetDirectory(0);
00560   fdData->Sumw2();
00561   TH1F* hfdPoTTorTGT = (TH1F*) fdDataFile.Get("hPottortgt");
00562   hfdPoTTorTGT->SetDirectory(0);
00563   TH1F* hfdSpillsPerFile = (TH1F*) fdDataFile.Get("hSpillsPerFile");
00564   hfdSpillsPerFile->SetDirectory(0);
00565   TH1F* hfdRun = (TH1F*) fdDataFile.Get("hRun");
00566   hfdRun->SetDirectory(0);
00567   fdDataFile.Close();
00568 
00569   //Get the ND PoT histograms
00570   TFile ndDataFile(fndDataFilename.c_str(),"READ");
00571   TH1F* hndPoTTorTGT = (TH1F*) ndDataFile.Get("hPottortgt");
00572   hndPoTTorTGT->SetDirectory(0);
00573   ndDataFile.Close();
00574 
00575   //Calculate the total PoT
00576   Double_t ndPoT = hndPoTTorTGT->GetMean()*hndPoTTorTGT->GetEntries()*1e12;
00577    Double_t fdPoT = hfdRun->GetEntries()*hfdPoTTorTGT->GetMean()*1e8*1e12;
00578    //   Double_t fdPoT = 2.5e20;
00579   MSG("NuMatrixFitter",Msg::kInfo)  << "ndPoT: " << ndPoT
00580        << "; fdPoT: " << fdPoT
00581        << endl;
00582 
00583   //Hard coded, unchecked fiducial masses. Should make configurable.
00584   Double_t ndFidMass = 45.128*Munits::tonne;
00585   //  Double_t fdFidMass = ndFidMass*(448/68)*(14.0-0.3*0.3/1.0);
00586   Double_t fdFidMass = 5.487*Munits::kilotonne*14.0*TMath::Pi()*448/(52.92*486);
00587   // Double_t fdFidMass = 5.404*Munits::kilotonne*14.0*TMath::Pi()*448/(52.92*486);
00588 
00589   NuMatrixMethod mmNuMu(NuBinningScheme::kUnknown,
00590                         ndPoT,
00591                         fdPoT,
00592                         ffluxPoT,
00593                         ndFidMass,
00594                         fdFidMass,
00595                         fhelperFilename,
00596                         fxsecFilename,
00597                         fndDataFilename,
00598                         foutputFilename);
00599   mmNuMu.SetExtrapolationScheme(NuMMExtrapolation::kCCNuMuNoOscBackground);
00600   //  mmNuMu.SetExtrapolationScheme(NuMMExtrapolation::kCCNoChargeCut);
00601   mmNuMu.PredictFDFluxUnosc();
00602   
00603   Double_t bestChi2 = 1.0e10;
00604   Double_t bestsn2 = 0.0;
00605   Double_t bestdm2 = 0.0;
00606   Double_t sn2 = 0.0;
00607   Double_t dm2 = 0.0;
00608 
00609   Int_t numSn2bins = (Int_t) round((fsn2NuHigh-fsn2NuLow)/fsn2NuGranularity);
00610   Int_t numDm2bins = (Int_t) round((fdm2NuHigh-fdm2NuLow)/fdm2NuGranularity);
00611   TH2D hChi2("hChi2",
00612              "",
00613              numSn2bins,
00614              fsn2NuLow,
00615              fsn2NuHigh,
00616              numDm2bins,
00617              fdm2NuLow,
00618              fdm2NuHigh);
00619    for (/*Double_t*/ sn2 = fsn2NuLow; sn2 < fsn2NuHigh+fsn2NuGranularity; sn2 += fsn2NuGranularity){
00620      for (/*Double_t*/ dm2 = fdm2NuLow; dm2 < fdm2NuHigh+fdm2NuGranularity; dm2 += fdm2NuGranularity){
00621       const TH1D* fdPred =
00622         mmNuMu.PredictFDSpectrumBackgroundNoOsc(dm2, sn2);
00623       Double_t chi2 = this->CalculateLikelihood(fdPred,fdData);
00624       Int_t xBin = hChi2.GetXaxis()->FindBin(sn2+fsn2NuGranularity/2.0);
00625       Int_t yBin = hChi2.GetYaxis()->FindBin(dm2+fdm2NuGranularity/2.0);
00626       Int_t bin2D = hChi2.GetBin(xBin,yBin);
00627       hChi2.SetBinContent(bin2D,chi2);
00628        MSG("NuMatrixFitter",Msg::kInfo)  << "Result: " << hChi2.GetBinContent(bin2D) << endl;
00629        MSG("NuMatrixFitter",Msg::kInfo)  << "dm2: " << dm2
00630            << "; sn2: " << sn2
00631            << "; chi2: " << chi2 << endl;
00632       if (chi2 < bestChi2){
00633         bestChi2 = chi2;
00634         bestsn2 = sn2;
00635         bestdm2 = dm2;
00636       }//if
00637      }//dm2
00638    }//sn2
00639   MSG("NuMatrixFitter",Msg::kInfo)  << "Best fit: dm2 = " << bestdm2
00640        << "; sn2 = " << bestsn2
00641        << "; bestChi2 = " << bestChi2 << endl;
00642   
00643   TFile *sfile = new TFile("$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutput.root","RECREATE");
00644   fdData->Write();
00645   hChi2.Write();
00646   sfile->Close();
00647   if (sfile){delete sfile; sfile = 0;}
00648   
00649 }

virtual void NuMatrixFitter::CCFitChargeCut ( const std::string  inputFilename,
const string  fdDataFilename,
const std::string  outputFilename,
const Double_t  nubarEcut = -1.0,
const NuModel_t  model = kOscillation 
) [virtual]
void NuMatrixFitter::CPTAnalysis (  )  [virtual]

Definition at line 2290 of file NuMatrixFitter.cxx.

References CalculateLikelihood(), MuELoss::e, ffdDataFilename, ffluxPoT, fhelperFilename, fndDataFilename, fxsecFilename, NuMatrixMethod::GetFDOtherCCTrueUnoscBackground(), NuMMExtrapolation::kCCNuMuBarOscBackgroundCombined, NuMMExtrapolation::kCCNuMuOscBackgroundCombined, Munits::kilotonne, Msg::kInfo, NuBinningScheme::kUnknown, MSG, NuMatrixMethod::PredictFDFluxUnosc(), NuMatrixMethod::PredictFDSpectrumBackgroundOscCombined(), NuMatrixMethod::SetExtrapolationScheme(), NuMatrixMethod::SetFDCCTrueUnoscBackground(), and Munits::tonne.

02291 {
02292   //Get the FD data & PoT histograms
02293   //Data:
02294   TFile fdDataFile(ffdDataFilename.c_str(),"READ");
02295   TH1D* numufdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
02296   numufdData->SetDirectory(0);
02297   numufdData->Sumw2();
02298   TH1D* numubarfdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
02299   numubarfdData->SetDirectory(0);
02300   numubarfdData->Sumw2();
02301   //PoT:
02302   TH1F* hfdPoTTorTGT = (TH1F*) fdDataFile.Get("hPottortgt");
02303   hfdPoTTorTGT->SetDirectory(0);
02304   TH1F* hfdSpillsPerFile = (TH1F*) fdDataFile.Get("hSpillsPerFile");
02305   hfdSpillsPerFile->SetDirectory(0);
02306   TH1F* hfdRun = (TH1F*) fdDataFile.Get("hRun");
02307   hfdRun->SetDirectory(0);
02308   fdDataFile.Close();
02309 
02310   //Get the ND PoT histograms
02311   TFile ndDataFile(fndDataFilename.c_str(),"READ");
02312   TH1F* hndPoTTorTGT = (TH1F*) ndDataFile.Get("hPottortgt");
02313   hndPoTTorTGT->SetDirectory(0);
02314   ndDataFile.Close();
02315 
02316   //Calculate the total PoT
02317   Double_t ndPoT = hndPoTTorTGT->GetMean()*hndPoTTorTGT->GetEntries()*1e12;
02318    Double_t fdPoT = hfdRun->GetEntries()*hfdPoTTorTGT->GetMean()*1e8*1e12;
02319    //   Double_t fdPoT = 6.5e20;
02320   MSG("NuMatrixFitter",Msg::kInfo)  << "ndPoT: " << ndPoT
02321        << "; fdPoT: " << fdPoT
02322        << endl;
02323 
02324   //Hard coded, unchecked fiducial masses. Should make configurable.
02325   Double_t ndFidMass = 45.128*Munits::tonne;
02326   Double_t fdFidMass = 5.487*Munits::kilotonne*14.0*TMath::Pi()*448/(52.92*486);
02327 
02328   NuMatrixMethod mmNuMu(NuBinningScheme::kUnknown,
02329                         ndPoT,
02330                         fdPoT,
02331                         ffluxPoT,
02332                         ndFidMass,
02333                         fdFidMass,
02334                         fhelperFilename,
02335                         fxsecFilename,
02336                         fndDataFilename,
02337                         "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutputNuMu.root");
02338   mmNuMu.SetExtrapolationScheme
02339     (NuMMExtrapolation::kCCNuMuOscBackgroundCombined);
02340 
02341   NuMatrixMethod mmNuMuBar(NuBinningScheme::kUnknown,
02342                            ndPoT,
02343                            fdPoT,
02344                            ffluxPoT,
02345                            ndFidMass,
02346                            fdFidMass,
02347                            fhelperFilename,
02348                            fxsecFilename,
02349                            fndDataFilename,
02350                            "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutputNuMuBar.root");
02351   mmNuMuBar.SetExtrapolationScheme
02352     (NuMMExtrapolation::kCCNuMuBarOscBackgroundCombined);
02353 
02354   mmNuMu.PredictFDFluxUnosc();
02355   mmNuMuBar.PredictFDFluxUnosc();
02356 
02357   mmNuMu.SetFDCCTrueUnoscBackground(mmNuMuBar.GetFDOtherCCTrueUnoscBackground());
02358   mmNuMuBar.SetFDCCTrueUnoscBackground(mmNuMu.GetFDOtherCCTrueUnoscBackground());
02359 
02360   // TODO: When we start using this code again, use the member
02361   // variables for the grid search parameters
02362   Double_t bestChi2 = 1.0e10;
02363   Double_t bestsn2 = 0.0;
02364   Double_t bestdm2 = 0.0;
02365   Double_t bestsn2bar = 0.0;
02366   Double_t bestdm2bar = 0.0;
02367   Double_t sn2 = 0;
02368   Double_t dm2 = 0;
02369   Double_t sn2bar = 0;
02370   Double_t dm2bar = 0;
02371   Double_t sn2Low = 0.8;
02372   Double_t sn2High = 1.0;
02373   Double_t sn2Granularity = 0.005;
02374   Double_t dm2Low = 1.8e-3;
02375   Double_t dm2High = 3.2e-3;
02376   Double_t dm2Granularity = 0.05e-3;
02377   Double_t sn2BarLow = 0.0;
02378   Double_t sn2BarHigh = 1.0;
02379   Double_t sn2BarGranularity = 0.01;
02380   Double_t dm2BarLow = 0e-3;
02381   Double_t dm2BarHigh = 10e-3;
02382   Double_t dm2BarGranularity = 0.1e-3;
02383   /*
02384   for (sn2 = sn2Low; sn2 <= sn2High; sn2 += sn2Granularity){
02385     for (dm2 = dm2Low; dm2 <= dm2High; dm2 += dm2Granularity){
02386       for (sn2bar = sn2BarLow; sn2bar <= sn2BarHigh; sn2bar += sn2BarGranularity){
02387         for (dm2bar = dm2BarLow; dm2bar <= dm2BarHigh; dm2bar += dm2BarGranularity){
02388           const TH1D* numufdPred =
02389             mmNuMu.PredictFDSpectrumBackgroundOscCombined(dm2,
02390                                                           sn2,
02391                                                           dm2bar,
02392                                                           sn2bar);
02393           const TH1D* numubarfdPred =
02394             mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(dm2bar,
02395                                                              sn2bar,
02396                                                              dm2,
02397                                                              sn2);
02398           Double_t chi2 = this->CalculateLikelihood(numufdPred,numufdData);
02399           chi2 += this->CalculateLikelihood(numubarfdPred,numubarfdData);
02400           MSG("NuMatrixFitter",Msg::kInfo)  << "dm2 " << dm2
02401                << "; sn2 " << sn2
02402                << "; dm2bar " << dm2bar
02403                << "; sn2bar " << sn2bar
02404                << "; chi2 " << chi2
02405                << endl;
02406           if (chi2 < bestChi2){
02407             bestChi2 = chi2;
02408             bestsn2 = sn2;
02409             bestdm2 = dm2;
02410             bestsn2bar = sn2bar;
02411             bestdm2bar = dm2bar;
02412           }//if
02413         }//dm2bar
02414       }//sn2bar
02415     }//dm2
02416   }//sn2
02417   */
02418  
02419   MSG("NuMatrixFitter",Msg::kInfo)  << "Best fit: dm2 = " << bestdm2
02420        << "; sn2 = " << bestsn2
02421        << "; dm2bar = " << bestdm2bar
02422        << "; sn2bar = " << bestsn2bar
02423        << "; bestChi2 = " << bestChi2 << endl;
02424 
02425   Int_t numSn2bins = (Int_t) round((sn2High-sn2Low)/sn2Granularity) + 1;
02426   Int_t numDm2bins = (Int_t) round((dm2High-dm2Low)/dm2Granularity) + 1;
02427   //Loop again to fill chi2 surfaces.
02428   TH2D hChi2("hChi2",
02429              "",
02430              numSn2bins,
02431              sn2Low-sn2Granularity/2.0,
02432              sn2High+sn2Granularity/2.0,
02433              numDm2bins,
02434              dm2Low-dm2Granularity/2.0,
02435              dm2High+dm2Granularity/2.0);
02436   
02437 //   sn2bar = bestsn2bar;
02438 //   dm2bar = bestdm2bar;
02439   sn2bar = 1.0;
02440   dm2bar = 6.0e-3;
02441   
02442   for (sn2 = sn2Low;
02443        sn2 < sn2High+sn2Granularity/2.0;
02444        sn2 += sn2Granularity){
02445     for (dm2 = dm2Low;
02446          dm2 < dm2High+dm2Granularity/2.0;
02447          dm2 += dm2Granularity){
02448       const TH1D* numufdPred =
02449         mmNuMu.PredictFDSpectrumBackgroundOscCombined(dm2,
02450                                                       sn2,
02451                                                       dm2bar,
02452                                                       sn2bar);
02453       const TH1D* numubarfdPred =
02454         mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(dm2bar,
02455                                                          sn2bar,
02456                                                          dm2,
02457                                                          sn2);
02458       Double_t chi2 = this->CalculateLikelihood(numufdPred,numufdData);
02459       chi2 += this->CalculateLikelihood(numubarfdPred,numubarfdData);
02460       Int_t xBin = hChi2.GetXaxis()->FindBin(sn2);
02461       Int_t yBin = hChi2.GetYaxis()->FindBin(dm2);
02462       Int_t bin2D = hChi2.GetBin(xBin,yBin);
02463       MSG("NuMatrixFitter",Msg::kInfo)  << "Setting bin " << sn2
02464            << " and " << dm2
02465            << " with " << chi2
02466            << endl;
02467       hChi2.SetBinContent(bin2D,chi2);
02468     }//dm2
02469   }//sn2
02470   
02471   Int_t numSn2BarBins =
02472     (Int_t) round((sn2BarHigh-sn2BarLow)/sn2BarGranularity) + 1;
02473   Int_t numDm2BarBins =
02474     (Int_t) round((dm2BarHigh-dm2BarLow)/dm2BarGranularity) + 1;
02475   //Loop again to fill chi2bar surfaces.
02476   
02477   TH2D hChi2Bar("hChi2Bar",
02478              "",
02479              numSn2BarBins,
02480              sn2BarLow-sn2BarGranularity/2.0,
02481              sn2BarHigh+sn2BarGranularity/2.0,
02482              numDm2BarBins,
02483              dm2BarLow-dm2BarGranularity/2.0,
02484              dm2BarHigh+dm2BarGranularity/2.0);
02485   
02486 //   sn2 = bestsn2;
02487 //   dm2 = bestdm2;
02488   sn2 = 1.0;
02489   dm2 = 2.3e-3;
02490   for (sn2bar = sn2BarLow;
02491        sn2bar < sn2BarHigh+sn2BarGranularity/2.0;
02492        sn2bar += sn2BarGranularity){
02493     for (dm2bar = dm2BarLow;
02494          dm2bar < dm2BarHigh+dm2BarGranularity/2.0;
02495          dm2bar += dm2BarGranularity){
02496       const TH1D* numufdPred =
02497         mmNuMu.PredictFDSpectrumBackgroundOscCombined(dm2,
02498                                                       sn2,
02499                                                       dm2bar,
02500                                                       sn2bar);
02501       const TH1D* numubarfdPred =
02502         mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(dm2bar,
02503                                                          sn2bar,
02504                                                          dm2,
02505                                                          sn2);
02506       Double_t chi2 = this->CalculateLikelihood(numufdPred,numufdData);
02507       chi2 += this->CalculateLikelihood(numubarfdPred,numubarfdData);
02508       Int_t xBin = hChi2Bar.GetXaxis()->FindBin(sn2bar);
02509       Int_t yBin = hChi2Bar.GetYaxis()->FindBin(dm2bar);
02510       MSG("NuMatrixFitter",Msg::kInfo)  << "Setting bin " << sn2bar
02511            << " and " << dm2bar
02512            << " with " << chi2
02513            << endl;
02514       Int_t bin2D = hChi2Bar.GetBin(xBin,yBin);
02515       hChi2Bar.SetBinContent(bin2D,chi2);
02516     }//dm2Bar
02517   }//sn2Bar
02518   
02519   
02520   TFile *sfile = new TFile("$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutputNuMu.root","RECREATE");
02521   numufdData->Write();
02522   hChi2.Write();
02523   hChi2Bar.Write();
02524   sfile->Close();
02525   if (sfile){delete sfile; sfile = 0;}
02526 
02527   TFile *sfilebar = new TFile("$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutputNuMuBar.root","RECREATE");
02528   numubarfdData->Write();
02529   sfilebar->Close();
02530   if (sfilebar){delete sfilebar; sfilebar = 0;}
02531 }

virtual void NuMatrixFitter::CPTFit ( const std::string  inputFilename,
const string  fdDataFilename,
const Double_t  dm2,
const Double_t  sn2,
const std::string  outputFilename 
) [virtual]
NuMatrixOutput * NuMatrixFitter::DoCCFitChargeCut ( const NuMatrixInput mmInput,
const TH1D *  numufdData,
const TH1D *  numubarfdData,
const Double_t  nubarEcut = -1.0,
const NuModel_t  model = kOscillation 
) [virtual]

Definition at line 777 of file NuMatrixFitter.cxx.

References NuMatrixOutput::BestFitPoint(), NCType::kDecay, NCType::kDecoherence, Msg::kError, Msg::kInfo, NuMMExtrapolation::kModularNuMuBarCPTFit, NuMMExtrapolation::kModularNuMuCPTFit, MSG, NuMatrixOutput::NuMuBarBestFitFDPrediction(), NuMatrixOutput::NuMuBarChi2Surface(), NuMatrixOutput::NuMuBarFDData(), NuMatrixOutput::NuMuBarNoOscPrediction(), NuMatrixOutput::NuMuBestFitFDPrediction(), NuMatrixOutput::NuMuChi2Surface(), NuMatrixOutput::NuMuFDData(), NuMatrixOutput::NuMuNoOscPrediction(), NuMatrixMethod::PredictFDDecaySpectrumBackgroundDecayCombined(), NuMatrixMethod::PredictFDDecoherenceSpectrumBackgroundDecohereCombined(), and NuMatrixMethod::PredictFDSpectrumBackgroundOscCombined().

00782 {
00783   NuMatrixMethod mmNuMu(mmInput,
00784                         NuMMExtrapolation::kModularNuMuCPTFit);
00785   NuMatrixMethod mmNuMuBar(mmInput,
00786                            NuMMExtrapolation::kModularNuMuBarCPTFit);
00787   
00788   Double_t bestChi2 = 1.0e10;
00789   Double_t bestsn2 = 0.0;
00790   Double_t bestdm2 = 0.0;
00791   
00792   //Create the Chi2 histogram
00793   Int_t numSn2Bins =
00794     (Int_t) round((fsn2NuHigh-fsn2NuLow)/fsn2NuGranularity) + 1;
00795   Int_t numDm2Bins =
00796     (Int_t) round((fdm2NuHigh-fdm2NuLow)/fdm2NuGranularity) + 1;
00797   TH2D hChi2("hChi2",
00798              "",
00799              numSn2Bins,
00800              fsn2NuLow-fsn2NuGranularity/2.0,
00801              fsn2NuHigh+fsn2NuGranularity/2.0,
00802              numDm2Bins,
00803              fdm2NuLow-fdm2NuGranularity/2.0,
00804              fdm2NuHigh+fdm2NuGranularity/2.0); 
00805   
00806   for (Double_t sn2 = fsn2NuLow;
00807        sn2 <= fsn2NuHigh+fsn2NuGranularity/2.0;
00808        sn2 += fsn2NuGranularity){
00809     MSG("NuMatrixFitter.cxx",Msg::kInfo)
00810       << "sn2 = " << sn2 << "; best chi2 " << bestChi2 << endl;
00811     for (Double_t dm2 = fdm2NuLow;
00812          dm2 <= fdm2NuHigh+fdm2NuGranularity/2.0;
00813          dm2 += fdm2NuGranularity){
00814       const TH1D* numufdPred = 0;
00815       const TH1D* numubarfdPred = 0;
00816       if (kOscillation == model){
00817         numufdPred =
00818           mmNuMu.PredictFDSpectrumBackgroundOscCombined(dm2,
00819                                                         sn2,
00820                                                         dm2,
00821                                                         sn2);
00822         numubarfdPred =
00823           mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(dm2,
00824                                                            sn2,
00825                                                            dm2,
00826                                                            sn2);
00827       }
00828       else if (kDecay == model){
00829         numufdPred =
00830           mmNuMu.PredictFDDecaySpectrumBackgroundDecayCombined(dm2,
00831                                                                sn2,
00832                                                                dm2,
00833                                                                sn2);
00834         numubarfdPred =
00835           mmNuMuBar.PredictFDDecaySpectrumBackgroundDecayCombined(dm2,
00836                                                                   sn2,
00837                                                                   dm2,
00838                                                                   sn2);
00839 //      MSG("NuMatrixFitter",Msg::kInfo)  << "Decayed; dm2 " << dm2
00840 //           << ", sn2 " << sn2
00841 //           << ", bin 4 contents: " << numufdPred->GetBinContent(4)
00842 //           << endl;
00843       }
00844       else if (kDecoherence == model){
00845         numufdPred =
00846           mmNuMu.PredictFDDecoherenceSpectrumBackgroundDecohereCombined(dm2,
00847                                                                         sn2,
00848                                                                         dm2,
00849                                                                         sn2);
00850         numubarfdPred =
00851           mmNuMuBar.PredictFDDecoherenceSpectrumBackgroundDecohereCombined(dm2,
00852                                                                            sn2,
00853                                                                            dm2,
00854                                                                            sn2);
00855       }
00856       else {
00857         MSG("NuMatrixFitter.cxx",Msg::kError)
00858           << "Unknown neutrino disappearance model!" << endl;
00859       }
00860       
00861       Double_t chi2 = this->CalculateLikelihood(numubarfdPred,numubarfdData,nubarEcut);
00862       chi2 += this->CalculateLikelihood(numufdPred,numufdData);
00863       Int_t xBin = hChi2.GetXaxis()->FindBin(sn2);
00864       Int_t yBin = hChi2.GetYaxis()->FindBin(dm2);
00865       Int_t bin2D = hChi2.GetBin(xBin,yBin);
00866       hChi2.SetBinContent(bin2D,chi2);
00867       if (chi2 < bestChi2){
00868         bestChi2 = chi2;
00869         bestsn2 = sn2;
00870         bestdm2 = dm2;
00871       }//if
00872     }//dm2
00873   }//sn2
00874   
00875   MSG("NuMatrixFitter.cxx",Msg::kInfo)
00876     << "Best fit:" << endl
00877     << "    sin2 = " << bestsn2 << endl
00878     << "    dm2 = " << bestdm2 << endl;
00879   
00880   NuMatrixOutput *mmOutput = new NuMatrixOutput();
00881 
00882   mmOutput->NuMuBarChi2Surface(hChi2);
00883   mmOutput->NuMuBarFDData(numubarfdData);
00884   if (kOscillation == model){  
00885     mmOutput->NuMuBarBestFitFDPrediction
00886       (*mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(bestdm2,
00887                                                          bestsn2,
00888                                                          bestdm2,
00889                                                          bestsn2));
00890     mmOutput->NuMuBarNoOscPrediction
00891       (*mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(0.0,
00892                                                          0.0,
00893                                                          0.0,
00894                                                          0.0));
00895     mmOutput->NuMuBestFitFDPrediction
00896       (*mmNuMu.PredictFDSpectrumBackgroundOscCombined(bestdm2,
00897                                                       bestsn2,
00898                                                       bestdm2,
00899                                                       bestsn2));
00900     mmOutput->NuMuNoOscPrediction
00901       (*mmNuMu.PredictFDSpectrumBackgroundOscCombined(0.0,
00902                                                       0.0,
00903                                                       0.0,
00904                                                       0.0));
00905   }
00906   else if (kDecay == model){
00907     mmOutput->NuMuBarBestFitFDPrediction
00908       (*mmNuMuBar.PredictFDDecaySpectrumBackgroundDecayCombined(bestdm2,
00909                                                                 bestsn2,
00910                                                                 bestdm2,
00911                                                                 bestsn2));
00912     mmOutput->NuMuBarNoOscPrediction
00913       (*mmNuMuBar.PredictFDDecaySpectrumBackgroundDecayCombined(0.0,
00914                                                                 0.0,
00915                                                                 0.0,
00916                                                                 0.0));
00917     mmOutput->NuMuBestFitFDPrediction
00918       (*mmNuMu.PredictFDDecaySpectrumBackgroundDecayCombined(bestdm2,
00919                                                              bestsn2,
00920                                                              bestdm2,
00921                                                              bestsn2));
00922     mmOutput->NuMuNoOscPrediction
00923       (*mmNuMu.PredictFDDecaySpectrumBackgroundDecayCombined(0.0,
00924                                                              0.0,
00925                                                              0.0,
00926                                                              0.0));
00927   }
00928   else if (kDecoherence == model){
00929     mmOutput->NuMuBarBestFitFDPrediction
00930       (*mmNuMuBar.
00931        PredictFDDecoherenceSpectrumBackgroundDecohereCombined(bestdm2,
00932                                                               bestsn2,
00933                                                               bestdm2,
00934                                                               bestsn2));
00935     mmOutput->NuMuBarNoOscPrediction
00936       (*mmNuMuBar.
00937        PredictFDDecoherenceSpectrumBackgroundDecohereCombined(0.0,
00938                                                               0.0,
00939                                                               0.0,
00940                                                               0.0));
00941     mmOutput->NuMuBestFitFDPrediction
00942       (*mmNuMu.
00943        PredictFDDecoherenceSpectrumBackgroundDecohereCombined(bestdm2,
00944                                                               bestsn2,
00945                                                               bestdm2,
00946                                                               bestsn2));
00947     mmOutput->NuMuNoOscPrediction
00948       (*mmNuMu.
00949        PredictFDDecoherenceSpectrumBackgroundDecohereCombined(0.0,
00950                                                               0.0,
00951                                                               0.0,
00952                                                               0.0));
00953   }
00954   else{
00955     MSG("NuMatrixFitter.cxx",Msg::kError)
00956       << "Unknown neutrino disappearance model!" << endl;
00957   }
00958   mmOutput->BestFitPoint(bestsn2, bestdm2);
00959 
00960   mmOutput->NuMuChi2Surface(hChi2);
00961   mmOutput->NuMuFDData(numufdData);  
00962   mmOutput->BestFitPoint(bestsn2, bestdm2);
00963   
00964   return mmOutput;
00965 }

NuMatrixOutput * NuMatrixFitter::DoCPTFit ( NuMatrixInput mmInput,
TH1D *  numubarfdData,
const Double_t  dm2,
const Double_t  sn2 
) [virtual]

Definition at line 707 of file NuMatrixFitter.cxx.

References NuMatrixOutput::BestFitPoint(), CalculateLikelihood(), fdm2BarGranularity, fdm2BarHigh, fdm2BarLow, fsn2BarGranularity, fsn2BarHigh, fsn2BarLow, Msg::kInfo, MSG, NuMatrixOutput::NuMuBarBestFitFDPrediction(), NuMatrixOutput::NuMuBarChi2Surface(), NuMatrixOutput::NuMuBarFDData(), NuMatrixOutput::NuMuBarNoOscPrediction(), and NuMatrixInput::PredictFDSpectrumNuMuBar().

00711 {
00712   
00713   Double_t bestChi2 = 1.0e10;
00714   Double_t bestsn2bar = 0.0;
00715   Double_t bestdm2bar = 0.0;
00716   
00717   //Create the Chi2 histogram
00718   Int_t numSn2BarBins = (Int_t) round((fsn2BarHigh-fsn2BarLow)/fsn2BarGranularity) + 1;
00719   Int_t numDm2BarBins = (Int_t) round((fdm2BarHigh-fdm2BarLow)/fdm2BarGranularity) + 1;
00720   TH2D hChi2Bar("hChi2Bar",
00721                 "",
00722                 numSn2BarBins,
00723                 fsn2BarLow-fsn2BarGranularity/2.0,
00724                 fsn2BarHigh+fsn2BarGranularity/2.0,
00725                 numDm2BarBins,
00726                 fdm2BarLow-fdm2BarGranularity/2.0,
00727                 fdm2BarHigh+fdm2BarGranularity/2.0); 
00728   
00729   for (Double_t sn2bar = fsn2BarLow;
00730        sn2bar <= fsn2BarHigh+fsn2BarGranularity/2.0;
00731        sn2bar += fsn2BarGranularity){
00732     MSG("NuMatrixFitter.cxx",Msg::kInfo)
00733       << "sn2bar = " << sn2bar << "; bestChi2 " << bestChi2 << endl;
00734     for (Double_t dm2bar = fdm2BarLow;
00735          dm2bar <= fdm2BarHigh+fdm2BarGranularity/2.0;
00736          dm2bar += fdm2BarGranularity){
00737       
00738       TH1D* numubarfdPred = mmInput->PredictFDSpectrumNuMuBar
00739         (dm2, sn2, dm2bar, sn2bar);
00740       
00741       Double_t chi2 = this->CalculateLikelihood(numubarfdPred,numubarfdData);
00742       
00743       //       chi2 += this->CalculateLikelihood(numufdPred,numufdData);
00744       Int_t xBin = hChi2Bar.GetXaxis()->FindBin(sn2bar);
00745       Int_t yBin = hChi2Bar.GetYaxis()->FindBin(dm2bar);
00746       Int_t bin2D = hChi2Bar.GetBin(xBin,yBin);
00747       hChi2Bar.SetBinContent(bin2D,chi2);
00748       if (chi2 < bestChi2){
00749         bestChi2 = chi2;
00750         bestsn2bar = sn2bar;
00751         bestdm2bar = dm2bar;
00752       }//if
00753     }//dm2bar
00754   }//sn2bar
00755   
00756   MSG("NuMatrixFitter.cxx",Msg::kInfo)
00757     << "Best fit:" << endl
00758     << "    sin2bar = " << bestsn2bar << endl
00759     << "    dm2bar = " << bestdm2bar << endl;
00760   
00761   NuMatrixOutput *mmOutput = new NuMatrixOutput();
00762   mmOutput->NuMuBarChi2Surface(hChi2Bar);
00763   mmOutput->NuMuBarFDData(numubarfdData);
00764   
00765   mmOutput->NuMuBarBestFitFDPrediction
00766     (*mmInput->PredictFDSpectrumNuMuBar(dm2, sn2, bestdm2bar, bestsn2bar));
00767   mmOutput->NuMuBarNoOscPrediction
00768     (*mmInput->PredictFDSpectrumNuMuBar(0., 0., 0., 0.));
00769   
00770   mmOutput->BestFitPoint(bestsn2bar, bestdm2bar);
00771   
00772         return mmOutput;
00773 }

const pair< vector< NuMatrixOutput * >, NuMatrixOutput * > NuMatrixFitter::DoMultiRunCPTFit ( const vector< NuMatrixInput * >  mmInput,
const vector< TH1D >  numubarfdData,
const Double_t  dm2,
const Double_t  sn2 
) const [virtual]

Definition at line 1502 of file NuMatrixFitter.cxx.

References NuMatrixOutput::BestFitPoint(), CalculateLikelihood(), fdm2BarGranularity, fdm2BarHigh, fdm2BarLow, fsn2BarGranularity, fsn2BarHigh, fsn2BarLow, Msg::kInfo, NuMMExtrapolation::kModularNuMuBarCPTFit, MSG, NuMatrixOutput::NuMuBarBestFitFDPrediction(), NuMatrixOutput::NuMuBarChi2Surface(), NuMatrixOutput::NuMuBarFDData(), and NuMatrixOutput::NuMuBarNoOscPrediction().

01506 {
01507   vector<NuMatrixMethod*> nubarExtrapolators;
01508   for (vector<NuMatrixInput*>::const_iterator inputIt = mmInput.begin();
01509        inputIt != mmInput.end();
01510        ++inputIt){
01511     nubarExtrapolators.push_back
01512       (new NuMatrixMethod(**inputIt,
01513                           NuMMExtrapolation::kModularNuMuBarCPTFit));
01514   }
01515   
01516   Double_t bestChi2 = 1.0e10;
01517   Double_t bestsn2bar = 0.0;
01518   Double_t bestdm2bar = 0.0;
01519   
01520   //Create the Chi2 histogram
01521   Int_t numSn2BarBins =
01522     (Int_t) round((fsn2BarHigh-fsn2BarLow)/fsn2BarGranularity) + 1;
01523   Int_t numDm2BarBins =
01524     (Int_t) round((fdm2BarHigh-fdm2BarLow)/fdm2BarGranularity) + 1;
01525   TH2D hChi2Bar("hChi2Bar",
01526                 "",
01527                 numSn2BarBins,
01528                 fsn2BarLow-fsn2BarGranularity/2.0,
01529                 fsn2BarHigh+fsn2BarGranularity/2.0,
01530                 numDm2BarBins,
01531                 fdm2BarLow-fdm2BarGranularity/2.0,
01532                 fdm2BarHigh+fdm2BarGranularity/2.0); 
01533   
01534   for (Double_t sn2bar = fsn2BarLow;
01535        sn2bar <= fsn2BarHigh+fsn2BarGranularity/2.0;
01536        sn2bar += fsn2BarGranularity){
01537     MSG("NuMatrixFitter.cxx",Msg::kInfo)
01538       << "sn2bar = " << sn2bar << "; bestChi2 " << endl;
01539     for (Double_t dm2bar = fdm2BarLow;
01540          dm2bar <= fdm2BarHigh+fdm2BarGranularity/2.0;
01541          dm2bar += fdm2BarGranularity){
01542 
01543       Double_t chi2 = 0.0;
01544       vector<TH1D>::const_iterator fdDataIt = numubarfdData.begin();
01545       vector<NuMatrixMethod*>::iterator mmNuBarIt = nubarExtrapolators.begin();
01546       for (;
01547            mmNuBarIt != nubarExtrapolators.end();
01548            ++mmNuBarIt, ++fdDataIt){
01549         const TH1D* numubarfdPred =
01550           (*mmNuBarIt)->PredictFDSpectrumBackgroundOscCombined(dm2bar,
01551                                                                sn2bar,
01552                                                                dm2,
01553                                                                sn2);
01554         chi2 += this->CalculateLikelihood(numubarfdPred,&(*fdDataIt));
01555       }
01556       
01557       Int_t xBin = hChi2Bar.GetXaxis()->FindBin(sn2bar);
01558       Int_t yBin = hChi2Bar.GetYaxis()->FindBin(dm2bar);
01559       Int_t bin2D = hChi2Bar.GetBin(xBin,yBin);
01560       hChi2Bar.SetBinContent(bin2D,chi2);
01561       if (chi2 < bestChi2){
01562         bestChi2 = chi2;
01563         bestsn2bar = sn2bar;
01564         bestdm2bar = dm2bar;
01565       }//if
01566     }//dm2bar
01567   }//sn2bar
01568 
01569   MSG("NuMatrixFitter.cxx",Msg::kInfo)
01570                 << "Best fit:" << endl
01571                 << "    sin2bar = " << bestsn2bar << endl
01572     << "    dm2bar = " << bestdm2bar << endl;
01573   
01574   TH1D* hbestFitPred = 0;
01575   TH1D* hnoOscPred = 0;
01576   TH1D* hcombinedFDData = 0;
01577   vector<NuMatrixOutput*> mmOuts;
01578   Bool_t firstTime = true;
01579   vector<TH1D>::const_iterator fdDataIt = numubarfdData.begin();
01580   vector<NuMatrixMethod*>::iterator mmNuBarIt = nubarExtrapolators.begin();
01581   for (;
01582        mmNuBarIt != nubarExtrapolators.end();
01583        ++mmNuBarIt, ++fdDataIt){
01584     NuMatrixOutput* mmOutput = new NuMatrixOutput();
01585     mmOutput->NuMuBarChi2Surface(hChi2Bar);
01586     mmOutput->NuMuBarFDData(*fdDataIt);
01587     mmOutput->NuMuBarBestFitFDPrediction(*(*mmNuBarIt)->
01588                                          PredictFDSpectrumBackgroundOscCombined(bestdm2bar,
01589                                                                                 bestsn2bar,
01590                                                                                 dm2,
01591                                                                                 sn2));
01592     mmOutput->NuMuBarNoOscPrediction(*(*mmNuBarIt)->
01593                                      PredictFDSpectrumBackgroundOscCombined(0.0,
01594                                                                             0.0,
01595                                                                             0.0,
01596                                                                             0.0));
01597     mmOutput->BestFitPoint(bestsn2bar, bestdm2bar);
01598     mmOuts.push_back(mmOutput);
01599     if (firstTime){
01600       hbestFitPred = new TH1D
01601         (*(*mmNuBarIt)->
01602          PredictFDSpectrumBackgroundOscCombined(bestdm2bar,
01603                                                 bestsn2bar,
01604                                                 dm2,
01605                                                 sn2));
01606       hnoOscPred = new TH1D
01607         (*(*mmNuBarIt)->
01608          PredictFDSpectrumBackgroundOscCombined(0.0,
01609                                                 0.0,
01610                                                 0.0,
01611                                                 0.0));
01612       hcombinedFDData = new TH1D(*fdDataIt);
01613       firstTime = false;
01614     }
01615     else{
01616       hbestFitPred->Add
01617         ((*mmNuBarIt)->
01618          PredictFDSpectrumBackgroundOscCombined(bestdm2bar,
01619                                                 bestsn2bar,
01620                                                 dm2,
01621                                                 sn2));
01622       hnoOscPred->Add
01623         ((*mmNuBarIt)->
01624          PredictFDSpectrumBackgroundOscCombined(0.0,
01625                                                 0.0,
01626                                                 0.0,
01627                                                 0.0));
01628       hcombinedFDData->Add(&(*fdDataIt));
01629     }
01630   }
01631         
01632   NuMatrixOutput *mmOutputCombined = new NuMatrixOutput();
01633   mmOutputCombined->NuMuBarChi2Surface(hChi2Bar);
01634   mmOutputCombined->NuMuBarFDData(*hcombinedFDData);
01635   mmOutputCombined->NuMuBarBestFitFDPrediction(*hbestFitPred);
01636   mmOutputCombined->NuMuBarNoOscPrediction(*hnoOscPred);
01637   mmOutputCombined->BestFitPoint(bestsn2bar, bestdm2bar);
01638 
01639   pair<vector<NuMatrixOutput*>,NuMatrixOutput*> allMMOuts(mmOuts,mmOutputCombined);
01640   return allMMOuts;
01641 }

const pair< vector< NuMatrixOutput * >, NuMatrixOutput * > NuMatrixFitter::DoMultiRunPRLCCFit ( const vector< NuMatrixInput * >  mmInput,
const vector< TH1D >  numufdData,
const NuModel_t  model = kOscillation 
) [virtual]

Definition at line 1674 of file NuMatrixFitter.cxx.

References NuMatrixOutput::BestFitPoint(), CalculateLikelihood(), fdm2NuGranularity, fdm2NuHigh, fdm2NuLow, fsn2NuGranularity, fsn2NuHigh, fsn2NuLow, kDecay, kDecoherence, Msg::kError, Msg::kInfo, NuMMExtrapolation::kModularPRLCCFit, kOscillation, MSG, NormalisationPenaltyTerm(), NuMatrixOutput::NuMuBestFitFDPrediction(), NuMatrixOutput::NuMuChi2Surface(), NuMatrixOutput::NuMuFDData(), NuMatrixOutput::NuMuNoOscPrediction(), and RebinForFit().

01677 {
01678   vector<NuMatrixMethod*> numuExtrapolators;
01679   for (vector<NuMatrixInput*>::const_iterator inputIt = mmInput.begin();
01680        inputIt != mmInput.end();
01681        ++inputIt){
01682     numuExtrapolators.push_back
01683       (new NuMatrixMethod(**inputIt,
01684                           NuMMExtrapolation::kModularPRLCCFit));
01685   }
01686   
01687   Double_t bestChi2 = 1.0e10;
01688   Double_t bestsn2 = 0.0;
01689   Double_t bestdm2 = 0.0;
01690   Double_t bestNorm = 0.0;
01691 
01692 //   fsn2NuLow = 0.0;
01693 //     fsn2NuHigh = 0.000;
01694 //     fsn2NuGranularity = 0.001;
01695 //     fdm2NuLow = 0.0e-3;
01696 //     fdm2NuHigh = 0.0e-3;
01697 //     fdm2NuGranularity = 0.005e-3;
01698     Double_t normLow =1.0;
01699     Double_t normHigh = 1.0;
01700     Double_t normGranularity = 0.0005;
01701 
01702   
01703   //Create the Chi2 histogram
01704   Int_t numSn2Bins =
01705     (Int_t) round((fsn2NuHigh-fsn2NuLow)/fsn2NuGranularity) + 1;
01706   Int_t numDm2Bins =
01707     (Int_t) round((fdm2NuHigh-fdm2NuLow)/fdm2NuGranularity) + 1;
01708   TH2D hChi2("hChi2",
01709              "",
01710              numSn2Bins,
01711              fsn2NuLow-fsn2NuGranularity/2.0,
01712              fsn2NuHigh+fsn2NuGranularity/2.0,
01713              numDm2Bins,
01714              fdm2NuLow-fdm2NuGranularity/2.0,
01715              fdm2NuHigh+fdm2NuGranularity/2.0); 
01716   
01717   for (Double_t sn2 = fsn2NuLow;
01718        sn2 <= fsn2NuHigh+fsn2NuGranularity/2.0;
01719        sn2 += fsn2NuGranularity){
01720     for (Double_t dm2 = fdm2NuLow;
01721          dm2 <= fdm2NuHigh+fdm2NuGranularity/2.0;
01722          dm2 += fdm2NuGranularity){
01723       
01724       Double_t chi2BestNorm = 1.0e10;
01725       Double_t bestNormThisLoop = 0.0;
01726 
01727       for (Double_t normalisation = normLow;
01728            normalisation <= normHigh+normGranularity/2.0;
01729            normalisation += normGranularity){
01730         
01731         Double_t chi2 = 0.0;
01732         
01733         vector<TH1D>::const_iterator fdDataIt = numufdData.begin();
01734         vector<NuMatrixMethod*>::iterator mmNuMuIt = numuExtrapolators.begin();
01735         for (;
01736              mmNuMuIt != numuExtrapolators.end();
01737              ++mmNuMuIt, ++fdDataIt){
01738           const TH1D* numufdPred = 0;
01739           if (kOscillation == model){
01740             numufdPred =
01741               (*mmNuMuIt)->PredictFDSpectrumBackgroundNoOsc(dm2,
01742                                                             sn2);
01743           }
01744           else if (kDecay == model){
01745             numufdPred =
01746               (*mmNuMuIt)->PredictFDDecaySpectrumBackgroundNoOsc(dm2,
01747                                                                  sn2);
01748           }
01749           else if (kDecoherence == model){
01750             numufdPred =
01751               (*mmNuMuIt)->PredictFDDecoherenceSpectrumBackgroundNoOsc(dm2,
01752                                                                        sn2);
01753           }
01754           else {
01755             MSG("NuMatrixFitter.cxx",Msg::kError)
01756               << "Invalid neutrino disappearance model" << endl;
01757           }
01758 
01759 
01760           TH1D numufdPredScaled = this->RebinForFit(numufdPred);
01761           numufdPredScaled.Scale(normalisation);
01762           //      TH1D* pPredScaled = numufdPredScaled.Rebin(numFitBins,"ScaledPred",fitBinsArray);
01763           TH1D rebinnedData = this->RebinForFit(&(*fdDataIt));
01764 //        TH1* pRebinnedData = rebinnedData.Rebin(numFitBins,"ScaledData",fitBinsArray);
01765           chi2 += this->CalculateLikelihood(&numufdPredScaled,&rebinnedData);
01766         }
01767         chi2 += this->NormalisationPenaltyTerm(normalisation);
01768         if (chi2 < bestChi2){
01769           bestChi2 = chi2;
01770           bestsn2 = sn2;
01771           bestdm2 = dm2;
01772           bestNorm = normalisation;
01773         }//if
01774         if (chi2 < chi2BestNorm){
01775           chi2BestNorm = chi2;
01776           bestNormThisLoop = normalisation;
01777         }//if
01778         MSG("NuMatrixFitter.cxx",Msg::kInfo)
01779           << "sn2 = " << sn2
01780           << ", dm2 = " << dm2
01781           << ", norm = " << normalisation
01782           << ", penalty = " << this->NormalisationPenaltyTerm(normalisation)
01783           << "; bestChi2 " << bestChi2 
01784           << "; bestsn2: " << bestsn2
01785           << "; bestdm2: " << bestdm2 
01786           << "; bestNorm: " << bestNorm << endl;
01787       }//normalisation
01788       Int_t xBin = hChi2.GetXaxis()->FindBin(sn2);
01789       Int_t yBin = hChi2.GetYaxis()->FindBin(dm2);
01790       Int_t bin2D = hChi2.GetBin(xBin,yBin);
01791       hChi2.SetBinContent(bin2D,chi2BestNorm);
01792     }//dm2
01793   }//sn2
01794 
01795   MSG("NuMatrixFitter.cxx",Msg::kInfo)
01796                 << "Best fit:" << endl
01797                 << "    sin2 = " << bestsn2 << endl
01798                 << "    dm2 = " << bestdm2 << endl
01799                 << "    norm = " << bestNorm << endl
01800                 << "    chi2 = " << bestChi2 << endl;
01801   
01802   TH1D* hbestFitPred = 0;
01803   TH1D* hnoOscPred = 0;
01804   TH1D* hcombinedFDData = 0;
01805   vector<NuMatrixOutput*> mmOuts;
01806   Bool_t firstTime = true;
01807   vector<TH1D>::const_iterator fdDataIt = numufdData.begin();
01808   vector<NuMatrixMethod*>::iterator mmNuMuIt = numuExtrapolators.begin();
01809   for (;
01810        mmNuMuIt != numuExtrapolators.end();
01811        ++mmNuMuIt, ++fdDataIt){
01812     NuMatrixOutput* mmOutput = new NuMatrixOutput();
01813     mmOutput->NuMuChi2Surface(hChi2);
01814     mmOutput->NuMuFDData(*fdDataIt);
01815     mmOutput->NuMuBestFitFDPrediction(*(*mmNuMuIt)->
01816                                       PredictFDSpectrumBackgroundNoOsc(bestdm2,
01817                                                                        bestsn2));
01818     mmOutput->NuMuNoOscPrediction(*(*mmNuMuIt)->
01819                                   PredictFDSpectrumBackgroundNoOsc(0.0,
01820                                                                    0.0));
01821     mmOutput->BestFitPoint(bestsn2, bestdm2);
01822     mmOuts.push_back(mmOutput);
01823     if (firstTime){
01824       if (kOscillation == model){
01825         hbestFitPred = new TH1D
01826           (*(*mmNuMuIt)->
01827            PredictFDSpectrumBackgroundNoOsc(bestdm2,
01828                                             bestsn2));
01829         hnoOscPred = new TH1D
01830           (*(*mmNuMuIt)->
01831            PredictFDSpectrumBackgroundNoOsc(0.0,
01832                                             0.0));
01833       }
01834       else if (kDecay == model){
01835         hbestFitPred = new TH1D
01836           (*(*mmNuMuIt)->
01837            PredictFDDecaySpectrumBackgroundNoOsc(bestdm2,
01838                                                  bestsn2));
01839         hnoOscPred = new TH1D
01840           (*(*mmNuMuIt)->
01841            PredictFDDecaySpectrumBackgroundNoOsc(0.0,
01842                                                  0.0));
01843       }
01844       else if (kDecoherence == model){
01845         hbestFitPred = new TH1D
01846           (*(*mmNuMuIt)->
01847            PredictFDDecoherenceSpectrumBackgroundNoOsc(bestdm2,
01848                                                  bestsn2));
01849         hnoOscPred = new TH1D
01850           (*(*mmNuMuIt)->
01851            PredictFDDecoherenceSpectrumBackgroundNoOsc(0.0,
01852                                                  0.0));
01853       }
01854       hcombinedFDData = new TH1D(*fdDataIt);
01855       firstTime = false;
01856     }
01857     else{
01858       if (kOscillation == model){
01859         hbestFitPred->Add
01860           ((*mmNuMuIt)->
01861            PredictFDSpectrumBackgroundNoOsc(bestdm2,
01862                                             bestsn2));
01863         hnoOscPred->Add
01864           ((*mmNuMuIt)->
01865            PredictFDSpectrumBackgroundNoOsc(0.0,
01866                                             0.0));
01867       }
01868       else if (kDecay == model){
01869         hbestFitPred->Add
01870           ((*mmNuMuIt)->
01871            PredictFDDecaySpectrumBackgroundNoOsc(bestdm2,
01872                                                  bestsn2));
01873         hnoOscPred->Add
01874           ((*mmNuMuIt)->
01875            PredictFDDecaySpectrumBackgroundNoOsc(0.0,
01876                                                  0.0));
01877       }
01878       else if (kDecoherence == model){
01879         hbestFitPred->Add
01880           ((*mmNuMuIt)->
01881            PredictFDDecoherenceSpectrumBackgroundNoOsc(bestdm2,
01882                                                  bestsn2));
01883         hnoOscPred->Add
01884           ((*mmNuMuIt)->
01885            PredictFDDecoherenceSpectrumBackgroundNoOsc(0.0,
01886                                                  0.0));
01887       }
01888       hcombinedFDData->Add(&(*fdDataIt));
01889     }
01890   }
01891   
01892   NuMatrixOutput *mmOutputCombined = new NuMatrixOutput();
01893   mmOutputCombined->NuMuChi2Surface(hChi2);
01894   hcombinedFDData = new TH1D(this->RebinForFit(hcombinedFDData));
01895   mmOutputCombined->NuMuFDData(*hcombinedFDData);
01896   hbestFitPred = new TH1D(this->RebinForFit(hbestFitPred));
01897   hbestFitPred->Scale(bestNorm);
01898   mmOutputCombined->NuMuBestFitFDPrediction(*hbestFitPred);
01899   hnoOscPred = new TH1D(this->RebinForFit(hnoOscPred));
01900   mmOutputCombined->NuMuNoOscPrediction(*hnoOscPred);
01901   mmOutputCombined->BestFitPoint(bestsn2, bestdm2);
01902 
01903   pair<vector<NuMatrixOutput*>,NuMatrixOutput*> allMMOuts(mmOuts,mmOutputCombined);
01904   return allMMOuts;
01905 }

const NuMatrixOutput * NuMatrixFitter::DoMultiRunTransitionFit ( const vector< NuMatrixInput mmInput,
const vector< TH1D >  numubarfdData,
const Double_t  dm2,
const Double_t  sn2 
) const [virtual]

Definition at line 1205 of file NuMatrixFitter.cxx.

References NuMatrixOutput::BestFitPoint(), Msg::kInfo, NuMMExtrapolation::kModularNuMuBarTransitionFit, MSG, NuMatrixOutput::NuMuBarBestFitFDPrediction(), NuMatrixOutput::NuMuBarChi2TransSurface(), NuMatrixOutput::NuMuBarFDData(), NuMatrixOutput::NuMuBarNoOscPrediction(), and NuMatrixOutput::NuMuBarNoTransPrediction().

01209 {
01210   vector<NuMatrixMethod*> nubarExtrapolators;
01211   for (vector<NuMatrixInput>::const_iterator inputIt = mmInput.begin();
01212        inputIt != mmInput.end();
01213        ++inputIt){
01214     nubarExtrapolators.push_back
01215       (new NuMatrixMethod(*inputIt,
01216                           NuMMExtrapolation::kModularNuMuBarTransitionFit));
01217   }
01218 
01219   Double_t bestChi2 = 1.0e10;
01220   Double_t bestTransitionProb = 0.0;
01221 
01222 //   Double_t transitionProbLow = 0.0;
01223 //   Double_t transitionProbHigh = 1.0;
01224 //   Double_t transitionProbGranularity = 0.05;
01225 
01226   //Create the chi2 histogram
01227   Int_t numTransitionProbBins =
01228     (Int_t) round((ftransitionProbHigh-ftransitionProbLow)/
01229                   ftransitionProbGranularity) + 1;
01230   TH1D hChi2Trans("hChi2Trans",
01231                   "",
01232                   numTransitionProbBins,
01233                   ftransitionProbLow-ftransitionProbGranularity/2.0,
01234                   ftransitionProbHigh+ftransitionProbGranularity/2.0);
01235   for (Double_t transitionProb = ftransitionProbLow;
01236        transitionProb <= ftransitionProbHigh+ftransitionProbGranularity/2.0;
01237        transitionProb += ftransitionProbGranularity){
01238     MSG("NuMatrixFitter.cxx",Msg::kInfo)
01239       << "Fitting transitionProb = " << transitionProb << endl;
01240 
01241       Double_t chi2 = 0.0;
01242       vector<TH1D>::const_iterator fdDataIt = numubarfdData.begin();
01243       vector<NuMatrixMethod*>::iterator mmNuBarIt = nubarExtrapolators.begin();
01244       for (;
01245            mmNuBarIt != nubarExtrapolators.end();
01246            ++mmNuBarIt, ++fdDataIt){
01247         const TH1D* numubarfdPred =
01248           (*mmNuBarIt)->PredictFDSpectrumAppearance(dm2,
01249                                                     sn2,
01250                                                     dm2,
01251                                                     sn2,
01252                                                     transitionProb);  
01253         chi2 += this->CalculateLikelihood(numubarfdPred,&(*fdDataIt));
01254       }
01255 
01256     Int_t xBin = hChi2Trans.GetXaxis()->FindBin(transitionProb);
01257     hChi2Trans.SetBinContent(xBin,chi2);
01258     if (chi2 < bestChi2){
01259       bestChi2 = chi2;
01260       bestTransitionProb = transitionProb;
01261     }
01262   }//transitionProb
01263 
01264   MSG("NuMatrixFitter.cxx",Msg::kInfo)
01265     << "Best fit: transitionProb = " << bestTransitionProb << endl;
01266   
01267   TH1D* hbestFitPred = 0;
01268   TH1D* hnoTransPred = 0;
01269   TH1D* hnoOscPred = 0;
01270   TH1D* hcombinedFDData = 0;
01271   Bool_t firstTime = true;
01272   vector<TH1D>::const_iterator fdDataIt = numubarfdData.begin();
01273   vector<NuMatrixMethod*>::iterator mmNuBarIt = nubarExtrapolators.begin();
01274   for (;
01275        mmNuBarIt != nubarExtrapolators.end();
01276        ++mmNuBarIt, ++fdDataIt){
01277     if (firstTime){
01278       hbestFitPred = new TH1D
01279         (*(*mmNuBarIt)->
01280          PredictFDSpectrumAppearance(dm2,
01281                                      sn2,
01282                                      dm2,
01283                                      sn2,
01284                                      bestTransitionProb));
01285       hnoTransPred = new TH1D
01286         (*(*mmNuBarIt)->
01287          PredictFDSpectrumAppearance(dm2,
01288                                      sn2,
01289                                      dm2,
01290                                      sn2,
01291                                      0.0));
01292       hnoOscPred = new TH1D
01293         (*(*mmNuBarIt)->
01294          PredictFDSpectrumAppearance(0.0,
01295                                      0.0,
01296                                      0.0,
01297                                      0.0,
01298                                      0.0));
01299       hcombinedFDData = new TH1D(*fdDataIt);
01300       firstTime = false;
01301     }
01302     else{
01303       hbestFitPred->Add
01304         ((*mmNuBarIt)->
01305          PredictFDSpectrumAppearance(dm2,
01306                                      sn2,
01307                                      dm2,
01308                                      sn2,
01309                                      bestTransitionProb));
01310       hnoTransPred->Add
01311         ((*mmNuBarIt)->
01312          PredictFDSpectrumAppearance(dm2,
01313                                      sn2,
01314                                      dm2,
01315                                      sn2,
01316                                      0.0));
01317       hnoOscPred->Add
01318         ((*mmNuBarIt)->
01319          PredictFDSpectrumAppearance(0.0,
01320                                      0.0,
01321                                      0.0,
01322                                      0.0,
01323                                      0.0));
01324       hcombinedFDData->Add(&(*fdDataIt));
01325     }
01326   }
01327         
01328   NuMatrixOutput *mmOutput = new NuMatrixOutput();
01329   mmOutput->NuMuBarChi2TransSurface(hChi2Trans);
01330   mmOutput->NuMuBarFDData(*hcombinedFDData);
01331   mmOutput->NuMuBarBestFitFDPrediction(*hbestFitPred);
01332   mmOutput->NuMuBarNoTransPrediction(*hnoTransPred);
01333   mmOutput->NuMuBarNoOscPrediction(*hnoOscPred);
01334   mmOutput->BestFitPoint(sn2, dm2, bestTransitionProb);
01335   return mmOutput;
01336 }

NuMatrixOutput * NuMatrixFitter::DoNoChargeCutFit ( NuMatrixInput mmInput,
TH1D *  numubarfdData,
const NuModel_t  model = kOscillation 
) [virtual]

Definition at line 968 of file NuMatrixFitter.cxx.

References NuMatrixOutput::BestFitPoint(), CalculateLikelihood(), fdm2NuGranularity, fdm2NuHigh, fdm2NuLow, fsn2NuGranularity, fsn2NuHigh, fsn2NuLow, kDecay, kDecoherence, Msg::kInfo, NuMMExtrapolation::kModularNoChargeCutFit, kOscillation, MSG, NuMatrixOutput::NuMuBarBestFitFDPrediction(), NuMatrixOutput::NuMuBarChi2Surface(), NuMatrixOutput::NuMuBarFDData(), NuMatrixOutput::NuMuBarNoOscPrediction(), NuMatrixMethod::PredictFDDecaySpectrumBackgroundNoOsc(), NuMatrixMethod::PredictFDDecoherenceSpectrumBackgroundNoOsc(), and NuMatrixMethod::PredictFDSpectrumBackgroundNoOsc().

00971 {
00972   NuMatrixMethod mmNuMuBar(mmInput,
00973                            NuMMExtrapolation::kModularNoChargeCutFit);
00974   
00975   Double_t bestChi2 = 1.0e10;
00976   Double_t bestsn2 = 0.0;
00977   Double_t bestdm2 = 0.0;
00978   
00979   //Create the Chi2 histogram
00980   Int_t numSn2Bins =
00981     (Int_t) round((fsn2NuHigh-fsn2NuLow)/fsn2NuGranularity) + 1;
00982   Int_t numDm2Bins =
00983     (Int_t) round((fdm2NuHigh-fdm2NuLow)/fdm2NuGranularity) + 1;
00984   TH2D hChi2("hChi2",
00985                 "",
00986                 numSn2Bins,
00987                 fsn2NuLow-fsn2NuGranularity/2.0,
00988                 fsn2NuHigh+fsn2NuGranularity/2.0,
00989                 numDm2Bins,
00990                 fdm2NuLow-fdm2NuGranularity/2.0,
00991                 fdm2NuHigh+fdm2NuGranularity/2.0); 
00992 
00993   for (Double_t sn2 = fsn2NuLow;
00994        sn2 <= fsn2NuHigh+fsn2NuGranularity/2.0;
00995        sn2 += fsn2NuGranularity){
00996     MSG("NuMatrixFitter.cxx",Msg::kInfo)
00997       << "sn2 = " << sn2 << "; bestChi2 " << bestChi2 
00998       << " (bestdm2 " << bestdm2
00999       << ", bestsn2 " << bestsn2 << ")" << endl;
01000     for (Double_t dm2 = fdm2NuLow;
01001                                  dm2 <= fdm2NuHigh+fdm2NuGranularity/2.0;
01002                                  dm2 += fdm2NuGranularity){
01003       const TH1D* fdPred = 0;
01004       if (kOscillation == model){
01005         fdPred =
01006           mmNuMuBar.PredictFDSpectrumBackgroundNoOsc(dm2,
01007                                                      sn2);
01008       }
01009       else if (kDecay == model){
01010         fdPred =
01011           mmNuMuBar.PredictFDDecaySpectrumBackgroundNoOsc(dm2,
01012                                                           sn2);
01013       }
01014       else if (kDecoherence == model){
01015         fdPred =
01016           mmNuMuBar.PredictFDDecoherenceSpectrumBackgroundNoOsc(dm2,
01017                                                                 sn2);
01018       }
01019       else {
01020         MSG("NuMatrixFitter.cxx",Msg::kInfo)
01021           << "Invalid neutrino disappearance model" << endl;
01022       }
01023       
01024       Double_t chi2 = this->CalculateLikelihood(fdPred,fdData);
01025       //       chi2 += this->CalculateLikelihood(numufdPred,numufdData);
01026       Int_t xBin = hChi2.GetXaxis()->FindBin(sn2);
01027       Int_t yBin = hChi2.GetYaxis()->FindBin(dm2);
01028       Int_t bin2D = hChi2.GetBin(xBin,yBin);
01029       hChi2.SetBinContent(bin2D,chi2);
01030       if (chi2 < bestChi2){
01031                                 bestChi2 = chi2;
01032                                 bestsn2 = sn2;
01033                                 bestdm2 = dm2;
01034       }//if
01035     }//dm2
01036   }//sn2
01037 
01038   MSG("NuMatrixFitter.cxx",Msg::kInfo)
01039                 << "Best fit:" << endl
01040                 << "    sin2 = " << bestsn2 << endl
01041     << "    dm2 = " << bestdm2 << endl;
01042         
01043   NuMatrixOutput *mmOutput = new NuMatrixOutput();
01044   mmOutput->NuMuBarChi2Surface(hChi2);
01045   mmOutput->NuMuBarFDData(fdData);
01046   if (kOscillation == model){
01047     mmOutput->NuMuBarBestFitFDPrediction
01048       (*mmNuMuBar.PredictFDSpectrumBackgroundNoOsc(bestdm2,
01049                                                    bestsn2));
01050     
01051     mmOutput->NuMuBarNoOscPrediction
01052       (*mmNuMuBar.PredictFDSpectrumBackgroundNoOsc(0.0,
01053                                                    0.0));
01054   }
01055   else if (kDecay == model){
01056     mmOutput->NuMuBarBestFitFDPrediction
01057       (*mmNuMuBar.PredictFDDecaySpectrumBackgroundNoOsc(bestdm2,
01058                                                         bestsn2));
01059     
01060     mmOutput->NuMuBarNoOscPrediction
01061       (*mmNuMuBar.PredictFDDecaySpectrumBackgroundNoOsc(0.0,
01062                                                         0.0));
01063   }
01064   else if (kDecoherence == model){
01065     mmOutput->NuMuBarBestFitFDPrediction
01066       (*mmNuMuBar.PredictFDDecoherenceSpectrumBackgroundNoOsc(bestdm2,
01067                                                         bestsn2));
01068     
01069     mmOutput->NuMuBarNoOscPrediction
01070       (*mmNuMuBar.PredictFDDecoherenceSpectrumBackgroundNoOsc(0.0,
01071                                                               0.0));
01072   }
01073   else{
01074     MSG("NuMatrixFitter.cxx",Msg::kInfo)
01075       << "Invalid neutrino disappearance model" << endl;
01076   }
01077   mmOutput->BestFitPoint(bestsn2, bestdm2);
01078   
01079         return mmOutput;
01080 }

NuMatrixOutput * NuMatrixFitter::DoPRLCCFit ( const NuMatrixInput mmInput,
const TH1D *  fdData,
const NuModel_t  model = kOscillation 
) [virtual]

Definition at line 1083 of file NuMatrixFitter.cxx.

References NuMatrixOutput::BestFitPoint(), CalculateLikelihood(), fdm2NuGranularity, fdm2NuHigh, fdm2NuLow, fsn2NuGranularity, fsn2NuHigh, fsn2NuLow, kDecay, kDecoherence, Msg::kError, Msg::kInfo, NuMMExtrapolation::kModularPRLCCFit, kOscillation, MSG, NuMatrixOutput::NuMuBarBestFitFDPrediction(), NuMatrixOutput::NuMuBarChi2Surface(), NuMatrixOutput::NuMuBarFDData(), NuMatrixOutput::NuMuBarNoOscPrediction(), NuMatrixMethod::PredictFDDecaySpectrumBackgroundNoOsc(), NuMatrixMethod::PredictFDDecoherenceSpectrumBackgroundNoOsc(), and NuMatrixMethod::PredictFDSpectrumBackgroundNoOsc().

Referenced by PRLCCFit().

01086 {
01087   NuMatrixMethod mmNuMuBar(mmInput,
01088                            NuMMExtrapolation::kModularPRLCCFit);
01089   
01090   Double_t bestChi2 = 1.0e10;
01091   Double_t bestsn2 = 0.0;
01092   Double_t bestdm2 = 0.0;
01093   
01094   //Create the Chi2 histogram
01095   Int_t numSn2Bins =
01096     (Int_t) round((fsn2NuHigh-fsn2NuLow)/fsn2NuGranularity) + 1;
01097   Int_t numDm2Bins =
01098     (Int_t) round((fdm2NuHigh-fdm2NuLow)/fdm2NuGranularity) + 1;
01099   TH2D hChi2("hChi2",
01100                 "",
01101                 numSn2Bins,
01102                 fsn2NuLow-fsn2NuGranularity/2.0,
01103                 fsn2NuHigh+fsn2NuGranularity/2.0,
01104                 numDm2Bins,
01105                 fdm2NuLow-fdm2NuGranularity/2.0,
01106                 fdm2NuHigh+fdm2NuGranularity/2.0); 
01107 
01108   for (Double_t sn2 = fsn2NuLow;
01109        sn2 <= fsn2NuHigh+fsn2NuGranularity/2.0;
01110        sn2 += fsn2NuGranularity){
01111     MSG("NuMatrixFitter.cxx",Msg::kInfo)
01112       << "sn2 = " << sn2 << "; bestChi2 " << bestChi2 
01113       << " (bestdm2 " << bestdm2
01114       << ", bestsn2 " << bestsn2 << ")" << endl;
01115     for (Double_t dm2 = fdm2NuLow;
01116          dm2 <= fdm2NuHigh+fdm2NuGranularity/2.0;
01117          dm2 += fdm2NuGranularity){
01118       const TH1D* fdPred = 0;
01119       if (kOscillation == model){
01120         fdPred =
01121           mmNuMuBar.PredictFDSpectrumBackgroundNoOsc(dm2,
01122                                                      sn2);
01123       }
01124       else if (kDecay == model){
01125         fdPred =
01126           mmNuMuBar.PredictFDDecaySpectrumBackgroundNoOsc(dm2,
01127                                                           sn2);
01128       }
01129       else if (kDecoherence == model){
01130         fdPred =
01131           mmNuMuBar.PredictFDDecoherenceSpectrumBackgroundNoOsc(dm2,
01132                                                                 sn2);
01133       }
01134       else {
01135         MSG("NuMatrixFitter.cxx",Msg::kError)
01136           << "Invalid neutrino disappearance model" << endl;
01137       }
01138       // 
01139       //TH1D* outputhist = (TH1D*)fdPred->Clone();
01140       //outputhist->SetName(Form("fdpred_s%.2f_m%.2f", sn2, dm2*1000));
01141       //outputhist->Write();
01142       Double_t chi2 = this->CalculateLikelihood(fdPred,fdData);
01143       Int_t xBin = hChi2.GetXaxis()->FindBin(sn2);
01144       Int_t yBin = hChi2.GetYaxis()->FindBin(dm2);
01145       Int_t bin2D = hChi2.GetBin(xBin,yBin);
01146       hChi2.SetBinContent(bin2D,chi2);
01147       //MSG("NuMatrixFitter.cxx",Msg::kInfo)
01148       //  << "sn2=" << sn2 << " dm2=" << dm2 << " chi2=" << chi2 << endl;
01149       if (chi2 < bestChi2){
01150                                 bestChi2 = chi2;
01151                                 bestsn2 = sn2;
01152                                 bestdm2 = dm2;
01153       }//if
01154     }//dm2
01155   }//sn2
01156 
01157   MSG("NuMatrixFitter.cxx",Msg::kInfo)
01158                 << "Best fit:" << endl
01159                 << "    sin2 = " << bestsn2 << endl
01160     << "    dm2 = " << bestdm2 << endl;
01161         
01162   NuMatrixOutput *mmOutput = new NuMatrixOutput();
01163   mmOutput->NuMuBarChi2Surface(hChi2);
01164   mmOutput->NuMuBarFDData(fdData);
01165   
01166   if (kOscillation == model){
01167     mmOutput->NuMuBarBestFitFDPrediction
01168       (*mmNuMuBar.PredictFDSpectrumBackgroundNoOsc(bestdm2,
01169                                                    bestsn2));
01170     
01171     mmOutput->NuMuBarNoOscPrediction
01172       (*mmNuMuBar.PredictFDSpectrumBackgroundNoOsc(0.0,
01173                                                    0.0));
01174   }
01175   else if (kDecay == model){
01176     mmOutput->NuMuBarBestFitFDPrediction
01177       (*mmNuMuBar.PredictFDDecaySpectrumBackgroundNoOsc(bestdm2,
01178                                                         bestsn2));
01179     
01180     mmOutput->NuMuBarNoOscPrediction
01181       (*mmNuMuBar.PredictFDDecaySpectrumBackgroundNoOsc(0.0,
01182                                                         0.0));
01183   }
01184   else if (kDecoherence == model){
01185     mmOutput->NuMuBarBestFitFDPrediction
01186       (*mmNuMuBar.PredictFDDecoherenceSpectrumBackgroundNoOsc(bestdm2,
01187                                                         bestsn2));
01188     
01189     mmOutput->NuMuBarNoOscPrediction
01190       (*mmNuMuBar.PredictFDDecoherenceSpectrumBackgroundNoOsc(0.0,
01191                                                               0.0));
01192   }
01193   else{
01194     MSG("NuMatrixFitter.cxx",Msg::kInfo)
01195       << "Invalid neutrino disappearance model" << endl;
01196   }
01197 
01198   mmOutput->BestFitPoint(bestsn2, bestdm2);
01199   MSG("NuMatrixFitter.cxx",Msg::kInfo) << "returning mmoutput" << endl;
01200   return mmOutput;
01201 }

NuMatrixOutput * NuMatrixFitter::DoTransitionFit ( const char *  inputFilename,
const TH1D *  numubarfdData,
const Double_t  dm2,
const Double_t  sn2 
) const [virtual]

Definition at line 652 of file NuMatrixFitter.cxx.

References NuMatrixOutput::BestFitPoint(), CalculateLikelihood(), Msg::kInfo, MSG, NuMatrixOutput::NuMuBarBestFitFDPrediction(), NuMatrixOutput::NuMuBarFDData(), NuMatrixOutput::NuMuBarNoOscPrediction(), NuMatrixOutput::NuMuBarNoTransPrediction(), NuMatrixOutput::NuMuChi2TransSurface(), and NuMatrixInput::PredictFDSpectrumTransition().

Referenced by TransitionFit().

00657 {
00658     static int calls = 0;
00659     TString hN = "";
00660     hN += calls;
00661     ++calls;
00662     NuMatrixInput *mmInput = new NuMatrixInput(inputFilename);
00663     
00664     Double_t bestChi2 = 1.0e10;
00665     Double_t bestTransitionProb = 0.0;
00666     
00667     float transitionProbLow = -0.5;
00668     float transitionProbHigh = 0.5;
00669     float transitionProbGranularity = 0.05;
00670 
00671     TH1D* numubarfdPred = 0;
00672     int bins = (int)((transitionProbHigh - transitionProbLow)/transitionProbGranularity+1);
00673     TH1D* surface = new TH1D("hLikelihood"+hN,"Likelihood Curve", bins, transitionProbLow, transitionProbHigh);
00674     int ibin = 1;
00675     for (Double_t transitionProb = transitionProbLow;
00676          transitionProb <= transitionProbHigh+transitionProbGranularity/2.0;
00677          transitionProb += transitionProbGranularity){
00678         numubarfdPred = mmInput->PredictFDSpectrumTransition(dm2, sn2, transitionProb);
00679         Double_t chi2 = this->CalculateLikelihood(numubarfdPred,numubarfdData);
00680         surface->SetBinContent(ibin++, chi2);
00681 
00682         
00683         if (chi2 < bestChi2){
00684             bestChi2 = chi2;
00685             bestTransitionProb = transitionProb;
00686         }
00687         delete numubarfdPred; numubarfdPred = 0;
00688     }//transitionProb
00689     
00690     MSG("NuMatrixFitter.cxx",Msg::kInfo) << "Best fit: transitionProb = " << bestTransitionProb << endl;
00691     
00692     NuMatrixOutput *mmOutput = new NuMatrixOutput();
00693     mmOutput->NuMuBarFDData(numubarfdData);
00694     
00695     mmOutput->NuMuBarBestFitFDPrediction(*mmInput->PredictFDSpectrumTransition(dm2, sn2, bestTransitionProb));
00696     mmOutput->NuMuBarNoTransPrediction(*mmInput->PredictFDSpectrumTransition(dm2, sn2, 0.0)); 
00697     mmOutput->NuMuBarNoOscPrediction(*mmInput->PredictFDSpectrumTransition(0.0, 0.0, 0.0));   
00698     mmOutput->NuMuChi2TransSurface(*surface);
00699     mmOutput->BestFitPoint(sn2, dm2, bestTransitionProb);
00700     
00701     delete mmInput;
00702     delete surface;
00703     return mmOutput;
00704 }

void NuMatrixFitter::FillGridSearchParams ( const NuXMLConfig xmlConfig  )  [private]

Definition at line 420 of file NuMatrixFitter.cxx.

References NuXMLConfig::DM2BarGranularity(), NuXMLConfig::DM2BarHigh(), NuXMLConfig::DM2BarLow(), NuXMLConfig::DM2NuGranularity(), NuXMLConfig::DM2NuHigh(), NuXMLConfig::DM2NuLow(), fdm2BarGranularity, fdm2BarHigh, fdm2BarLow, fdm2NuGranularity, fdm2NuHigh, fdm2NuLow, fmu2NuGranularity, fmu2NuHigh, fmu2NuLow, fsn2BarGranularity, fsn2BarHigh, fsn2BarLow, fsn2NuGranularity, fsn2NuHigh, fsn2NuLow, ftransitionProbGranularity, ftransitionProbHigh, ftransitionProbLow, NuXMLConfig::MU2NuGranularity(), NuXMLConfig::MU2NuHigh(), NuXMLConfig::MU2NuLow(), NuXMLConfig::SN2BarGranularity(), NuXMLConfig::SN2BarHigh(), NuXMLConfig::SN2BarLow(), NuXMLConfig::SN2NuGranularity(), NuXMLConfig::SN2NuHigh(), NuXMLConfig::SN2NuLow(), NuXMLConfig::TransitionProbGranularity(), NuXMLConfig::TransitionProbHigh(), and NuXMLConfig::TransitionProbLow().

Referenced by NuMatrixFitter(), PRLCCFit(), SetXMLConfig(), and TransitionFit().

00421 {
00422    fdm2NuLow = xmlConfig->DM2NuLow();
00423    fdm2NuHigh = xmlConfig->DM2NuHigh();
00424    fdm2NuGranularity = xmlConfig->DM2NuGranularity();
00425 
00426    fsn2NuLow = xmlConfig->SN2NuLow();
00427    fsn2NuHigh = xmlConfig->SN2NuHigh();
00428    //MSG("NuMatrixFitter",Msg::kInfo)  << "SN2NuGranularity = " << xmlConfig->SN2NuGranularity() << endl;
00429    fsn2NuGranularity = xmlConfig->SN2NuGranularity();
00430 
00431    fdm2BarLow = xmlConfig->DM2BarLow();
00432    fdm2BarHigh = xmlConfig->DM2BarHigh();
00433    fdm2BarGranularity = xmlConfig->DM2BarGranularity();
00434 
00435    fsn2BarLow = xmlConfig->SN2BarLow();
00436    fsn2BarHigh = xmlConfig->SN2BarHigh();
00437    fsn2BarGranularity = xmlConfig->SN2BarGranularity();
00438 
00439    ftransitionProbLow = xmlConfig->TransitionProbLow();
00440    ftransitionProbHigh = xmlConfig->TransitionProbHigh();
00441    ftransitionProbGranularity = xmlConfig->TransitionProbGranularity();
00442 
00443    fmu2NuLow = xmlConfig->MU2NuLow();
00444    fmu2NuHigh = xmlConfig->MU2NuHigh();
00445    fmu2NuGranularity = xmlConfig->MU2NuGranularity();
00446 }

virtual void NuMatrixFitter::MultiRunCPTFit ( const std::vector< std::string >  inputFilenames,
const std::vector< std::string >  fdDataFilenames,
const Double_t  dm2,
const Double_t  sn2,
const std::vector< std::string >  outputFilenames,
const std::string  combinedOutputFilename 
) [virtual]
virtual void NuMatrixFitter::MultiRunPRLCCFit ( const vector< string >  inputFilenames,
const vector< string >  fdDataFilenames,
const std::vector< std::string >  outputFilenames,
const std::string  combinedOutputFilename,
const NuModel_t  model = kOscillation 
) [virtual]
virtual void NuMatrixFitter::NoChargeCutFit ( const std::string  inputFilename,
const string  fdDataFilename,
const std::string  outputFilename,
const NuModel_t  model = kOscillation 
) [virtual]
Double_t NuMatrixFitter::NormalisationPenaltyTerm ( const Double_t  normalisation  )  const

Definition at line 519 of file NuMatrixFitter.cxx.

Referenced by DoMultiRunPRLCCFit().

00520 {
00521   Double_t normOneSigma = 0.04;
00522   return ((normalisation - 1.0)*(normalisation-1.0))/
00523     (normOneSigma*normOneSigma);
00524 }

void NuMatrixFitter::NuMuBarAppearanceAnalysis (  )  [virtual]

Definition at line 3050 of file NuMatrixFitter.cxx.

03051 {
03052 
03053 }

void NuMatrixFitter::PRLCCFit ( const string  inputFilename,
const string  fdDataFilename,
const string  outputFilename,
const NuModel_t  model = kOscillation 
)

Definition at line 2066 of file NuMatrixFitter.cxx.

References DoPRLCCFit(), FillGridSearchParams(), NuXMLConfig::FullTitle(), fxmlConfig, Msg::kInfo, MSG, and NuMatrixOutput::XMLConfig().

02070 {  
02071   MSG("NuMatrixFitter",Msg::kInfo)  << "NuMatrixFitter::PRLCCFit()" << endl;
02072   TFile fdDataFile(fdDataFilename.c_str(),"READ");
02073   TH1D* fdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
02074   fdData->SetDirectory(0);
02075   fdData->Sumw2();
02076   fdDataFile.Close();
02077   
02078   TFile *infile = new TFile(inputFilename.c_str(), "READ");
02079   
02080   NuMatrixInput mmInput(infile);
02081 
02082   TString objName = "NuMatrixOutput";
02083   if (!fxmlConfig) {
02084     fxmlConfig = (NuXMLConfig*)gROOT->FindObject("NuXMLConfig");
02085   }
02086   
02087   MSG("NuMatrixFitter",Msg::kInfo)  << "fxmlConfig=" << fxmlConfig << endl;
02088   if (fxmlConfig) {
02089     FillGridSearchParams(fxmlConfig);
02090     MSG("NuMatrixFitter",Msg::kInfo)  << "This file fit is for systematic: " << fxmlConfig->FullTitle() << endl;
02091 
02092     //objName = fxmlConfig->FullName();
02093   }
02094   else {
02095     MSG("NuMatrixFitter",Msg::kInfo)  << "No NuXMLConfig object -- nominal fit assumed." << endl;
02096     //objName = "Nominal";
02097   }
02098 
02099   NuMatrixOutput *mmOutput = this->DoPRLCCFit(mmInput, fdData, model);
02100   if (fxmlConfig) {
02101      mmOutput->XMLConfig(fxmlConfig);
02102   }
02103   MSG("NuMatrixFitter",Msg::kInfo)  << "Recreating file " << outputFilename << endl;
02104   TFile *outputFile = new TFile(outputFilename.c_str(),"RECREATE");
02105   outputFile->cd();
02106   mmOutput->Write(objName);
02107   outputFile->Close();
02108   infile->cd();
02109   MSG("NuMatrixFitter",Msg::kInfo)  << "Don't mind the segfault.  Everything is done and it appears to be harmless." << endl;
02110   //    infile->Close();
02111 }

TH1D NuMatrixFitter::RebinForFit ( const TH1D *  fineHist  ) 

Definition at line 1645 of file NuMatrixFitter.cxx.

References Form().

Referenced by DoMultiRunPRLCCFit().

01646 {
01647 //   NuBinningScheme::NuBinningScheme_t binningScheme = NuBinningScheme::k0_5GeVTo200;
01648 //   const NuUtilities utils;
01649 //   const vector<Double_t> recoBins = utils.RecoBins(binningScheme);
01650 //   const Int_t numFitBins = recoBins.size()-1;
01651 //   Double_t *fitBinsArray;
01652 //   fitBinsArray=new Double_t[numFitBins+1];
01653 //   {
01654 //     Int_t i=0;
01655 //     for (vector<Double_t>::const_iterator itBin = recoBins.begin();
01656 //       itBin != recoBins.end();
01657 //       ++itBin, ++i){
01658 //       fitBinsArray[i] = *itBin;
01659 //     }
01660 //   }
01661 
01662   static Int_t stupidRoot = 0;
01663   TH1D newHist(Form("NewHist%d",stupidRoot++),"",800,0,200);
01664   newHist.SetBinContent(1,0);
01665   for (Int_t oldBin = 1; oldBin<fineHist->GetNbinsX()+1; ++oldBin){
01666     newHist.SetBinContent(oldBin+1,fineHist->GetBinContent(oldBin));
01667   }
01668   newHist.Rebin(2);
01669   return newHist;
01670 }

void NuMatrixFitter::SetGridParamsDefault (  )  [private]

Definition at line 393 of file NuMatrixFitter.cxx.

References fdm2BarGranularity, fdm2BarHigh, fdm2BarLow, fdm2NuGranularity, fdm2NuHigh, fdm2NuLow, fmu2NuGranularity, fmu2NuHigh, fmu2NuLow, fsn2BarGranularity, fsn2BarHigh, fsn2BarLow, fsn2NuGranularity, fsn2NuHigh, fsn2NuLow, ftransitionProbGranularity, ftransitionProbHigh, and ftransitionProbLow.

Referenced by NuMatrixFitter().

00394 {
00395   fdm2NuLow=0.001;
00396   fdm2NuHigh=0.01;
00397   fdm2NuGranularity=0.001;
00398 
00399   fsn2NuLow=0;
00400   fsn2NuHigh=1;
00401   fsn2NuGranularity=0.05;
00402 
00403   fdm2BarLow=0.001;
00404   fdm2BarHigh=0.01;
00405   fdm2BarGranularity=0.001;
00406 
00407   fsn2BarLow=0;
00408   fsn2BarHigh=1;
00409   fsn2BarGranularity=0.05;
00410 
00411   ftransitionProbLow=0;
00412   ftransitionProbHigh=1;
00413   ftransitionProbGranularity=0.05;
00414 
00415   fmu2NuLow=0.001;
00416   fmu2NuHigh=0.01;
00417   fmu2NuGranularity=0.001; 
00418 }

void NuMatrixFitter::SetXMLConfig ( const NuXMLConfig xmlConfig  ) 

Definition at line 386 of file NuMatrixFitter.cxx.

References FillGridSearchParams(), and fxmlConfig.

00387 {
00388   fxmlConfig=xmlConfig;
00389   this->FillGridSearchParams(xmlConfig);
00390 }

void NuMatrixFitter::TransitionFit ( const string  inputFilename,
const string  fdDataFilename,
const string  outputFilename,
const Int_t  experiments 
)

Definition at line 2114 of file NuMatrixFitter.cxx.

References NuXMLConfig::DM2Nu(), DoTransitionFit(), FillGridSearchParams(), NuXMLConfig::FullTitle(), fxmlConfig, Msg::kInfo, MSG, NuXMLConfig::SN2Nu(), and NuMatrixOutput::XMLConfig().

02118 {  
02119     TFile fdDataFile(fdDataFilename.c_str(),"READ");
02120     TH1D* fdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
02121     fdData->SetDirectory(0);
02122     fdData->Sumw2();
02123     fdDataFile.Close();
02124     
02125     TFile *infile = new TFile(inputFilename.c_str());
02126     TString objName = "NuMatrixOutput";
02127     if (!fxmlConfig) {
02128         fxmlConfig = (NuXMLConfig*)gROOT->FindObject("NuXMLConfig");
02129     }
02130     infile->Close();
02131     
02132     
02133     MSG("NuMatrixFitter",Msg::kInfo) << "fxmlConfig=" << fxmlConfig << endl;
02134     if (fxmlConfig) {
02135         FillGridSearchParams(fxmlConfig);
02136         MSG("NuMatrixFitter",Msg::kInfo) << "This file fit is for systematic: " << fxmlConfig->FullTitle() << endl;
02137         
02138         //objName = fxmlConfig->FullName();
02139     }
02140     else {
02141         MSG("NuMatrixFitter",Msg::kInfo) << "No NuXMLConfig object -- nominal fit assumed." << endl;
02142         //objName = "Nominal";
02143     }
02144     
02145     
02146     MSG("NuMatrixFitter",Msg::kInfo) << "Recreating file " << outputFilename << endl;
02147     TFile *outputFile = new TFile(outputFilename.c_str(),"RECREATE");
02148 
02149     NuFluctuator *fl = new NuFluctuator();
02150     TH1D *fd_sub_Data = 0;
02151     TH1F *bestFits = new TH1F("bestFits","Best Fit Points",40,0.0,1.0);
02152     NuMatrixOutput *mmOutput = DoTransitionFit(inputFilename.c_str()/*mmInput*/, 
02153                                                fdData, 
02154                                                fxmlConfig->DM2Nu(), 
02155                                                fxmlConfig->SN2Nu());
02156     
02157     if (fxmlConfig) mmOutput->XMLConfig(fxmlConfig);
02158     
02159     TString name = objName;
02160     name += "_";
02161     name += "nofluct";
02162     outputFile->WriteObject(mmOutput,name);
02163     delete fd_sub_Data;
02164     
02165     for (int i = 0; i < experiments; i++) {
02166         fd_sub_Data = new TH1D(fl->Fluctuate(fdData));
02167         NuMatrixOutput *mmOutput = DoTransitionFit(inputFilename.c_str()/*mmInput*/, 
02168                                                    fd_sub_Data, 
02169                                                    fxmlConfig->DM2Nu(), 
02170                                                    fxmlConfig->SN2Nu());
02171         
02172         if (fxmlConfig) mmOutput->XMLConfig(fxmlConfig);
02173 
02174         bestFits->Fill(mmOutput->bestTransitionProb);
02175         TString name = objName;
02176         name += "_";
02177         name += i;
02178         outputFile->WriteObject(mmOutput,name);
02179         delete fd_sub_Data;
02180     }
02181     outputFile->WriteObject(bestFits,bestFits->GetName());
02182     outputFile->Close();
02183     MSG("NuMatrixFitter",Msg::kInfo) << "Don't mind the segfault.  Everything is done and it appears to be harmless." << endl;
02184 }

void NuMatrixFitter::WriteHistosForFitter ( TString  outputFileName = "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutputForFitter.root"  ) 

Definition at line 2534 of file NuMatrixFitter.cxx.

References NuXMLConfig::BinningScheme(), ffdDataFilename, ffluxPoT, fhelperFilename, NuCuts::FiducialMass(), fndDataFilename, fxmlConfig, fxsecFilename, NuMatrixMethod::GetFDOtherCCTrueUnoscBackground(), NuMatrixMethod::GetFDPotentialAppearanceFidFlux(), NuMMExtrapolation::kCCNuMuBarOscBackgroundCombined, NuMMExtrapolation::kCCNuMuOscBackgroundCombined, Detector::kFar, Msg::kInfo, NuCuts::kJJE1, SimFlag::kMC, Detector::kNear, MSG, NuXMLConfig::OverriddenFDPoTs(), NuMatrixMethod::PredictFDFluxUnosc(), NuMatrixMethod::PredictFDSpectrumBackgroundOscCombined(), NuXMLConfig::ScaledFDPoTs(), NuMatrixMethod::SetExtrapolationScheme(), NuMatrixMethod::SetFDAppearedFidFlux(), NuMatrixMethod::SetFDCCTrueUnoscBackground(), and NuMatrixMethod::WriteFilesForFitter().

02535 {
02536 
02537   //Get the FD data & PoT histograms
02538   //Data:
02539   TFile fdDataFile(ffdDataFilename.c_str(),"READ");
02540   TH1D* numufdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
02541   numufdData->SetDirectory(0);
02542   numufdData->Sumw2();
02543   TH1D* numubarfdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
02544   numubarfdData->SetDirectory(0);
02545   numubarfdData->Sumw2();
02546   //PoT:
02547   TH1F* hfdTotalPot = (TH1F*) fdDataFile.Get("hTotalPot");
02548   hfdTotalPot->SetDirectory(0);
02549   fdDataFile.Close();
02550 
02551   //Get the ND PoT histograms
02552   TFile ndDataFile(fndDataFilename.c_str(),"READ");
02553   TH1F* hndTotalPot = (TH1F*) ndDataFile.Get("hTotalPot");
02554   hndTotalPot->SetDirectory(0);
02555   
02556   //Open the helper file.
02557   TFile helperFile(fhelperFilename.c_str(),"READ");
02558 
02559   // Get the NuXMLConfig Object. Look in the ND data file first. If
02560   // it's not in there, look in the helper file.
02561   if (!fxmlConfig) {
02562     fxmlConfig = (NuXMLConfig*) ndDataFile.Get("NuXMLConfig");
02563   }
02564   if (!fxmlConfig){
02565      fxmlConfig = (NuXMLConfig*)helperFile.Get("NuXMLConfig");
02566    }
02567 
02568   ndDataFile.Close();
02569 
02570   //Calculate the total PoT
02571   Double_t ndPoT = hndTotalPot->Integral();
02572   Double_t fdPoT = hfdTotalPot->Integral();
02573   if (fxmlConfig->OverriddenFDPoTs()>0) {
02574     fdPoT = fxmlConfig->OverriddenFDPoTs();
02575   }
02576 
02577   // If scaling is used, the scaling should have been done by this
02578   // stage, so use the scaled pots instead of the overridden pots
02579   if (fxmlConfig->ScaledFDPoTs()>0) {
02580     fdPoT = fxmlConfig->ScaledFDPoTs();
02581   }
02582 
02583   MSG("NuMatrixFitter",Msg::kInfo) << "ndPoT: " << ndPoT
02584        << "; fdPoT: " << fdPoT
02585        << endl;
02586 
02587   NuCuts nuCuts;
02588   Double_t ndFidMass = nuCuts.FiducialMass(Detector::kNear,
02589                                            SimFlag::kMC,
02590                                            NuCuts::kJJE1);
02591   Double_t fdFidMass = nuCuts.FiducialMass(Detector::kFar,
02592                                            SimFlag::kMC,
02593                                            NuCuts::kJJE1);
02594   MSG("NuMatrixFitter",Msg::kInfo)  << "Near mass: " << ndFidMass
02595        << "; far mas: " << fdFidMass << endl;
02596 
02597   //Get the binning scheme to use
02598   NuBinningScheme::NuBinningScheme_t binningScheme =
02599     static_cast<NuBinningScheme::NuBinningScheme_t>(fxmlConfig->BinningScheme());
02600   MSG("NuMatrixFitter",Msg::kInfo)  << "Using binningScheme " << binningScheme << endl;
02601 
02602   NuMatrixMethod mmNuMu(binningScheme,
02603                         ndPoT,
02604                         fdPoT,
02605                         ffluxPoT,
02606                         ndFidMass,
02607                         fdFidMass,
02608                         fhelperFilename,
02609                         fxsecFilename,
02610                         fndDataFilename,
02611                         "");
02612 
02613   //pull in all the histograms from the file
02614   mmNuMu.SetExtrapolationScheme
02615     (NuMMExtrapolation::kCCNuMuOscBackgroundCombined);
02616 
02617   //set all the parameters
02618   //pass in the filenames
02619   NuMatrixMethod mmNuMuBar(binningScheme,
02620                            ndPoT,
02621                            fdPoT,
02622                            ffluxPoT,
02623                            ndFidMass,
02624                            fdFidMass,
02625                            fhelperFilename,
02626                            fxsecFilename,
02627                            fndDataFilename,
02628                            "");
02629 
02630   //pull in all the histograms from the file
02631   mmNuMuBar.SetExtrapolationScheme
02632     (NuMMExtrapolation::kCCNuMuBarOscBackgroundCombined);
02633   cout << "Debug Start" << endl;
02634   //do lots and lots of steps
02635   mmNuMu.PredictFDFluxUnosc();
02636   mmNuMuBar.PredictFDFluxUnosc();
02637   cout << "Debug End" << endl;
02638   //use extrapolated numu spectrum to predict numu contamination in 
02639   //numubar spectrum (and vice versa)
02640   //The "K" spectrum to pass out to external fitter is now created
02641   mmNuMu.SetFDCCTrueUnoscBackground(mmNuMuBar.GetFDOtherCCTrueUnoscBackground());
02642   mmNuMuBar.SetFDCCTrueUnoscBackground(mmNuMu.GetFDOtherCCTrueUnoscBackground());
02643 
02644   //For a nubar appearance analyses: Use the exptrapolated numu
02645   //spectrum to create the spectrum of nubars that could appear. This
02646   //then needs to be oscillated (numu DISAPPEARANCE) to get the actual
02647   //spectrum that could appear. Oscillation is done in the fitting code.
02648   mmNuMu.SetFDAppearedFidFlux(mmNuMuBar.GetFDPotentialAppearanceFidFlux());
02649   mmNuMuBar.SetFDAppearedFidFlux(mmNuMu.GetFDPotentialAppearanceFidFlux());
02650 
02651   mmNuMu.PredictFDSpectrumBackgroundOscCombined(0.0025,1,
02652                                                 0.0025,1);
02653   mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(0.0025,1,
02654                                                    0.0025,1);
02655 
02656 
02657   //write out files
02658   //create the file
02659   TFile* sfile = new TFile(outputFileName,"RECREATE");
02660   MSG("NuMatrixFitter",Msg::kInfo) <<"Created file:"<<endl;
02661   sfile->Print();
02662 
02663   fxmlConfig->Write();
02664   
02665   //write files
02666   mmNuMu.WriteFilesForFitter("NM");
02667   mmNuMuBar.WriteFilesForFitter("NMB");
02668 
02669   //close file
02670   sfile->Close();
02671   if (sfile){delete sfile; sfile = 0;}  
02672 
02673 
02674 /*
02675           const TH1D* numufdPred =
02676             mmNuMu.PredictFDSpectrumBackgroundOscCombined(dm2,
02677                                                           sn2,
02678                                                           dm2bar,
02679                                                           sn2bar);
02680           const TH1D* numubarfdPred =
02681             mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(dm2bar,
02682                                                              sn2bar,
02683                                                              dm2,
02684                                                              sn2);
02685 */
02686 
02687 }

void NuMatrixFitter::WriteMultiRunHistosForFitter ( const vector< TString >  outputFileNames  )  [virtual]

Definition at line 2691 of file NuMatrixFitter.cxx.

References NuMMExtrapolation::kCCNuMuBarOscBackgroundCombined, NuMMExtrapolation::kCCNuMuOscBackgroundCombined, Munits::kilotonne, Msg::kInfo, MSG, NuMatrixMethod::SetExtrapolationScheme(), and Munits::tonne.

02692 {
02693   //Hard coded, unchecked fiducial masses. Should make configurable.
02694   Double_t ndFidMass = 45.128*Munits::tonne;
02695   Double_t fdFidMass = 5.487*Munits::kilotonne*14.0*TMath::Pi()/52.92;
02696 
02697   //FD data
02698   vector<TH1D> numufdData;
02699   vector<TH1D> nubarfdData;
02700   //NuMu extrapolators
02701   vector<NuMatrixMethod*> numuExtrapolators;
02702   vector<NuMatrixMethod*> nubarExtrapolators;
02703   //XML configs
02704   vector<NuXMLConfig> xmlConfigs;
02705   
02706   vector<string>::const_iterator helperFileIt = fvhelperFilenames.begin();
02707   vector<string>::const_iterator fdFileIt = fvfdDataFilenames.begin();
02708   for (vector<string>:: const_iterator ndFileIt = fvndDataFilenames.begin();
02709        ndFileIt != fvndDataFilenames.end();
02710        ++ndFileIt, ++fdFileIt, ++helperFileIt){
02711     //Calculate ND PoT
02712     TFile ndDataFile((*ndFileIt).c_str(),"READ");
02713     TH1F* hndPoTTorTGT = (TH1F*) ndDataFile.Get("hPottortgt");
02714     Double_t ndPoT = hndPoTTorTGT->GetMean()*hndPoTTorTGT->GetEntries()*1e12;
02715     ndDataFile.Close();
02716 
02717     //Get FD data
02718     TFile fdDataFile((*fdFileIt).c_str(),"READ");
02719     TH1D* numuHist = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
02720     numuHist->SetDirectory(0);
02721     numuHist->Sumw2();
02722     numufdData.push_back(*numuHist);
02723     TH1D* numubarHist = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
02724     numubarHist->SetDirectory(0);
02725     numubarHist->Sumw2();
02726     nubarfdData.push_back(*numubarHist);
02727 
02728     //Get FD PoT
02729     TH1F* hfdPoTTorTGT = (TH1F*) fdDataFile.Get("hPottortgt");
02730     TH1F* hfdRun = (TH1F*) fdDataFile.Get("hRun");
02731     Double_t fdPoT = hfdRun->GetEntries()*hfdPoTTorTGT->GetMean()*1e8*1e12;
02732     fdDataFile.Close();
02733 
02734     //Get NuXMLConfig objects
02735     TFile helperFile((*helperFileIt).c_str(),"READ");
02736     NuXMLConfig *xmlConfig = (NuXMLConfig*) helperFile.Get("NuXMLConfig");
02737     xmlConfigs.push_back(*xmlConfig);
02738 
02739     //Binning scheme
02740     NuBinningScheme::NuBinningScheme_t binningScheme =
02741       static_cast<NuBinningScheme::NuBinningScheme_t>(xmlConfig->BinningScheme());
02742     MSG("NuMatrixFitter",Msg::kInfo)  << "Using binning scheme" << binningScheme << endl;
02743 
02744     //Construct the extrapolators
02745     NuMatrixMethod* numuExtrap = new NuMatrixMethod(binningScheme,
02746                                                     ndPoT,
02747                                                     fdPoT,
02748                                                     ffluxPoT,
02749                                                     ndFidMass,
02750                                                     fdFidMass,
02751                                                     *helperFileIt,
02752                                                     fxsecFilename,
02753                                                     *ndFileIt,
02754                                                     "");
02755     numuExtrap->SetExtrapolationScheme
02756       (NuMMExtrapolation::kCCNuMuOscBackgroundCombined);
02757     numuExtrapolators.push_back(numuExtrap);
02758     
02759     NuMatrixMethod* nubarExtrap = new NuMatrixMethod(binningScheme,
02760                                                      ndPoT,
02761                                                      fdPoT,
02762                                                      ffluxPoT,
02763                                                      ndFidMass,
02764                                                      fdFidMass,
02765                                                      *helperFileIt,
02766                                                      fxsecFilename,
02767                                                      *ndFileIt,
02768                                                      "");
02769     nubarExtrap->SetExtrapolationScheme
02770       (NuMMExtrapolation::kCCNuMuBarOscBackgroundCombined);
02771     nubarExtrapolators.push_back(nubarExtrap);
02772   }
02773 
02774   //Now do the extrapolation steps
02775   vector<NuMatrixMethod*>::iterator numuExtrapIt = numuExtrapolators.begin();
02776   vector<NuMatrixMethod*>::iterator nubarExtrapIt = nubarExtrapolators.begin();
02777   vector<TString>::const_iterator outfileIt = outputFileNames.begin();
02778   vector<NuXMLConfig>::iterator xmlConfigIt = xmlConfigs.begin();
02779   for (;
02780        numuExtrapIt != numuExtrapolators.end();
02781        ++numuExtrapIt, ++nubarExtrapIt, ++outfileIt){
02782     (*numuExtrapIt)->PredictFDFluxUnosc();
02783     (*nubarExtrapIt)->PredictFDFluxUnosc();
02784 
02785     (*numuExtrapIt)->SetFDCCTrueUnoscBackground
02786       ((*nubarExtrapIt)->GetFDOtherCCTrueUnoscBackground());
02787     (*nubarExtrapIt)->SetFDCCTrueUnoscBackground
02788       ((*numuExtrapIt)->GetFDOtherCCTrueUnoscBackground());
02789 
02790     //For nubar appearance analysis
02791     (*numuExtrapIt)->SetFDAppearedFidFlux
02792       ((*nubarExtrapIt)->GetFDPotentialAppearanceFidFlux());
02793     (*nubarExtrapIt)->SetFDAppearedFidFlux
02794       ((*numuExtrapIt)->GetFDPotentialAppearanceFidFlux());
02795 
02796     //Don't think we need to run this stage:
02797     (*numuExtrapIt)->PredictFDSpectrumBackgroundOscCombined(0.0025,1,
02798                                                           0.0025,1);
02799     (*nubarExtrapIt)->PredictFDSpectrumBackgroundOscCombined(0.0025,1,
02800                                                            0.0025,1);
02801 
02802     //Write the output files
02803     TFile* sfile = new TFile(*outfileIt,"RECREATE");
02804     MSG("NuMatrixFitter",Msg::kInfo) <<"Created file:"<<endl;
02805     sfile->Print();
02806     (*xmlConfigIt).Write();
02807     (*numuExtrapIt)->WriteFilesForFitter("NM");
02808     (*nubarExtrapIt)->WriteFilesForFitter("NMB");
02809     sfile->Close();
02810     if (sfile){delete sfile; sfile = 0;} 
02811   }
02812 }

void NuMatrixFitter::WriteNoChargeCutHistosForFitter ( std::string  outputFileName = ""  ) 

Definition at line 2817 of file NuMatrixFitter.cxx.

References NuCuts::FiducialMass(), NuMMExtrapolation::kCCNoChargeCut, Detector::kFar, Msg::kInfo, NuCuts::kJJE1, SimFlag::kMC, Detector::kNear, MSG, NuMatrixMethod::PredictFDFluxUnosc(), NuMatrixMethod::SetExtrapolationScheme(), and NuMatrixMethod::WriteFilesForFitter().

02818 {
02819   //Get the FD data & PoT histograms
02820   //Data:
02821   TFile fdDataFile(ffdDataFilename.c_str(),"READ");
02822   TH1D* numufdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
02823   numufdData->SetDirectory(0);
02824   numufdData->Sumw2();
02825   TH1D* numubarfdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
02826   numubarfdData->SetDirectory(0);
02827   numubarfdData->Sumw2();
02828   //PoT:
02829   TH1F* hfdTotalPot = (TH1F*) fdDataFile.Get("hTotalPot");
02830   hfdTotalPot->SetDirectory(0);
02831   fdDataFile.Close();
02832 
02833   //Get the ND PoT histograms
02834   TFile ndDataFile(fndDataFilename.c_str(),"READ");
02835   TH1F* hndTotalPot = (TH1F*) ndDataFile.Get("hTotalPot");
02836   hndTotalPot->SetDirectory(0);
02837   
02838   //Open the helper file.
02839   TFile helperFile(fhelperFilename.c_str(),"READ");
02840 
02841   // Get the NuXMLConfig Object. Look in the ND data file first. If
02842   // it's not in there, look in the helper file.
02843   if (!fxmlConfig) {
02844     fxmlConfig = (NuXMLConfig*) ndDataFile.Get("NuXMLConfig");
02845   }
02846   if (!fxmlConfig){
02847     fxmlConfig = (NuXMLConfig*)helperFile.Get("NuXMLConfig");
02848   }
02849 
02850   ndDataFile.Close();
02851 
02852   //Calculate the total PoT
02853   Double_t ndPoT = hndTotalPot->Integral();
02854   Double_t fdPoT = hfdTotalPot->Integral();
02855   if (fxmlConfig->OverriddenFDPoTs()>0) {
02856     fdPoT = fxmlConfig->OverriddenFDPoTs();
02857   }
02858 
02859   // If scaling is used, the scaling should have been done by this
02860   // stage, so use the scaled pots instead of the overridden pots
02861   if (fxmlConfig->ScaledFDPoTs()>0) {
02862     fdPoT = fxmlConfig->ScaledFDPoTs();
02863   }
02864 
02865   MSG("NuMatrixFitter",Msg::kInfo)  << "ndPoT: " << ndPoT
02866        << "; fdPoT: " << fdPoT
02867        << endl;
02868 
02869   NuCuts nuCuts;
02870   Double_t ndFidMass = nuCuts.FiducialMass(Detector::kNear,
02871                                            SimFlag::kMC,
02872                                            NuCuts::kJJE1);
02873   Double_t fdFidMass = nuCuts.FiducialMass(Detector::kFar,
02874                                            SimFlag::kMC,
02875                                            NuCuts::kJJE1);
02876 
02877   MSG("NuMatrixFitter",Msg::kInfo)  << "Near mass: " << ndFidMass
02878        << "; far mas: " << fdFidMass << endl;
02879 
02880   //Binning scheme
02881   NuBinningScheme::NuBinningScheme_t binningScheme =
02882     static_cast<NuBinningScheme::NuBinningScheme_t>(fxmlConfig->BinningScheme());
02883   MSG("NuMatrixFitter",Msg::kInfo)  << "Using binning scheme " << binningScheme << endl;
02884 
02885   NuMatrixMethod mmNuMu(binningScheme,
02886                         ndPoT,
02887                         fdPoT,
02888                         ffluxPoT,
02889                         ndFidMass,
02890                         fdFidMass,
02891                         fhelperFilename,
02892                         fxsecFilename,
02893                         fndDataFilename,
02894                         "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutputNuMu.root");
02895 
02896   //pull in all the histograms from the file
02897   mmNuMu.SetExtrapolationScheme(NuMMExtrapolation::kCCNoChargeCut);
02898   mmNuMu.PredictFDFluxUnosc();
02899   
02900 
02901   //write out files
02902   //create the file
02903   TFile* sfile=0;
02904   if (outputFileName=="") {
02905     sfile=new TFile("$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMNoCutOutputForFitter.root","RECREATE");
02906   }
02907   else {
02908     sfile=new TFile(outputFileName.c_str(),"RECREATE");
02909   }
02910   MSG("NuMatrixFitter",Msg::kInfo) <<"Created file:"<<endl;
02911   sfile->Print();
02912   
02913   fxmlConfig->Write();
02914   
02915   //write files
02916   mmNuMu.WriteFilesForFitter("All");
02917   
02918   //close file
02919   sfile->Close();
02920   if (sfile){delete sfile; sfile = 0;}  
02921 
02922 
02923 /*
02924           const TH1D* numufdPred =
02925             mmNuMu.PredictFDSpectrumBackgroundOscCombined(dm2,
02926                                                           sn2,
02927                                                           dm2bar,
02928                                                           sn2bar);
02929           const TH1D* numubarfdPred =
02930             mmNuMuBar.PredictFDSpectrumBackgroundOscCombined(dm2bar,
02931                                                              sn2bar,
02932                                                              dm2,
02933                                                              sn2);
02934 */
02935 }

void NuMatrixFitter::WritePRLCCHistosForFitter ( std::string  outputFileName  )  [virtual]

Definition at line 2939 of file NuMatrixFitter.cxx.

References NuCuts::FiducialMass(), SimFlag::kData, Detector::kFar, Msg::kInfo, NuCuts::kJJE1, Detector::kNear, NuMMExtrapolation::kPRLCCAnalysis, MSG, NuMatrixMethod::PredictFDFluxUnosc(), NuMatrixMethod::SetExtrapolationScheme(), and NuMatrixMethod::WriteFilesForFitter().

02940 {
02941   //Get the FD data & PoT histograms
02942   //Data:
02943   TFile fdDataFile(ffdDataFilename.c_str(),"READ");
02944   TH1D* numufdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
02945   numufdData->SetDirectory(0);
02946   numufdData->Sumw2();
02947   /*
02948   TH1D* numubarfdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
02949   numubarfdData->SetDirectory(0);
02950   numubarfdData->Sumw2();
02951   */
02952   //PoT:
02953   TH1F* hfdTotalPot = (TH1F*) fdDataFile.Get("hTotalPot");
02954   hfdTotalPot->SetDirectory(0);
02955   fdDataFile.Close();
02956 
02957   //Get the ND PoT histograms
02958   TFile ndDataFile(fndDataFilename.c_str(),"READ");
02959   TH1F* hndTotalPot = (TH1F*) ndDataFile.Get("hTotalPot");
02960   hndTotalPot->SetDirectory(0);
02961   
02962   //Open the helper file.
02963   TFile helperFile(fhelperFilename.c_str(),"READ");
02964 
02965   // Get the NuXMLConfig Object. Look in the ND data file first. If
02966   // it's not in there, look in the helper file.
02967   if (!fxmlConfig) {
02968     fxmlConfig = (NuXMLConfig*) ndDataFile.Get("NuXMLConfig");
02969   }
02970   if (!fxmlConfig){
02971     fxmlConfig = (NuXMLConfig*)helperFile.Get("NuXMLConfig");
02972   }
02973 
02974   ndDataFile.Close();
02975 
02976   //Calculate the total PoT
02977   Double_t ndPoT = hndTotalPot->Integral();
02978   Double_t fdPoT = hfdTotalPot->Integral();
02979   if (fxmlConfig->OverriddenFDPoTs()>0) {
02980     fdPoT = fxmlConfig->OverriddenFDPoTs();
02981   }
02982 
02983   // If scaling is used, the scaling should have been done by this
02984   // stage, so use the scaled pots instead of the overridden pots
02985   if (fxmlConfig->ScaledFDPoTs()>0) {
02986     fdPoT = fxmlConfig->ScaledFDPoTs();
02987   }
02988 
02989   MSG("NuMatrixFitter",Msg::kInfo)  << "ndPoT: " << ndPoT
02990        << "; fdPoT: " << fdPoT
02991        << endl;
02992 
02993   NuCuts nuCuts;
02994   Double_t ndFidMass = nuCuts.FiducialMass(Detector::kNear,
02995                                            SimFlag::kData,
02996                                            NuCuts::kJJE1);
02997   Double_t fdFidMass = nuCuts.FiducialMass(Detector::kFar,
02998                                            SimFlag::kData,
02999                                            NuCuts::kJJE1);
03000   
03001   MSG("NuMatrixFitter",Msg::kInfo)  << "Near mass: " << ndFidMass
03002        << "; far mas: " << fdFidMass << endl;
03003 
03004   //Binning scheme
03005   NuBinningScheme::NuBinningScheme_t binningScheme =
03006     static_cast<NuBinningScheme::NuBinningScheme_t>(fxmlConfig->BinningScheme());
03007   MSG("NuMatrixFitter",Msg::kInfo)  << "Using binning scheme " << binningScheme << endl;
03008 
03009   NuMatrixMethod mmNuMu(binningScheme,
03010                         ndPoT,
03011                         fdPoT,
03012                         ffluxPoT,
03013                         ndFidMass,
03014                         fdFidMass,
03015                         fhelperFilename,
03016                         fxsecFilename,
03017                         fndDataFilename,
03018                         "$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMOutputNuMu.root");
03019 
03020   //pull in all the histograms from the file
03021   mmNuMu.SetExtrapolationScheme(NuMMExtrapolation::kPRLCCAnalysis);
03022   mmNuMu.PredictFDFluxUnosc();
03023   
03024 
03025   //write out files
03026   //create the file
03027   TFile* sfile=0;
03028   if (outputFileName=="") {
03029     sfile=new TFile("$SRT_PRIVATE_CONTEXT/NtupleUtils/results/MMNoCutOutputForFitter.root","RECREATE");
03030   }
03031   else {
03032     sfile=new TFile(outputFileName.c_str(),"RECREATE");
03033   }
03034   MSG("NuMatrixFitter",Msg::kInfo) <<"Created file:"<<endl;
03035   sfile->Print();
03036   
03037   fxmlConfig->Write();
03038   
03039   //write files
03040   mmNuMu.WriteFilesForFitter("NM");
03041   
03042   //close file
03043   sfile->Close();
03044   if (sfile){delete sfile; sfile = 0;}  
03045 
03046 }


Member Data Documentation

Float_t NuMatrixFitter::fdm2BarHigh [private]
Float_t NuMatrixFitter::fdm2BarLow [private]
Float_t NuMatrixFitter::fdm2NuHigh [private]
Float_t NuMatrixFitter::fdm2NuLow [private]
std::string NuMatrixFitter::ffdDataFilename [private]
Double_t NuMatrixFitter::ffluxPoT [private]
std::string NuMatrixFitter::fhelperFilename [private]

Definition at line 378 of file NuMatrixFitter.h.

Referenced by FillGridSearchParams(), and SetGridParamsDefault().

Float_t NuMatrixFitter::fmu2NuHigh [private]

Definition at line 377 of file NuMatrixFitter.h.

Referenced by FillGridSearchParams(), and SetGridParamsDefault().

Float_t NuMatrixFitter::fmu2NuLow [private]

Definition at line 376 of file NuMatrixFitter.h.

Referenced by FillGridSearchParams(), and SetGridParamsDefault().

std::string NuMatrixFitter::fndDataFilename [private]
std::string NuMatrixFitter::foutputFilename [private]

Definition at line 348 of file NuMatrixFitter.h.

Referenced by CCAnalysis(), and NuMatrixFitter().

Float_t NuMatrixFitter::fsn2BarHigh [private]
Float_t NuMatrixFitter::fsn2BarLow [private]
Float_t NuMatrixFitter::fsn2NuHigh [private]
Float_t NuMatrixFitter::fsn2NuLow [private]

Definition at line 374 of file NuMatrixFitter.h.

Referenced by FillGridSearchParams(), and SetGridParamsDefault().

Definition at line 373 of file NuMatrixFitter.h.

Referenced by FillGridSearchParams(), and SetGridParamsDefault().

Definition at line 372 of file NuMatrixFitter.h.

Referenced by FillGridSearchParams(), and SetGridParamsDefault().

std::vector<std::string> NuMatrixFitter::fvfdDataFilenames [private]

Definition at line 351 of file NuMatrixFitter.h.

std::vector<std::string> NuMatrixFitter::fvhelperFilenames [private]

Definition at line 352 of file NuMatrixFitter.h.

std::vector<std::string> NuMatrixFitter::fvndDataFilenames [private]

Definition at line 350 of file NuMatrixFitter.h.

std::string NuMatrixFitter::fxsecFilename [private]

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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1