Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

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 const Double_t CalculateChi2 (const TH1D *fdPred, const TH1D *fdData) const
virtual const Double_t CalculateLikelihood (const TH1D *fdPred, const TH1D *fdData) const
virtual const Double_t CalculateLikelihoodNormPenalty (const TH1D *fdPred, const TH1D *fdData, const Double_t normalisation) const
virtual const Double_t CalculateLikelihood (const TH1D *fdPred, const TH1D *fdData, const Double_t nubarEcut) const
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


Member Enumeration Documentation

enum NuMatrixFitter::NuModel_t
 

Enumeration values:
kUndefined 
kOscillation 
kDecay 
kDecoherence 

Definition at line 203 of file NuMatrixFitter.h.

00203               {
00204     kUndefined = 0,
00205       kOscillation = 1,
00206       kDecay = 2,
00207       kDecoherence = 3
00208   } NuModel_t;


Constructor & Destructor Documentation

NuMatrixFitter::NuMatrixFitter  ) 
 

Definition at line 312 of file NuMatrixFitter.cxx.

00313 {
00314   fndDataFilename = "";
00315   ffdDataFilename = "";
00316   fhelperFilename = "";
00317   fxsecFilename = "";
00318 
00319   SetGridParamsDefault();
00320   fxmlConfig=0;
00321 }

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 324 of file NuMatrixFitter.cxx.

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

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

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 375 of file NuMatrixFitter.cxx.

00376 {
00377 }


Member Function Documentation

void NuMatrixFitter::AppearanceAnalysis  )  [virtual]
 

Definition at line 2174 of file NuMatrixFitter.cxx.

References ffdDataFilename, ffluxPoT, fhelperFilename, fndDataFilename, fxsecFilename, NuMatrixMethod::GetFDOtherCCTrueUnoscBackground(), NuMatrixMethod::GetFDPotentialAppearanceFidFlux(), MSG, NuMatrixMethod::PredictFDFluxUnosc(), NuMatrixMethod::PredictFDSpectrumAppearance(), NuMatrixMethod::SetExtrapolationScheme(), NuMatrixMethod::SetFDAppearedFidFlux(), and NuMatrixMethod::SetFDCCTrueUnoscBackground().

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

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

Definition at line 436 of file NuMatrixFitter.cxx.

00438 {
00439   //Return Chi2.
00440 
00441   Double_t chi2 = 0;
00442 
00443   for (Int_t i=1; i<=fdPred->GetNbinsX(); ++i){
00444     //Bizarre limits because root histograms are silly
00445     Double_t mnu = fdPred->GetBinContent(i);
00446     Double_t dnu = fdData->GetBinContent(i);
00447 //     MSG("NuMatrixFitter",Msg::kInfo)  << "Bin " << i
00448 //       << ": MC " << mnu
00449 //       << "; data " << dnu
00450 //       << endl;
00451     if (mnu<0.0001){mnu = 0.0001;}
00452     chi2 += (dnu-mnu)*(dnu-mnu)/mnu;
00453   }
00454   return chi2;
00455 }

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

Definition at line 515 of file NuMatrixFitter.cxx.

00518 {
00519   //Aim to minimise -2lnL. That's what I'm returning.
00520 
00521   Double_t like = 0;
00522   Int_t i = 1;
00523   if (i>0.0){
00524     i = fdPred->GetXaxis()->FindBin(Ecut);
00525   }
00526   for (; i<=fdPred->GetNbinsX(); ++i){
00527     //Bizarre limits because root histograms are silly
00528     Double_t mnu = fdPred->GetBinContent(i);
00529     Double_t dnu = fdData->GetBinContent(i);
00530     if (mnu<0.0001){mnu = 0.0001;}
00531     if (dnu){
00532       like += 2*(mnu - dnu + dnu*log(dnu/mnu));
00533     }
00534     else{like += 2*(mnu - dnu);}
00535   }
00536   return like;
00537 }

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

Definition at line 458 of file NuMatrixFitter.cxx.

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

00460 {
00461   //Aim to minimise -2lnL. That's what I'm returning.
00462 
00463   Double_t like = 0;
00464 
00465   for (Int_t i=1; i<=fdPred->GetNbinsX(); ++i){
00466     //Bizarre limits because root histograms are silly
00467     Double_t mnu = fdPred->GetBinContent(i);
00468     Double_t dnu = fdData->GetBinContent(i);
00469     if (mnu<0.0001){mnu = 0.0001;}
00470     if (dnu){
00471       like += 2*(mnu - dnu + dnu*log(dnu/mnu));
00472     }
00473     else{like += 2*(mnu - dnu);}
00474   }
00475 //   MSG("NuMatrixFitter",Msg::kInfo)  << "Like = " << like << endl;
00476   return like;
00477 }

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

Definition at line 480 of file NuMatrixFitter.cxx.

00483 {
00484   //Aim to minimise -2lnL. That's what I'm returning.
00485 
00486   Double_t normOneSigma = 0.04;
00487   Double_t like = 0;
00488 
00489   for (Int_t i=1; i<=fdPred->GetNbinsX(); ++i){
00490     //Bizarre limits because root histograms are silly
00491     Double_t mnu = fdPred->GetBinContent(i);
00492     Double_t dnu = fdData->GetBinContent(i);
00493     if (mnu<0.0001){mnu = 0.0001;}
00494     if (dnu){
00495       like += 2*(mnu - dnu + dnu*log(dnu/mnu)) +
00496         ((normalisation - 1.0)*(normalisation-1.0))/
00497         (normOneSigma*normOneSigma);
00498     }
00499     else{like += 2*(mnu - dnu);}
00500   }
00501 //   MSG("NuMatrixFitter",Msg::kInfo)  << "Like = " << like << endl;
00502   return like;
00503 }

void NuMatrixFitter::CCAnalysis  )  [virtual]
 

Definition at line 540 of file NuMatrixFitter.cxx.

References CalculateLikelihood(), fdm2NuGranularity, fdm2NuHigh, fdm2NuLow, ffdDataFilename, ffluxPoT, fhelperFilename, fndDataFilename, foutputFilename, fsn2NuGranularity, fsn2NuHigh, fsn2NuLow, fxsecFilename, MSG, NuMatrixMethod::PredictFDFluxUnosc(), NuMatrixMethod::PredictFDSpectrumBackgroundNoOsc(), and NuMatrixMethod::SetExtrapolationScheme().

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

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]
 

Definition at line 1950 of file NuMatrixFitter.cxx.

References DoCCFitChargeCut(), FillGridSearchParams(), NuXMLConfig::FullTitle(), fxmlConfig, infile, MSG, and NuMatrixOutput::XMLConfig().

01955 {  
01956   TFile fdDataFile(fdDataFilename.c_str(),"READ");
01957   TH1D* numubarfdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
01958   numubarfdData->SetDirectory(0);
01959   numubarfdData->Sumw2();
01960   TH1D* numufdData = (TH1D*) fdDataFile.Get("RecoEnergy_FD");
01961   numufdData->SetDirectory(0);
01962   numufdData->Sumw2();
01963   fdDataFile.Close();
01964   
01965   TFile *infile = new TFile(inputFilename.c_str(), "READ");
01966   
01967   NuMatrixInput mmInput(infile);
01968   
01969   TString objName = "NuMatrixOutput";
01970   if (!fxmlConfig) {
01971     fxmlConfig = (NuXMLConfig*)gROOT->FindObject("NuXMLConfig");
01972   }
01973   
01974   if (fxmlConfig) {
01975     FillGridSearchParams(fxmlConfig);
01976     MSG("NuMatrixFitter",Msg::kInfo)  << "This file fit is for systematic: " << fxmlConfig->FullTitle() << endl;
01977 
01978     //objName = xmlConfig->FullName();
01979   }
01980   else {
01981     MSG("NuMatrixFitter",Msg::kInfo)  << "No NuXMLConfig object -- nominal fit assumed." << endl;
01982     //objName = "Nominal";
01983   }
01984 
01985   NuMatrixOutput *mmOutput = this->DoCCFitChargeCut(mmInput, 
01986                                                     numufdData, 
01987                                                     numubarfdData, 
01988                                                     nubarEcut,
01989                                                     model);
01990 
01991   if (fxmlConfig) {
01992     mmOutput->XMLConfig(fxmlConfig);
01993   }
01994   MSG("NuMatrixFitter",Msg::kInfo)  << "Recreating file " << outputFilename << endl;
01995   TFile *outputFile = new TFile(outputFilename.c_str(),"RECREATE");
01996   outputFile->cd();
01997   mmOutput->Write(objName);
01998   outputFile->Close();
01999   infile->cd();
02000   MSG("NuMatrixFitter",Msg::kInfo)  << "Don't mind the segfault.  Everything is done and it appears to be harmless." << endl;
02001   //    infile->Close();
02002 }

void NuMatrixFitter::CPTAnalysis  )  [virtual]
 

Definition at line 2277 of file NuMatrixFitter.cxx.

References CalculateLikelihood(), ffdDataFilename, ffluxPoT, fhelperFilename, fndDataFilename, fxsecFilename, NuMatrixMethod::GetFDOtherCCTrueUnoscBackground(), MSG, NuMatrixMethod::PredictFDFluxUnosc(), NuMatrixMethod::PredictFDSpectrumBackgroundOscCombined(), NuMatrixMethod::SetExtrapolationScheme(), and NuMatrixMethod::SetFDCCTrueUnoscBackground().

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

void NuMatrixFitter::CPTFit const std::string  inputFilename,
const string  fdDataFilename,
const Double_t  dm2,
const Double_t  sn2,
const std::string  outputFilename
[virtual]
 

Definition at line 1895 of file NuMatrixFitter.cxx.

References DoCPTFit(), FillGridSearchParams(), NuXMLConfig::FullTitle(), fxmlConfig, infile, MSG, and NuMatrixOutput::XMLConfig().

01900 {
01901   TFile fdDataFile(fdDataFilename.c_str(),"READ");
01902   TH1D* numubarfdData = (TH1D*) fdDataFile.Get("RecoEnergyPQ_FD");
01903   numubarfdData->SetDirectory(0);
01904   numubarfdData->Sumw2();
01905   fdDataFile.Close();
01906 
01907   TFile *infile = new TFile(inputFilename.c_str(), "READ");
01908         
01909 
01910   NuMatrixInput *mmInput = 0;
01911   mmInput = (NuMatrixInput*) infile->Get("NuMatrixInput");
01912   if (!mmInput) {
01913     mmInput = new NuMatrixInput(infile);
01914   }
01915   
01916 
01917   
01918   TString objName = "NuMatrixOutput";
01919   if (!fxmlConfig) {
01920     fxmlConfig = (NuXMLConfig*)gROOT->FindObject("NuXMLConfig");
01921   }
01922   
01923   if (fxmlConfig) {
01924     FillGridSearchParams(fxmlConfig);
01925     MSG("NuMatrixFitter",Msg::kInfo)  << "This file fit is for systematic: " << fxmlConfig->FullTitle() << endl;
01926 
01927     //objName = fxmlConfig->FullName();
01928   }
01929   else {
01930     MSG("NuMatrixFitter",Msg::kInfo)  << "No NuXMLConfig object -- nominal fit assumed." << endl;
01931     //objName = "Nominal";
01932   }
01933 
01934   NuMatrixOutput *mmOutput = this->DoCPTFit(mmInput, numubarfdData, dm2, sn2);
01935   if (fxmlConfig) {
01936     mmOutput->XMLConfig(fxmlConfig);
01937   }
01938   
01939   MSG("NuMatrixFitter",Msg::kInfo)  << "Recreating file " << outputFilename << endl;
01940   TFile *outputFile = new TFile(outputFilename.c_str(),"RECREATE");
01941   outputFile->cd();
01942   mmOutput->Write(objName);
01943   outputFile->Close();
01944   infile->cd();
01945   MSG("NuMatrixFitter",Msg::kInfo)  << "Don't mind the segfault.  Everything is done and it appears to be harmless." << endl;
01946   //    infile->Close();
01947 }

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 764 of file NuMatrixFitter.cxx.

References NuMatrixOutput::BestFitPoint(), MSG, NuMatrixOutput::NuMuBarBestFitFDPrediction(), NuMatrixOutput::NuMuBarChi2Surface(), NuMatrixOutput::NuMuBarFDData(), NuMatrixOutput::NuMuBarNoOscPrediction(), NuMatrixOutput::NuMuBestFitFDPrediction(), NuMatrixOutput::NuMuChi2Surface(), NuMatrixOutput::NuMuFDData(), NuMatrixOutput::NuMuNoOscPrediction(), NuMatrixMethod::PredictFDDecaySpectrumBackgroundDecayCombined(), NuMatrixMethod::PredictFDDecoherenceSpectrumBackgroundDecohereCombined(), and NuMatrixMethod::PredictFDSpectrumBackgroundOscCombined().

Referenced by CCFitChargeCut().

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

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

Definition at line 694 of file NuMatrixFitter.cxx.

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

Referenced by CPTFit().

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

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 1489 of file NuMatrixFitter.cxx.

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

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

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

Definition at line 1661 of file NuMatrixFitter.cxx.

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

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

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 1192 of file NuMatrixFitter.cxx.

References NuMatrixOutput::BestFitPoint(), MSG, NuMatrixOutput::NuMuBarBestFitFDPrediction(), NuMatrixOutput::NuMuBarChi2TransSurface(), NuMatrixOutput::NuMuBarFDData(), NuMatrixOutput::NuMuBarNoOscPrediction(), and NuMatrixOutput::NuMuBarNoTransPrediction().

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

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

Definition at line 955 of file NuMatrixFitter.cxx.

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

Referenced by NoChargeCutFit().

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

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

Definition at line 1070 of file NuMatrixFitter.cxx.

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

Referenced by PRLCCFit().

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

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

Definition at line 639 of file NuMatrixFitter.cxx.

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

Referenced by TransitionFit().

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

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

Definition at line 411 of file NuMatrixFitter.cxx.

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

Referenced by CCFitChargeCut(), CPTFit(), NoChargeCutFit(), NuMatrixFitter(), PRLCCFit(), SetXMLConfig(), and TransitionFit().

00412 {
00413    fdm2NuLow = xmlConfig->DM2NuLow();
00414    fdm2NuHigh = xmlConfig->DM2NuHigh();
00415    fdm2NuGranularity = xmlConfig->DM2NuGranularity();
00416 
00417    fsn2NuLow = xmlConfig->SN2NuLow();
00418    fsn2NuHigh = xmlConfig->SN2NuHigh();
00419    //MSG("NuMatrixFitter",Msg::kInfo)  << "SN2NuGranularity = " << xmlConfig->SN2NuGranularity() << endl;
00420    fsn2NuGranularity = xmlConfig->SN2NuGranularity();
00421 
00422    fdm2BarLow = xmlConfig->DM2BarLow();
00423    fdm2BarHigh = xmlConfig->DM2BarHigh();
00424    fdm2BarGranularity = xmlConfig->DM2BarGranularity();
00425 
00426    fsn2BarLow = xmlConfig->SN2BarLow();
00427    fsn2BarHigh = xmlConfig->SN2BarHigh();
00428    fsn2BarGranularity = xmlConfig->SN2BarGranularity();
00429 
00430    ftransitionProbLow = xmlConfig->TransitionProbLow();
00431    ftransitionProbHigh = xmlConfig->TransitionProbHigh();
00432    ftransitionProbGranularity = xmlConfig->TransitionProbGranularity();
00433 }

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]
 

void NuMatrixFitter::NoChargeCutFit const std::string  inputFilename,
const string  fdDataFilename,
const std::string  outputFilename,
const NuModel_t  model = kOscillation
[virtual]
 

Definition at line 2005 of file NuMatrixFitter.cxx.

References DoNoChargeCutFit(), FillGridSearchParams(), NuXMLConfig::FullTitle(), fxmlConfig, infile, MSG, and NuMatrixOutput::XMLConfig().

02009 {  
02010   TFile fdDataFile(fdDataFilename.c_str(),"READ");
02011   TH1D* fdData = (TH1D*) fdDataFile.Get("RecoEnergyAll_FD");
02012   fdData->SetDirectory(0);
02013   fdData->Sumw2();
02014   fdDataFile.Close();
02015   
02016   TFile *infile = new TFile(inputFilename.c_str(), "READ");
02017   
02018   NuMatrixInput mmInput(infile);
02019 
02020   
02021   TString objName = "NuMatrixOutput";
02022   if (!fxmlConfig) {
02023     fxmlConfig = (NuXMLConfig*)gROOT->FindObject("NuXMLConfig");
02024   }
02025   
02026   if (fxmlConfig) {
02027     FillGridSearchParams(fxmlConfig);
02028     MSG("NuMatrixFitter",Msg::kInfo)  << "This file fit is for systematic: " << fxmlConfig->FullTitle() << endl;
02029 
02030     //objName = fxmlConfig->FullName();
02031   }
02032   else {
02033     MSG("NuMatrixFitter",Msg::kInfo)  << "No NuXMLConfig object -- nominal fit assumed." << endl;
02034     //objName = "Nominal";
02035   }
02036 
02037   NuMatrixOutput *mmOutput = this->DoNoChargeCutFit(mmInput, fdData, model);
02038 
02039   if (fxmlConfig) {
02040     mmOutput->XMLConfig(fxmlConfig);
02041   }
02042   MSG("NuMatrixFitter",Msg::kInfo)  << "Recreating file " << outputFilename << endl;
02043   TFile *outputFile = new TFile(outputFilename.c_str(),"RECREATE");
02044   outputFile->cd();
02045   mmOutput->Write(objName);
02046   outputFile->Close();
02047   infile->cd();
02048   MSG("NuMatrixFitter",Msg::kInfo)  << "Don't mind the segfault.  Everything is done and it appears to be harmless." << endl;
02049   //    infile->Close();
02050 }

const Double_t NuMatrixFitter::NormalisationPenaltyTerm const Double_t  normalisation  )  const
 

Definition at line 506 of file NuMatrixFitter.cxx.

Referenced by DoMultiRunPRLCCFit().

00507 {
00508   Double_t normOneSigma = 0.04;
00509   return ((normalisation - 1.0)*(normalisation-1.0))/
00510     (normOneSigma*normOneSigma);
00511 }

void NuMatrixFitter::NuMuBarAppearanceAnalysis  )  [virtual]
 

Definition at line 3037 of file NuMatrixFitter.cxx.

03038 {
03039 
03040 }

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

Definition at line 2053 of file NuMatrixFitter.cxx.

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

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

TH1D NuMatrixFitter::RebinForFit const TH1D *  fineHist  ) 
 

Definition at line 1632 of file NuMatrixFitter.cxx.

References Form().

Referenced by DoMultiRunPRLCCFit().

01633 {
01634 //   NuBinningScheme::NuBinningScheme_t binningScheme = NuBinningScheme::k0_5GeVTo200;
01635 //   const NuUtilities utils;
01636 //   const vector<Double_t> recoBins = utils.RecoBins(binningScheme);
01637 //   const Int_t numFitBins = recoBins.size()-1;
01638 //   Double_t *fitBinsArray;
01639 //   fitBinsArray=new Double_t[numFitBins+1];
01640 //   {
01641 //     Int_t i=0;
01642 //     for (vector<Double_t>::const_iterator itBin = recoBins.begin();
01643 //       itBin != recoBins.end();
01644 //       ++itBin, ++i){
01645 //       fitBinsArray[i] = *itBin;
01646 //     }
01647 //   }
01648 
01649   static Int_t stupidRoot = 0;
01650   TH1D newHist(Form("NewHist%d",stupidRoot++),"",800,0,200);
01651   newHist.SetBinContent(1,0);
01652   for (Int_t oldBin = 1; oldBin<fineHist->GetNbinsX()+1; ++oldBin){
01653     newHist.SetBinContent(oldBin+1,fineHist->GetBinContent(oldBin));
01654   }
01655   newHist.Rebin(2);
01656   return newHist;
01657 }

void NuMatrixFitter::SetGridParamsDefault  )  [private]
 

Definition at line 387 of file NuMatrixFitter.cxx.

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

Referenced by NuMatrixFitter().

00388 {
00389   fdm2NuLow=0.001;
00390   fdm2NuHigh=0.01;
00391   fdm2NuGranularity=0.001;
00392 
00393   fsn2NuLow=0;
00394   fsn2NuHigh=1;
00395   fsn2NuGranularity=0.05;
00396 
00397   fdm2BarLow=0.001;
00398   fdm2BarHigh=0.01;
00399   fdm2BarGranularity=0.001;
00400 
00401   fsn2BarLow=0;
00402   fsn2BarHigh=1;
00403   fsn2BarGranularity=0.05;
00404 
00405   ftransitionProbLow=0;
00406   ftransitionProbHigh=1;
00407   ftransitionProbGranularity=0.05;
00408 
00409 }

void NuMatrixFitter::SetXMLConfig const NuXMLConfig xmlConfig  ) 
 

Definition at line 380 of file NuMatrixFitter.cxx.

References FillGridSearchParams(), and fxmlConfig.

00381 {
00382   fxmlConfig=xmlConfig;
00383   this->FillGridSearchParams(xmlConfig);
00384 }

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

Definition at line 2101 of file NuMatrixFitter.cxx.

References NuMatrixOutput::bestTransitionProb, NuXMLConfig::DM2Nu(), DoTransitionFit(), FillGridSearchParams(), NuFluctuator::Fluctuate(), NuXMLConfig::FullTitle(), fxmlConfig, infile, MSG, NuXMLConfig::SN2Nu(), and NuMatrixOutput::XMLConfig().

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

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

Definition at line 2521 of file NuMatrixFitter.cxx.

References NuXMLConfig::BinningScheme(), ffdDataFilename, ffluxPoT, fhelperFilename, NuCuts::FiducialMass(), fndDataFilename, fxmlConfig, fxsecFilename, NuMatrixMethod::GetFDOtherCCTrueUnoscBackground(), NuMatrixMethod::GetFDPotentialAppearanceFidFlux(), MSG, NuXMLConfig::OverriddenFDPoTs(), NuMatrixMethod::PredictFDFluxUnosc(), NuMatrixMethod::PredictFDSpectrumBackgroundOscCombined(), NuXMLConfig::ScaledFDPoTs(), NuMatrixMethod::SetExtrapolationScheme(), NuMatrixMethod::SetFDAppearedFidFlux(), NuMatrixMethod::SetFDCCTrueUnoscBackground(), and NuMatrixMethod::WriteFilesForFitter().

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

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

Definition at line 2678 of file NuMatrixFitter.cxx.

References NuXMLConfig::BinningScheme(), MSG, NuMatrixMethod::PredictFDFluxUnosc(), NuMatrixMethod::PredictFDSpectrumBackgroundOscCombined(), NuMatrixMethod::SetExtrapolationScheme(), NuMatrixMethod::SetFDAppearedFidFlux(), NuMatrixMethod::SetFDCCTrueUnoscBackground(), and NuMatrixMethod::WriteFilesForFitter().

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

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

Definition at line 2804 of file NuMatrixFitter.cxx.

References NuCuts::FiducialMass(), MSG, NuMatrixMethod::PredictFDFluxUnosc(), NuMatrixMethod::SetExtrapolationScheme(), and NuMatrixMethod::WriteFilesForFitter().

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

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

Definition at line 2926 of file NuMatrixFitter.cxx.

References NuCuts::FiducialMass(), MSG, NuMatrixMethod::PredictFDFluxUnosc(), NuMatrixMethod::SetExtrapolationScheme(), and NuMatrixMethod::WriteFilesForFitter().

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


Member Data Documentation

Float_t NuMatrixFitter::fdm2BarGranularity [private]
 

Definition at line 348 of file NuMatrixFitter.h.

Referenced by DoCPTFit(), DoMultiRunCPTFit(), FillGridSearchParams(), and SetGridParamsDefault().

Float_t NuMatrixFitter::fdm2BarHigh [private]
 

Definition at line 347 of file NuMatrixFitter.h.

Referenced by DoCPTFit(), DoMultiRunCPTFit(), FillGridSearchParams(), and SetGridParamsDefault().

Float_t NuMatrixFitter::fdm2BarLow [private]
 

Definition at line 346 of file NuMatrixFitter.h.

Referenced by DoCPTFit(), DoMultiRunCPTFit(), FillGridSearchParams(), and SetGridParamsDefault().

Float_t NuMatrixFitter::fdm2NuGranularity [private]
 

Definition at line 340 of file NuMatrixFitter.h.

Referenced by CCAnalysis(), DoMultiRunPRLCCFit(), DoNoChargeCutFit(), DoPRLCCFit(), FillGridSearchParams(), and SetGridParamsDefault().

Float_t NuMatrixFitter::fdm2NuHigh [private]
 

Definition at line 339 of file NuMatrixFitter.h.

Referenced by CCAnalysis(), DoMultiRunPRLCCFit(), DoNoChargeCutFit(), DoPRLCCFit(), FillGridSearchParams(), and SetGridParamsDefault().

Float_t NuMatrixFitter::fdm2NuLow [private]
 

Definition at line 338 of file NuMatrixFitter.h.

Referenced by CCAnalysis(), DoMultiRunPRLCCFit(), DoNoChargeCutFit(), DoPRLCCFit(), FillGridSearchParams(), and SetGridParamsDefault().

std::string NuMatrixFitter::ffdDataFilename [private]
 

Definition at line 327 of file NuMatrixFitter.h.

Referenced by AppearanceAnalysis(), CCAnalysis(), CPTAnalysis(), NuMatrixFitter(), and WriteHistosForFitter().

Double_t NuMatrixFitter::ffluxPoT [private]
 

Definition at line 325 of file NuMatrixFitter.h.

Referenced by AppearanceAnalysis(), CCAnalysis(), CPTAnalysis(), NuMatrixFitter(), and WriteHistosForFitter().

std::string NuMatrixFitter::fhelperFilename [private]
 

Definition at line 328 of file NuMatrixFitter.h.

Referenced by AppearanceAnalysis(), CCAnalysis(), CPTAnalysis(), NuMatrixFitter(), and WriteHistosForFitter().

std::string NuMatrixFitter::fndDataFilename [private]
 

Definition at line 326 of file NuMatrixFitter.h.

Referenced by AppearanceAnalysis(), CCAnalysis(), CPTAnalysis(), NuMatrixFitter(), and WriteHistosForFitter().

std::string NuMatrixFitter::foutputFilename [private]
 

Definition at line 330 of file NuMatrixFitter.h.

Referenced by CCAnalysis(), and NuMatrixFitter().

Float_t NuMatrixFitter::fsn2BarGranularity [private]
 

Definition at line 352 of file NuMatrixFitter.h.

Referenced by DoCPTFit(), DoMultiRunCPTFit(), FillGridSearchParams(), and SetGridParamsDefault().

Float_t NuMatrixFitter::fsn2BarHigh [private]
 

Definition at line 351 of file NuMatrixFitter.h.

Referenced by DoCPTFit(), DoMultiRunCPTFit(), FillGridSearchParams(), and SetGridParamsDefault().

Float_t NuMatrixFitter::fsn2BarLow [private]
 

Definition at line 350 of file NuMatrixFitter.h.

Referenced by DoCPTFit(), DoMultiRunCPTFit(), FillGridSearchParams(), and SetGridParamsDefault().

Float_t NuMatrixFitter::fsn2NuGranularity [private]
 

Definition at line 344 of file NuMatrixFitter.h.

Referenced by CCAnalysis(), DoMultiRunPRLCCFit(), DoNoChargeCutFit(), DoPRLCCFit(), FillGridSearchParams(), and SetGridParamsDefault().

Float_t NuMatrixFitter::fsn2NuHigh [private]
 

Definition at line 343 of file NuMatrixFitter.h.

Referenced by CCAnalysis(), DoMultiRunPRLCCFit(), DoNoChargeCutFit(), DoPRLCCFit(), FillGridSearchParams(), and SetGridParamsDefault().

Float_t NuMatrixFitter::fsn2NuLow [private]
 

Definition at line 342 of file NuMatrixFitter.h.

Referenced by CCAnalysis(), DoMultiRunPRLCCFit(), DoNoChargeCutFit(), DoPRLCCFit(), FillGridSearchParams(), and SetGridParamsDefault().

Float_t NuMatrixFitter::ftransitionProbGranularity [private]
 

Definition at line 356 of file NuMatrixFitter.h.

Referenced by FillGridSearchParams(), and SetGridParamsDefault().

Float_t NuMatrixFitter::ftransitionProbHigh [private]
 

Definition at line 355 of file NuMatrixFitter.h.

Referenced by FillGridSearchParams(), and SetGridParamsDefault().

Float_t NuMatrixFitter::ftransitionProbLow [private]
 

Definition at line 354 of file NuMatrixFitter.h.

Referenced by FillGridSearchParams(), and SetGridParamsDefault().

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

Definition at line 333 of file NuMatrixFitter.h.

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

Definition at line 334 of file NuMatrixFitter.h.

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

Definition at line 332 of file NuMatrixFitter.h.

const NuXMLConfig* NuMatrixFitter::fxmlConfig [private]
 

Definition at line 336 of file NuMatrixFitter.h.

Referenced by CCFitChargeCut(), CPTFit(), NoChargeCutFit(), NuMatrixFitter(), PRLCCFit(), SetXMLConfig(), TransitionFit(), and WriteHistosForFitter().

std::string NuMatrixFitter::fxsecFilename [private]
 

Definition at line 329 of file NuMatrixFitter.h.

Referenced by AppearanceAnalysis(), CCAnalysis(), CPTAnalysis(), NuMatrixFitter(), and WriteHistosForFitter().


The documentation for this class was generated from the following files:
Generated on Mon Feb 15 11:09:53 2010 for loon by  doxygen 1.3.9.1