NCExtrapolationBeamMatrix Class Reference

Extrapolation using the BeamMatrix method. More...

#include <NCExtrapolationBeamMatrix.h>

Inheritance diagram for NCExtrapolationBeamMatrix:
NCExtrapolation

List of all members.

Public Member Functions

 NCExtrapolationBeamMatrix ()
virtual ~NCExtrapolationBeamMatrix ()
void Prepare (const Registry &r)
void WriteResources (const NC::OscProb::OscPars *trueOscPars)
virtual TString GetShortName () const
 This is the name used to name things in the output file etc.
virtual TString GetLongName () const
 This is the name the extrapolation is known under on plots and such.
virtual bool EnableNearToFarExtrapolation (bool enable)
 Enable or disable corrections to FD spectra based on ND spectra.
std::vector< const TH1 * > GetNearMCSpectra (std::vector< NCBeam::Info > beamsToUse)

Static Public Member Functions

static const RegistryDefaultConfig ()
 Return a default config that will be merged with the NCExtrapolationModule DefaultConfig.

Private Types

typedef std::map< NCBeam::Info,
int > 
BeamMatrixMap

Private Member Functions

void ConstructFarSpectrum (const NC::OscProb::OscPars *pars, BeamMatrix *bmtouse)
void FillDataMCHistogramsNear (NCBeam *beam, NCType::EEventType nccc, BeamMatrix *bmtouse, bool systFit)
void FillDataMCHistogramsFar (NCBeam *beam, NCType::EEventType nccc, BeamMatrix *bmtouse)
void FillResultHistograms (NCBeam::Info beamInfo, NCType::EEventType nccc, std::vector< std::vector< double > > &nominal, std::vector< std::vector< double > > &osc)
void DoneFilling ()
 Called after the final AddEvent() but before any calls to FindSpectraForPars().
void DoPurityCorrectionND (BeamMatrix *bmtouse)
void DoRecoToTrueND (BeamMatrix *bmtouse, int beamindex=-1)
void DoEfficencyCorrectionND (BeamMatrix *bmtouse)
void DoXsectionCorrectionND (BeamMatrix *bmtouse)
void DoBeamMatrixExtrapolation (BeamMatrix *bmtouse)
void DoXsectionCorrectionFD (BeamMatrix *bmtouse)
void DoEfficencyCorrectionFD (BeamMatrix *bmtouse)
void DoUnoscTrueToRecoFD (BeamMatrix *bmtouse)
void DoUnoscPurityCorrectionFD (BeamMatrix *bmtouse)
void DoOscillations (const NC::OscProb::OscPars *pars, BeamMatrix *bmtouse)
void GetNCFromCCFlux (BeamMatrix *bmtouse)
void SetFiducialMasses (NCType::ECuts cutsType)
void rebinRecoTrueToRecoRecoHist (TH2D *t2r, TH2D *t2rRebin)
void normaliseRecoToTrue (TH2D *t2rRebin, TH1D *norm)
NC::RunUtil::ERunType GetRunTypeFromBM (BeamMatrix *BM)
void GetEfficiencyHistograms (NCBeam *beam, BeamMatrix *bmtouse)
void GetBeamMatrix (NCBeam *beam, BeamMatrix *bmtouse)
template<class T >
void HistSmoothHelper (T *h, bool is2D, double weight, double x, double y=-1)
void FillHistogramSmoothed (TH1D *h, double x, double weight)
void FillHistogramSmoothed2D (TH2D *h, double x, double y, double weight)
virtual void FindSpectraForPars (const NC::OscProb::OscPars *oscPars, const NC::SystPars &systPars, std::vector< NCBeam::Info > beamsToUse, std::vector< TH1 * > &exps, std::vector< TH1 * > &obss)
 Override this in the derived class.
virtual void CleanupSpectra (std::vector< TH1 * > exps, std::vector< TH1 * > obss)
 Called after FindSpectraForPars() to delete necessary spectra.
void doNDXsectionFit (int PrintLevel, NCBeam *beam, int ibm)
void Fit (int PrintLevel)
int WhichFitParam (double trueEnergy)
int WhichFluxCorrParam (double trueEnergy)
void CopyNDDataForSystFit (BeamMatrix *bmtouse, bool doSyst, int bmind)
void FillNDHistsForXSectionFit (NCBeam *beam_nd, NCType::EEventType nccc, bool doSyst)
void ScaleBeamsToSameNumberOfEvents (int beamType, bool doSyst)
void correctNDNCFlux (BeamMatrix *bmtouse, int beamIndex)
void GetNCFromCCFluxND (BeamMatrix *bmtouse)
void DoEfficencyCorrectionForFluxCorrND (BeamMatrix *bmtouse)
void makeFluxCorrectionVector (BeamMatrix *bmtouse, int beamindex)

Static Private Member Functions

static void LogLikelihoodFunc (Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag)

Private Attributes

double fNDFidVolMass
double fFDFidVolMass
double fNDPOTData
double fFDPOTData
double fNDPOTMC
double fFDPOTMC
double fFluxPOT
TString fSelEffFilePath
TH1D * dataNDhistCCLE
TH1D * dataNDhistCCME
TH1D * dataNDhistCCHE
TH1D * dataNDhistNCLE
TH1D * dataNDhistNCME
TH1D * dataNDhistNCHE
TH2D * t2rCCnd
TH2D * t2rNCnd
TH2D * t2rCC
TH2D * t2rNC
TH2D * t2rCCRebin
TH2D * t2rNCRebin
TString fBeamMatrixPath
const TH2D * fBeamMatrix
int fNCalls
 Counts number of calls to the method.
std::vector< std::vector
< double > > 
fNominal
std::vector< std::vector
< double > > 
fOsc
std::vector< std::vector
< double > > 
fNominalBM
std::vector< std::vector
< double > > 
fOscBM
double fTrue_bin_width
std::vector< double > true_bins
double fUE3Sqr
std::vector< BeamMatrix * > BMbybeams
BeamMatrixMap fBMbybeams
std::vector< NCBeam * > beamsNear
std::vector< std::vector< TH1D * > > ND_XsecByParams
std::vector< std::vector< TH1D * > > dataNDhist
std::vector< TH1D * > ND_XSectionFit
std::vector< TH1D * > ND_XSectionFitRew
std::vector< NCBeam * > beamsFit
std::vector< std::vector
< BeamMatrix * > > 
BMforFitTypes
std::vector< BeamMatrix * > BMforFit
std::vector< Double_t > mInputScale
std::vector< Double_t > mScaleNC
std::vector< Double_t > mScaleNCErr
std::vector< Double_t > m0ScaleNC
std::vector< Double_t > m0ScaleNCErr
std::vector< std::vector
< Double_t > > 
fluxCorrectionsNDNC
std::vector< std::vector
< Double_t > > 
fluxCorrectionsNDNCbkg
std::vector< Double_t > fluxCorrectionsFDCCinNC
double mScaleNeg [4]
double mScalePos [4]
double corr [4]
std::ofstream * logFile
int fExtractionType
TH1D * tempExtractNCfromCCflux
TH1D * fFD_data [2]
double fSmoothWidth
bool fDoExtrapolation

Static Private Attributes

static NCExtrapolationBeamMatrixbmFit = NULL

Detailed Description

Extrapolation using the BeamMatrix method.

Definition at line 32 of file NCExtrapolationBeamMatrix.h.


Member Typedef Documentation

typedef std::map<NCBeam::Info, int> NCExtrapolationBeamMatrix::BeamMatrixMap [private]

Definition at line 183 of file NCExtrapolationBeamMatrix.h.


Constructor & Destructor Documentation

NCExtrapolationBeamMatrix::NCExtrapolationBeamMatrix (  ) 

Definition at line 218 of file NCExtrapolationBeamMatrix.cxx.

References bmFit, dataNDhist, fDoExtrapolation, fluxCorrectionsFDCCinNC, fluxCorrectionsNDNC, fluxCorrectionsNDNCbkg, fNCalls, fSmoothWidth, kEnergyBinsFar, kNumEnergyBinsFar, kNumTrueEnergyBins, ND_XsecByParams, ND_XSectionFit, ND_XSectionFitRew, t2rCC, t2rCCnd, t2rCCRebin, t2rNC, t2rNCnd, t2rNCRebin, and true_bins.

00219 {
00220   fNCalls = 0;
00221   fDoExtrapolation = true;
00222   
00223   int numBins = kNumEnergyBinsFar;
00224   double bins[kNumEnergyBinsFar+1];
00225   
00226   
00227   //use NCType to fill array of bins as it has the right binning already
00228   for( int i = 0; i < numBins+1; ++i ){
00229     bins[i] = kEnergyBinsFar[i];
00230   }
00231   
00232   // Create a true energy binning scheme for the
00233   // FD prediction to oscillate
00234   
00235   //1000
00236   for( int i = 0; i < 200; ++i )
00237     true_bins.push_back(i*0.01);
00238   for( int i = 80; i < 801; ++i )
00239     true_bins.push_back(i*0.025);
00240   for( int i = 41; i < 81; ++i )
00241     true_bins.push_back(i*0.5);
00242   for( int i = 42; i < 120; i+=2 )
00243     true_bins.push_back(i);
00244   true_bins.push_back(120) ;
00245   
00246   
00247   // needed for ND fitting
00248   std::vector< TH1D* > NDByEnergyForXsecLE; 
00249   NDByEnergyForXsecLE.push_back(new TH1D(" NC0to4LE" , " NC0to4LE" , numBins, bins) );
00250   NDByEnergyForXsecLE.push_back(new TH1D(" NC5to8LE" , " NC5to8LE" , numBins, bins) );
00251   NDByEnergyForXsecLE.push_back(new TH1D(" NC8to15LE" , "NC8to15LE" , numBins, bins) );
00252   NDByEnergyForXsecLE.push_back(new TH1D(" NCover15LE" , "NCover15LE" , numBins, bins) );
00253   NDByEnergyForXsecLE.push_back(new TH1D(" NCnotMuLE" , "NCnotMuLE" , numBins, bins) );
00254 
00255   std::vector< TH1D* > NDByEnergyForXsecME; 
00256   NDByEnergyForXsecME.push_back(new TH1D(" NC0to4ME" , " NC0to4ME" , numBins, bins) );
00257   NDByEnergyForXsecME.push_back(new TH1D(" NC5to8ME" , " NC5to8ME" , numBins, bins) );
00258   NDByEnergyForXsecME.push_back(new TH1D(" NC8to15ME" , "NC8to15ME" , numBins, bins) );
00259   NDByEnergyForXsecME.push_back(new TH1D(" NCover15ME" , "NCover15ME" , numBins, bins) );
00260   NDByEnergyForXsecME.push_back(new TH1D(" NCnotMuME" , "NCnotMuME" , numBins, bins) );
00261   
00262   std::vector< TH1D* > NDByEnergyForXsecHE; 
00263   NDByEnergyForXsecHE.push_back(new TH1D(" NC0to4HE" , " NC0to4HE" , numBins, bins) );
00264   NDByEnergyForXsecHE.push_back(new TH1D(" NC5to8HE" , " NC5to8HE" , numBins, bins) );
00265   NDByEnergyForXsecHE.push_back(new TH1D(" NC8to15HE" , "NC8to15HE" , numBins, bins) );
00266   NDByEnergyForXsecHE.push_back(new TH1D(" NCover15HE" , "NCover15HE" , numBins, bins) );
00267   NDByEnergyForXsecHE.push_back(new TH1D(" NCnotMuHE" , "NCnotMuHE" , numBins, bins) );
00268   
00269   ND_XsecByParams.push_back(NDByEnergyForXsecLE);
00270   ND_XsecByParams.push_back(NDByEnergyForXsecME);
00271   ND_XsecByParams.push_back(NDByEnergyForXsecHE);
00272 
00273   std::vector<Double_t> fluxCorrectionsNDLE;  
00274   std::vector<Double_t> fluxCorrectionsNDME;  
00275   std::vector<Double_t> fluxCorrectionsNDHE;  
00276 
00277   std::vector<Double_t> fluxCorrectionsNDLEbkg;  
00278   std::vector<Double_t> fluxCorrectionsNDMEbkg;  
00279   std::vector<Double_t> fluxCorrectionsNDHEbkg;  
00280 
00281   for(int i = 0 ; i < numBins; i++){
00282     fluxCorrectionsNDLE.push_back(1.0);
00283     fluxCorrectionsNDME.push_back(1.0);
00284     fluxCorrectionsNDHE.push_back(1.0);
00285     fluxCorrectionsNDLEbkg.push_back(1.0);
00286     fluxCorrectionsNDMEbkg.push_back(1.0);
00287     fluxCorrectionsNDHEbkg.push_back(1.0);
00288     fluxCorrectionsFDCCinNC.push_back(1.0);
00289   }
00290 
00291   fluxCorrectionsNDNC.push_back(fluxCorrectionsNDLE);
00292   fluxCorrectionsNDNC.push_back(fluxCorrectionsNDME);
00293   fluxCorrectionsNDNC.push_back(fluxCorrectionsNDHE);
00294 
00295   fluxCorrectionsNDNCbkg.push_back(fluxCorrectionsNDLEbkg);
00296   fluxCorrectionsNDNCbkg.push_back(fluxCorrectionsNDMEbkg);
00297   fluxCorrectionsNDNCbkg.push_back(fluxCorrectionsNDHEbkg);
00298 
00299   ND_XSectionFitRew.push_back(new TH1D("NCDataLE" , "NCDataLE" , numBins, bins) );
00300   ND_XSectionFitRew.push_back(new TH1D("NCMCLE" , "NCMCLE" , numBins, bins) );
00301   ND_XSectionFitRew.push_back(new TH1D("NCDataME" , "NCDataME" , numBins, bins) );
00302   ND_XSectionFitRew.push_back(new TH1D("NCMCME" , "NCMCME" , numBins, bins) );
00303   ND_XSectionFitRew.push_back(new TH1D("NCDataHE" , "NCDataHE" , numBins, bins) );
00304   ND_XSectionFitRew.push_back(new TH1D("NCMCHE" , "NCMCHE" , numBins, bins) );
00305   ND_XSectionFitRew.push_back(new TH1D("CCDataLE" , "CCDataLE" , numBins, bins) );
00306   ND_XSectionFitRew.push_back(new TH1D("CCMCLE" , "CCMCLE" , numBins, bins) );
00307   ND_XSectionFitRew.push_back(new TH1D("CCDataME" , "CCDataME" , numBins, bins) );
00308   ND_XSectionFitRew.push_back(new TH1D("CCMCME" , "CCMCME" , numBins, bins) );
00309   ND_XSectionFitRew.push_back(new TH1D("CCDataHE" , "CCDataHE" , numBins, bins) );
00310   ND_XSectionFitRew.push_back(new TH1D("CCMCHE" , "CCMCHE" , numBins, bins) );
00311 
00312   ND_XSectionFit.push_back(new TH1D("NCDataLEOrig" , "NCDataLEOrig" , numBins, bins) );
00313   ND_XSectionFit.push_back(new TH1D("NCMCLEOrig" , "NCMCLEOrig" , numBins, bins) );
00314   ND_XSectionFit.push_back(new TH1D("NCDataMEOrig" , "NCDataMEOrig" , numBins, bins) );
00315   ND_XSectionFit.push_back(new TH1D("NCMCMEOrig" , "NCMCMEOrig" , numBins, bins) );
00316   ND_XSectionFit.push_back(new TH1D("NCDataHEOrig" , "NCDataHEOrig" , numBins, bins) );
00317   ND_XSectionFit.push_back(new TH1D("NCMCHEOrig" , "NCMCHEOrig" , numBins, bins) );
00318   ND_XSectionFit.push_back(new TH1D("CCDataLEOrig" , "CCDataLEOrig" , numBins, bins) );
00319   ND_XSectionFit.push_back(new TH1D("CCMCLEOrig" , "CCMCLEOrig" , numBins, bins) );
00320   ND_XSectionFit.push_back(new TH1D("CCDataMEOrig" , "CCDataMEOrig" , numBins, bins) );
00321   ND_XSectionFit.push_back(new TH1D("CCMCMEOrig" , "CCMCMEOrig" , numBins, bins) );
00322   ND_XSectionFit.push_back(new TH1D("CCDataHEOrig" , "CCDataHEOrig" , numBins, bins) );
00323   ND_XSectionFit.push_back(new TH1D("CCMCHEOrig" , "CCMCHEOrig" , numBins, bins) );
00324 
00325 
00326   std::vector< TH1D* > dataNDhistLE; 
00327   dataNDhistLE.push_back(new TH1D("dataNDhistNCLE","dataNDhistNCLE", numBins, bins));
00328   dataNDhistLE.push_back(new TH1D("dataNDhistCCLE","dataNDhistCCLE", numBins, bins));
00329 
00330   std::vector< TH1D* > dataNDhistME; 
00331   dataNDhistME.push_back(new TH1D("dataNDhistNCME","dataNDhistNCME", numBins, bins));
00332   dataNDhistME.push_back(new TH1D("dataNDhistCCME","dataNDhistCCME", numBins, bins));
00333 
00334   std::vector< TH1D* > dataNDhistHE; 
00335   dataNDhistHE.push_back(new TH1D("dataNDhistCCHE","dataNDhistCCHE", numBins, bins));
00336   dataNDhistHE.push_back(new TH1D("dataNDhistNCHE","dataNDhistNCHE", numBins, bins));
00337 
00338   dataNDhist.push_back(dataNDhistLE);
00339   dataNDhist.push_back(dataNDhistME);
00340   dataNDhist.push_back(dataNDhistHE);
00341   
00342   //temporary histograms needed for true to reco 
00343   t2rCCnd = new TH2D("t2rCCnd","t2rCCnd", numBins, bins, numBins, bins );
00344   t2rNCnd = new TH2D("t2rNCnd","t2rNCnd", numBins, bins, numBins, bins );
00345 
00346   t2rCC = new TH2D("t2rCC","t2rCC",kNumTrueEnergyBins, &true_bins.at(0) , numBins, bins );
00347   t2rNC = new TH2D("t2rNC","t2rNC",kNumTrueEnergyBins, &true_bins.at(0) , numBins, bins );
00348 
00349   t2rCCRebin = new TH2D("t2rCCRebin","t2rCCRebin",numBins, bins , numBins, bins );
00350   t2rNCRebin = new TH2D("t2rNCRebin","t2rNCRebin",numBins, bins , numBins, bins );
00351   
00352   fSmoothWidth = 0; // Don't do any smoothing
00353   
00354   //needed for x section fitting
00355   bmFit = this ;
00356 }

NCExtrapolationBeamMatrix::~NCExtrapolationBeamMatrix (  )  [virtual]

Definition at line 362 of file NCExtrapolationBeamMatrix.cxx.

00362 {}


Member Function Documentation

virtual void NCExtrapolationBeamMatrix::CleanupSpectra ( std::vector< TH1 * >  exps,
std::vector< TH1 * >  obss 
) [private, virtual]

Called after FindSpectraForPars() to delete necessary spectra.

Parameters:
exps The expected spectra just returned from FindSpectraForPars()
obss The observed spectra just returned from FindSpectraForPars()

Reimplemented from NCExtrapolation.

void NCExtrapolationBeamMatrix::ConstructFarSpectrum ( const NC::OscProb::OscPars pars,
BeamMatrix bmtouse 
) [private]
void NCExtrapolationBeamMatrix::CopyNDDataForSystFit ( BeamMatrix bmtouse,
bool  doSyst,
int  bmind 
) [private]

Definition at line 2444 of file NCExtrapolationBeamMatrix.cxx.

References dataNDhist, BeamMatrix::GetDataPredHists(), kCCDataRecoEnergy_ND, and kNCDataRecoEnergy_ND.

Referenced by ScaleBeamsToSameNumberOfEvents().

02444                                                                  {
02445 
02446   // this is the nominal beam and therefore has data in it. Copy to relevant dataNDHist
02447 
02448   if(!doSyst)
02449     for(int i = 1 ; i < kNumEnergyBinsFar+1 ; i++){
02450       dataNDhist[bmind][0]->SetBinContent(i,bmtouse->GetDataPredHists().at(kNCDataRecoEnergy_ND)->GetBinContent(i));
02451       dataNDhist[bmind][1]->SetBinContent(i,bmtouse->GetDataPredHists().at(kCCDataRecoEnergy_ND)->GetBinContent(i));
02452     }
02453   //this is a systematically shifted beam and has no data in it. Copy from nominal beam
02454   else
02455     for(int i = 1 ; i < kNumEnergyBinsFar+1 ; i++){   
02456       bmtouse->GetDataPredHists().at(kNCDataRecoEnergy_ND)->SetBinContent(i,dataNDhist[bmind][0]->GetBinContent(i) );
02457       bmtouse->GetDataPredHists().at(kCCDataRecoEnergy_ND)->SetBinContent(i,dataNDhist[bmind][1]->GetBinContent(i) );
02458     } 
02459 }

void NCExtrapolationBeamMatrix::correctNDNCFlux ( BeamMatrix bmtouse,
int  beamIndex 
) [private]

Definition at line 2647 of file NCExtrapolationBeamMatrix.cxx.

References DoEfficencyCorrectionForFluxCorrND(), DoEfficencyCorrectionND(), DoPurityCorrectionND(), DoRecoToTrueND(), BeamMatrix::GetDataPredHists(), BeamMatrix::GetMCHists(), GetNCFromCCFluxND(), kCCDataTrueEnergyEffCorr_ND, kCCTrueEnergyEffCorr_ND, kCCTrueEnergyEffCorrFluxCorr_ND, kFluxCorrection_ND, kNCTrueEnergyEffCorrFluxCorr_ND, kNCTrueEnergyEffCorrFromCC_ND, and makeFluxCorrectionVector().

Referenced by ScaleBeamsToSameNumberOfEvents().

02647                                                    {
02648 
02649   DoPurityCorrectionND(bmtouse);
02650   DoRecoToTrueND(bmtouse, beamIndex);
02651   DoEfficencyCorrectionND(bmtouse);
02652   GetNCFromCCFluxND(bmtouse);
02653  
02654   //making correction histogram
02655   for(int i = 0; i < kNumEnergyBinsFar+1; ++i ){
02656     if(bmtouse->GetMCHists().at(kCCTrueEnergyEffCorr_ND)->GetBinContent(i) > 0)
02657       bmtouse->GetMCHists().at(kFluxCorrection_ND)->SetBinContent(i,bmtouse->GetDataPredHists().at(kCCDataTrueEnergyEffCorr_ND)->GetBinContent(i)
02658                                                                   / bmtouse->GetMCHists().at(kCCTrueEnergyEffCorr_ND)->GetBinContent(i)) ;
02659     else bmtouse->GetMCHists().at(kFluxCorrection_ND)->SetBinContent(i,0);
02660     
02661   }
02662   //applying correction
02663   for(int i = 0; i < kNumEnergyBinsFar+1; ++i ){
02664     bmtouse->GetMCHists().at(kNCTrueEnergyEffCorrFluxCorr_ND)->SetBinContent(i, bmtouse->GetMCHists().at(kNCTrueEnergyEffCorrFromCC_ND)->GetBinContent(i) 
02665                                                                              * bmtouse->GetMCHists().at(kFluxCorrection_ND)->GetBinContent(i));
02666     bmtouse->GetMCHists().at(kCCTrueEnergyEffCorrFluxCorr_ND)->SetBinContent(i, bmtouse->GetMCHists().at(kCCTrueEnergyEffCorr_ND)->GetBinContent(i) 
02667                                                                              * bmtouse->GetMCHists().at(kFluxCorrection_ND)->GetBinContent(i));
02668   }
02669   DoEfficencyCorrectionForFluxCorrND(bmtouse);
02670   makeFluxCorrectionVector(bmtouse, beamIndex);
02671 }

const Registry & NCExtrapolationBeamMatrix::DefaultConfig ( void   )  [static]

Return a default config that will be merged with the NCExtrapolationModule DefaultConfig.

Reimplemented from NCExtrapolation.

Definition at line 365 of file NCExtrapolationBeamMatrix.cxx.

References Registry::LockValues(), Registry::Set(), and Registry::UnLockValues().

00366 {
00367   static Registry r;
00368 
00369   r.UnLockValues();
00370 
00371   r.Set("FarNearSmoothingWidth", 0.00);
00372   r.Set("BeamMatrixFileLocation","/data/minos/pittam/beamdata/NewFluxHelpersRunIAllSkzpPiMinusNCUtilsBin.root");
00373   r.Set("SelectionEfficiencyFileLocation","/home/pittam/minos/devtest/EffMomCorrectionsWithNCNCUtilsBin.root");
00374 
00375   r.LockValues();
00376   return r;
00377 }

void NCExtrapolationBeamMatrix::DoBeamMatrixExtrapolation ( BeamMatrix bmtouse  )  [private]

Definition at line 1543 of file NCExtrapolationBeamMatrix.cxx.

References BeamMatrix::GetBeamMatrix(), BeamMatrix::GetDataPredHists(), BeamMatrix::GetMCHists(), kCCDataTrueEnergyFlux_ND, kCCDataTrueEnergyMatrix_FD, kCCTrueEnergyFlux_ND, and kCCTrueEnergyMatrix_FD.

Referenced by ConstructFarSpectrum().

01543                                               {
01544 
01545   int beamMatrixNearBins = bmtouse->GetBeamMatrix()->GetXaxis()->GetNbins();
01546   int beamMatrixFarBins = bmtouse->GetBeamMatrix()->GetYaxis()->GetNbins();
01547 
01548 
01549   vector<double> temp(beamMatrixNearBins+1,0) ;
01550   vector<double> tempData(beamMatrixNearBins+1,0) ;
01551   
01552   for(int i = 1; i < beamMatrixNearBins+1; ++i ){
01553     for(int j = 1; j <beamMatrixFarBins+1 ; ++j ){      
01554       temp[j] += (bmtouse->GetMCHists().at(kCCTrueEnergyFlux_ND)->GetBinContent(i) * bmtouse->GetBeamMatrix()->GetBinContent( i,j ));
01555       tempData[j] += (bmtouse->GetDataPredHists().at(kCCDataTrueEnergyFlux_ND)->GetBinContent(i) * bmtouse->GetBeamMatrix()->GetBinContent( i,j ));
01556     }
01557   }
01558   
01559   for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
01560     bmtouse->GetMCHists().at(kCCTrueEnergyMatrix_FD)->SetBinContent(i,temp[i]);
01561     bmtouse->GetDataPredHists().at(kCCDataTrueEnergyMatrix_FD)->SetBinContent(i,tempData[i]);
01562   }  
01563   temp.clear() ;
01564   tempData.clear() ;
01565 }

void NCExtrapolationBeamMatrix::DoEfficencyCorrectionFD ( BeamMatrix bmtouse  )  [private]

Definition at line 1636 of file NCExtrapolationBeamMatrix.cxx.

References BMbybeams, fBMbybeams, fluxCorrectionsFDCCinNC, BeamMatrix::GetDataPredHists(), BeamMatrix::GetEffCorrHists(), BeamMatrix::GetMCHists(), GetRunTypeFromBM(), kCCDataTrueEnergyEffCorr_FD, kCCDataTrueEnergyFlux_FD, kCCDataTrueEnergySelEffCorr_FD, kCCinNCDataTrueEnergySelEffCorr_FD, kCCinNCTrueEnergySelEffCorr_FD, kCCTrueEnergyEffCorr_FD, kCCTrueEnergyFlux_FD, kCCTrueEnergySelEffCorr_FD, BeamType::kL010z185i, kNCDataTrueEnergyEffCorr_FD, kNCDataTrueEnergyFlux_FD, kNCDataTrueEnergySelEffCorr_FD, kNCTrueEnergyBkg_FD, kNCTrueEnergyEffCorr_FD, kNCTrueEnergyFlux_FD, kNCTrueEnergySelEffCorr_FD, kRecoCCEffFDHist, kRecoNCEffFDHist, kSelCCEffFDHist, kSelCCinNCEffFDHist, and kSelNCEffFDHist.

Referenced by ConstructFarSpectrum().

01636                                             {
01637  
01638   //
01639   // Do reconstruction then selection efficiency in the Far Detector
01640   //    
01641 
01642   NC::RunUtil::ERunType rt = GetRunTypeFromBM(bmtouse) ;
01643   const NCBeam::Info beamInfoLE = NCBeam::Info(BeamType::kL010z185i, rt);
01644   int leNomBeam = fBMbybeams.find(beamInfoLE)->second ;    
01645   
01646   for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
01647 
01648     // CC correction
01649     if(bmtouse->GetEffCorrHists().at(kRecoCCEffFDHist)->GetBinContent(i) > 0){
01650       
01651       bmtouse->GetMCHists().at(kCCTrueEnergyEffCorr_FD)->SetBinContent(i,bmtouse->GetMCHists().at(kCCTrueEnergyFlux_FD)->GetBinContent(i) 
01652                                                                        * bmtouse->GetEffCorrHists().at(kRecoCCEffFDHist)->GetBinContent(i));
01653       bmtouse->GetDataPredHists().at(kCCDataTrueEnergyEffCorr_FD)->SetBinContent(i,bmtouse->GetDataPredHists().at(kCCDataTrueEnergyFlux_FD)->GetBinContent(i) 
01654                                                                                  * bmtouse->GetEffCorrHists().at(kRecoCCEffFDHist)->GetBinContent(i));
01655       
01656     }
01657     else{
01658       bmtouse->GetMCHists().at(kCCTrueEnergyEffCorr_FD)->SetBinContent(i,0);
01659       bmtouse->GetDataPredHists().at(kCCDataTrueEnergyEffCorr_FD)->SetBinContent(i,0);
01660     }
01661     
01662     // NC correction
01663     if(bmtouse->GetEffCorrHists().at(kRecoNCEffFDHist)->GetBinContent(i) > 0){
01664       bmtouse->GetMCHists().at(kNCTrueEnergyEffCorr_FD)->SetBinContent(i,bmtouse->GetMCHists().at(kNCTrueEnergyFlux_FD)->GetBinContent(i) 
01665                                                                        * bmtouse->GetEffCorrHists().at(kRecoNCEffFDHist)->GetBinContent(i));
01666       bmtouse->GetDataPredHists().at(kNCDataTrueEnergyEffCorr_FD)->SetBinContent(i,bmtouse->GetDataPredHists().at(kNCDataTrueEnergyFlux_FD)->GetBinContent(i) 
01667                                                                                  * bmtouse->GetEffCorrHists().at(kRecoNCEffFDHist)->GetBinContent(i));
01668     }
01669     else{
01670       bmtouse->GetMCHists().at(kNCTrueEnergyEffCorr_FD)->SetBinContent(i,0);
01671       bmtouse->GetDataPredHists().at(kNCDataTrueEnergyEffCorr_FD)->SetBinContent(i,0);
01672     }
01673 
01675     
01676     //CC
01677     if( bmtouse->GetEffCorrHists().at(kSelCCEffFDHist)->GetBinContent(i) > 0){
01678       bmtouse->GetMCHists().at(kCCTrueEnergySelEffCorr_FD)->SetBinContent(i,bmtouse->GetMCHists().at(kCCTrueEnergyEffCorr_FD)->GetBinContent(i) 
01679                                                                           * bmtouse->GetEffCorrHists().at(kSelCCEffFDHist)->GetBinContent(i));
01680       bmtouse->GetDataPredHists().at(kCCDataTrueEnergySelEffCorr_FD)->SetBinContent(i,bmtouse->GetDataPredHists().at(kCCDataTrueEnergyEffCorr_FD)->GetBinContent(i) 
01681                                                                                     * bmtouse->GetEffCorrHists().at(kSelCCEffFDHist)->GetBinContent(i));    
01682     }
01683     else{
01684       bmtouse->GetMCHists().at(kCCTrueEnergySelEffCorr_FD)->SetBinContent(i,0);
01685       bmtouse->GetDataPredHists().at(kCCDataTrueEnergySelEffCorr_FD)->SetBinContent(i,0);
01686     }
01687     
01688     
01689     //CC in NC
01690     if( bmtouse->GetEffCorrHists().at(kSelCCinNCEffFDHist)->GetBinContent(i) > 0){
01691       bmtouse->GetMCHists().at(kCCinNCTrueEnergySelEffCorr_FD)->SetBinContent(i,bmtouse->GetMCHists().at(kCCTrueEnergyEffCorr_FD)->GetBinContent(i) 
01692                                                                               * bmtouse->GetEffCorrHists().at(kSelCCinNCEffFDHist)->GetBinContent(i));
01693       bmtouse->GetDataPredHists().at(kCCinNCDataTrueEnergySelEffCorr_FD)->SetBinContent(i,bmtouse->GetDataPredHists().at(kCCDataTrueEnergyEffCorr_FD)->GetBinContent(i) 
01694                                                                                         * bmtouse->GetEffCorrHists().at(kSelCCinNCEffFDHist)->GetBinContent(i));    
01695     }
01696     else{
01697       bmtouse->GetMCHists().at(kCCinNCTrueEnergySelEffCorr_FD)->SetBinContent(i,0);
01698       bmtouse->GetDataPredHists().at(kCCinNCDataTrueEnergySelEffCorr_FD)->SetBinContent(i,0);
01699     }
01700     
01701     
01702     //true to reco this for CC in NC backgournd prediction
01703     fluxCorrectionsFDCCinNC.at(i-1) =  bmtouse->GetDataPredHists().at(kCCinNCDataTrueEnergySelEffCorr_FD)->GetBinContent(i)
01704       / BMbybeams[leNomBeam]->GetMCHists().at(kNCTrueEnergyBkg_FD)->GetBinContent(i);
01705     
01706     //NC
01707     if( bmtouse->GetEffCorrHists().at(kSelNCEffFDHist)->GetBinContent(i) > 0){
01708       bmtouse->GetMCHists().at(kNCTrueEnergySelEffCorr_FD)->SetBinContent(i,bmtouse->GetMCHists().at(kNCTrueEnergyEffCorr_FD)->GetBinContent(i) 
01709                                                                           * bmtouse->GetEffCorrHists().at(kSelNCEffFDHist)->GetBinContent(i));
01710       bmtouse->GetDataPredHists().at(kNCDataTrueEnergySelEffCorr_FD)->SetBinContent(i,bmtouse->GetDataPredHists().at(kNCDataTrueEnergyEffCorr_FD)->GetBinContent(i) 
01711                                                                                     * bmtouse->GetEffCorrHists().at(kSelNCEffFDHist)->GetBinContent(i));    
01712     }
01713     else{
01714       bmtouse->GetMCHists().at(kNCTrueEnergySelEffCorr_FD)->SetBinContent(i,0);
01715       bmtouse->GetDataPredHists().at(kNCDataTrueEnergySelEffCorr_FD)->SetBinContent(i,0);
01716     }    
01717   }  
01718 }

void NCExtrapolationBeamMatrix::DoEfficencyCorrectionForFluxCorrND ( BeamMatrix bmtouse  )  [private]

Definition at line 2682 of file NCExtrapolationBeamMatrix.cxx.

References BeamMatrix::GetEffCorrHists(), BeamMatrix::GetMCHists(), kCCTrueEnergyEffCorrFluxCorr_ND, kCCTrueEnergyPurCorrFluxCorr_ND, kCCTrueEnergySelEffCorrFluxCorr_ND, kNCTrueEnergyEffCorrFluxCorr_ND, kNCTrueEnergyPurCorrFluxCorr_ND, kNCTrueEnergySelEffCorrFluxCorr_ND, kRecoCCEffNDHist, kRecoNCEffNDHist, kSelCCinNCEffNDHist, and kSelNCEffNDHist.

Referenced by correctNDNCFlux().

02682                                                        {
02683  
02684   //
02685   // Do reconstruction then selection efficiency in the Near Detector to get correction due to flux
02686   //
02687   
02688   for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
02689 
02690     // NC correction
02691     if(bmtouse->GetEffCorrHists().at(kRecoNCEffNDHist)->GetBinContent(i) > 0){
02692       bmtouse->GetMCHists().at(kNCTrueEnergySelEffCorrFluxCorr_ND)->SetBinContent(i,bmtouse->GetMCHists().at(kNCTrueEnergyEffCorrFluxCorr_ND)->GetBinContent(i) 
02693                                                                                     * bmtouse->GetEffCorrHists().at(kRecoNCEffNDHist)->GetBinContent(i));
02694     }//sel
02695     else
02696       bmtouse->GetMCHists().at(kNCTrueEnergySelEffCorrFluxCorr_ND)->SetBinContent(i,0);
02697 
02698     //CC correction
02699     if(bmtouse->GetEffCorrHists().at(kRecoCCEffNDHist)->GetBinContent(i) > 0){
02700       bmtouse->GetMCHists().at(kCCTrueEnergySelEffCorrFluxCorr_ND)->SetBinContent(i,bmtouse->GetMCHists().at(kCCTrueEnergyEffCorrFluxCorr_ND)->GetBinContent(i) 
02701                                                                                   * bmtouse->GetEffCorrHists().at(kRecoCCEffNDHist)->GetBinContent(i));
02702     }//sel
02703     else
02704       bmtouse->GetMCHists().at(kCCTrueEnergySelEffCorrFluxCorr_ND)->SetBinContent(i,0);
02705 
02706         
02707     //NC
02708     if( bmtouse->GetEffCorrHists().at(kSelNCEffNDHist)->GetBinContent(i) > 0){
02709       bmtouse->GetMCHists().at(kNCTrueEnergyPurCorrFluxCorr_ND)->SetBinContent(i,bmtouse->GetMCHists().at(kNCTrueEnergySelEffCorrFluxCorr_ND)->GetBinContent(i) 
02710                                                                                        * bmtouse->GetEffCorrHists().at(kSelNCEffNDHist)->GetBinContent(i));
02711     }
02712     else
02713       bmtouse->GetMCHists().at(kNCTrueEnergyPurCorrFluxCorr_ND)->SetBinContent(i,0);
02714 
02715     //CC
02716     if( bmtouse->GetEffCorrHists().at(kSelCCinNCEffNDHist)->GetBinContent(i) > 0){
02717       bmtouse->GetMCHists().at(kCCTrueEnergyPurCorrFluxCorr_ND)->SetBinContent(i,bmtouse->GetMCHists().at(kCCTrueEnergySelEffCorrFluxCorr_ND)->GetBinContent(i) 
02718                                                                                * bmtouse->GetEffCorrHists().at(kSelCCinNCEffNDHist)->GetBinContent(i));
02719     }
02720     else
02721       bmtouse->GetMCHists().at(kCCTrueEnergyPurCorrFluxCorr_ND)->SetBinContent(i,0);
02722     
02723   }  
02724 }

void NCExtrapolationBeamMatrix::DoEfficencyCorrectionND ( BeamMatrix bmtouse  )  [private]

Definition at line 1455 of file NCExtrapolationBeamMatrix.cxx.

References BeamMatrix::GetDataPredHists(), BeamMatrix::GetEffCorrHists(), BeamMatrix::GetMCHists(), kCCDataTrueEnergyEffCorr_ND, kCCDataTrueEnergyPurCorr_ND, kCCDataTrueEnergySelEffCorr_ND, kCCTrueEnergyEffCorr_ND, kCCTrueEnergyPurCorr_ND, kCCTrueEnergySelEffCorr_ND, kRecoCCEffNDHist, and kSelCCEffNDHist.

Referenced by ConstructFarSpectrum(), and correctNDNCFlux().

01455                                             {
01456   
01457   //do selection then reconstruction in the ND
01458      
01459   for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
01460     
01461     // selection efficiency correction
01462     
01463     if(bmtouse->GetEffCorrHists().at(kSelCCEffNDHist)->GetBinContent(i) > 0){
01464       bmtouse->GetMCHists().at(kCCTrueEnergySelEffCorr_ND)->SetBinContent(i,bmtouse->GetMCHists().at(kCCTrueEnergyPurCorr_ND)->GetBinContent(i) 
01465                                                                           / bmtouse->GetEffCorrHists().at(kSelCCEffNDHist)->GetBinContent(i));
01466       bmtouse->GetDataPredHists().at(kCCDataTrueEnergySelEffCorr_ND)->SetBinContent(i,bmtouse->GetDataPredHists().at(kCCDataTrueEnergyPurCorr_ND)->GetBinContent(i) 
01467                                                                                     / bmtouse->GetEffCorrHists().at(kSelCCEffNDHist)->GetBinContent(i));
01468     }
01469     else{
01470       bmtouse->GetMCHists().at(kCCTrueEnergySelEffCorr_ND)->SetBinContent(i,0);
01471       bmtouse->GetDataPredHists().at(kCCDataTrueEnergySelEffCorr_ND)->SetBinContent(i,0);
01472     }
01473     
01474     if(bmtouse->GetEffCorrHists().at(kRecoCCEffNDHist)->GetBinContent(i) > 0){
01475       bmtouse->GetMCHists().at(kCCTrueEnergyEffCorr_ND)->SetBinContent(i,bmtouse->GetMCHists().at(kCCTrueEnergySelEffCorr_ND)->GetBinContent(i) 
01476                                                                        / bmtouse->GetEffCorrHists().at(kRecoCCEffNDHist)->GetBinContent(i));
01477       bmtouse->GetDataPredHists().at(kCCDataTrueEnergyEffCorr_ND)->SetBinContent(i,bmtouse->GetDataPredHists().at(kCCDataTrueEnergySelEffCorr_ND)->GetBinContent(i) 
01478                                                                                  / bmtouse->GetEffCorrHists().at(kRecoCCEffNDHist)->GetBinContent(i));
01479     }
01480     else{
01481       bmtouse->GetMCHists().at(kCCTrueEnergyEffCorr_ND)->SetBinContent(i,0);
01482       bmtouse->GetDataPredHists().at(kCCDataTrueEnergyEffCorr_ND)->SetBinContent(i,0);
01483     }
01484   }
01485 }

void NCExtrapolationBeamMatrix::doNDXsectionFit ( int  PrintLevel,
NCBeam beam,
int  ibm 
) [private]

Definition at line 2052 of file NCExtrapolationBeamMatrix.cxx.

References beamsFit, beamsNear, BMbybeams, bmFit, BMforFit, BMforFitTypes, corr, fBMbybeams, FillDataMCHistogramsNear(), Fit(), GetBeamMatrix(), GetEfficiencyHistograms(), NCBeam::GetInfo(), NCBeam::GetRunType(), NCType::kCC, Msg::kDebug, Msg::kInfo, BeamType::kL010z185i, NCType::kNC, kNCDataHE, kNCDataLE, kNCDataME, kNCMCHE, kNCMCLE, kNCMCME, NC::RunUtil::kRunAll, m0ScaleNC, m0ScaleNCErr, mInputScale, mScaleNC, mScaleNCErr, mScaleNeg, mScalePos, MSG, ND_XsecByParams, ND_XSectionFit, ND_XSectionFitRew, ScaleBeamsToSameNumberOfEvents(), and size.

Referenced by DoneFilling().

02052                                                               {
02053 
02054   //only need to do thisfor "All beams" beam
02055 
02056   if(beam_ndLE->GetRunType() !=  NC::RunUtil::kRunAll){ 
02057     
02058     // the way the beam vector is made means this will work
02059     NC::RunUtil::ERunType runNum = beam_ndLE->GetRunType();
02060     const NCBeam::Info beamInfo = NCBeam::Info(BeamType::kL010z185i, runNum);
02061     int nomBeam = fBMbybeams.find(beamInfo)->second ;    
02062     int allrunget = bemequ - nomBeam ;
02063     
02064     for(int i = 0 ; i < 4; i++)
02065       BMbybeams[bemequ]->SetmScaleNCparam(i,BMbybeams[allrunget]->GetmScaleNC().at(i) );   
02066     
02067     MSG("NCExtrapolationBeamMatrix",Msg::kDebug) << "Copying NDFit parameters for bemequ " << bemequ << " From " << allrunget  << endl;
02068     return;
02069   }
02070 
02071   for (Int_t i=0; i<4; i++) {
02072     m0ScaleNC.at(i) = 1.0;
02073     m0ScaleNCErr.at(i) = 0.005;
02074     mScaleNC.at(i) = 1.0;
02075     mScaleNCErr.at(i) = 0.005;
02076     mScaleNeg[i] = 0.0 ;
02077     mScalePos[i] = 0.0 ;
02078     corr[i] = 0.0 ;
02079   }
02080 
02081   //flag for if a systematic beam or nominal  
02082   bool doSyst = true ;
02083   if(bemequ == 0) doSyst = false ;
02084 
02085   // the way the beams vector is constructed mean this will work
02086   int diff = beamsNear.size()/3 ;
02087   MSG("NCExtrapolationBeamMatrix",Msg::kDebug) << "LE beamsNear.size() " << beamsNear.size() << " bemequ " << bemequ << " diff " << diff <<endl;
02088   NCBeam* beam_ndME = beamsNear.at(bemequ+diff);
02089   NCBeam* beam_ndHE = beamsNear.at(bemequ+(2*diff));
02090  
02091   MSG("NCExtrapolationBeamMatrix",Msg::kDebug) << "is LE? " << beam_ndLE->GetInfo().GetDescription() 
02092                                               << " is ME? " << beam_ndME->GetInfo().GetDescription() 
02093                                               << " is HE? " << beam_ndHE->GetInfo().GetDescription() 
02094                                               << endl ;
02095 
02096   bmFit->beamsFit.clear();
02097   bmFit->beamsFit.push_back(beam_ndLE);
02098   bmFit->beamsFit.push_back(beam_ndME);
02099   bmFit->beamsFit.push_back(beam_ndHE);
02100 
02101   BMforFit.clear();
02102   BMforFit.push_back(new BeamMatrix("beamMatrixForXsec_"+(string)beam_ndLE->GetInfo().GetDescription(),
02103                                     "titleForXsec_"+(string)beam_ndLE->GetInfo().GetDescription(),
02104                                     "ForXsec"+(string)beam_ndLE->GetInfo().GetDescription()));
02105   BMforFit.push_back(new BeamMatrix("beamMatrixForXsec_"+(string)beam_ndME->GetInfo().GetDescription(),
02106                                     "titleForXsec_"+(string)beam_ndME->GetInfo().GetDescription(),
02107                                     "ForXsec"+(string)beam_ndME->GetInfo().GetDescription()));
02108   BMforFit.push_back(new BeamMatrix("beamMatrixForXsec_"+(string)beam_ndHE->GetInfo().GetDescription(),
02109                                     "titleForXsec_"+(string)beam_ndHE->GetInfo().GetDescription(),
02110                                     "ForXsec"+(string)beam_ndHE->GetInfo().GetDescription()));
02111   
02112   // needed to write out all beams in case with systematic and normalise eff corrections properly
02113   BMforFitTypes.push_back(BMforFit) ;
02114   
02115   for(UInt_t i = 0; i < (BMforFit.size()); ++i){
02116     GetEfficiencyHistograms(bmFit->beamsFit.at(i),BMforFit.at(i));
02117     GetBeamMatrix(bmFit->beamsFit.at(i),BMforFit.at(i));
02118     FillDataMCHistogramsNear(bmFit->beamsFit.at(i), NCType::kNC, BMforFit.at(i),true);
02119     FillDataMCHistogramsNear(bmFit->beamsFit.at(i), NCType::kCC, BMforFit.at(i),true);
02120   }
02121   
02122   // may want to pass this in as part of the registry
02123   // fills the histogrmas and calculates the scales needed. 
02124 
02125   ScaleBeamsToSameNumberOfEvents(kNCDataME, doSyst);
02126 
02127   for(UInt_t i = 0; i < (bmFit->beamsFit).size(); ++i){
02128     int data = -1 ;
02129     int mc = -1;
02130 
02131     if(i == 0){data = kNCDataLE; mc = kNCMCLE; } 
02132     if(i == 1){data = kNCDataME; mc = kNCMCME; } 
02133     if(i == 2){data = kNCDataHE; mc = kNCMCHE; } 
02134         
02135     //scaling actually done here
02136     ND_XSectionFit.at(mc)->Scale(bmFit->mInputScale.at(mc));
02137     ND_XSectionFit.at(data)->Scale(bmFit->mInputScale.at(data));
02138     ND_XSectionFitRew.at(data)->Scale(bmFit->mInputScale.at(data));
02139     for(UInt_t ipar = 0 ; ipar < ND_XsecByParams.at(0).size() ; ipar++)
02140       ND_XsecByParams[i][ipar]->Scale(bmFit->mInputScale.at(mc));
02141   }
02142 
02143   Fit(PrintLevel);
02144   for(int i = 0 ; i < 4; i++)
02145     BMbybeams[bemequ]->SetmScaleNCparam(i,mScaleNC.at(i));
02146     
02147   MSG("NCExtrapolationBeamMatrix",Msg::kInfo) 
02148     << "Fitted Parameters: 0-4 " << BMbybeams[bemequ]->GetmScaleNC().at(0)  
02149     << " 4-8 " << BMbybeams[bemequ]->GetmScaleNC().at(1)  
02150     << " 8-15 " << BMbybeams[bemequ]->GetmScaleNC().at(2) 
02151     << " over 15 " << BMbybeams[bemequ]->GetmScaleNC().at(3) 
02152     << " bemequ " << bemequ 
02153     << endl;
02154   
02155 }

void NCExtrapolationBeamMatrix::DoneFilling (  )  [private, virtual]

Called after the final AddEvent() but before any calls to FindSpectraForPars().

Reimplemented from NCExtrapolation.

Definition at line 841 of file NCExtrapolationBeamMatrix.cxx.

References beamsNear, BMbybeams, det, doNDXsectionFit(), NCExtrapolation::fBeams, fBMbybeams, fFDPOTData, fFDPOTMC, FillDataMCHistogramsFar(), FillDataMCHistogramsNear(), fNDPOTData, fNDPOTMC, NCExtrapolation::fUseCC, NCExtrapolation::fUseNC, GetBeamMatrix(), NCBeam::GetBeamType(), NCBeam::Info::GetDescription(), NCBeam::GetDetector(), GetEfficiencyHistograms(), NCBeam::GetInfo(), it, NCType::kCC, kCCSig, Detector::kFar, Msg::kInfo, BeamType::kL010z185i, NCType::kNC, kNCSig, Detector::kNear, m0ScaleNC, m0ScaleNCErr, mScaleNC, mScaleNCErr, and MSG.

00842 {
00843   //only need to do nd xsec fit and fill hists once
00844   MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Done Filling" << endl;
00845 
00846   // this is cleverly taken care of in the beam construction
00847   fFDPOTData = 1.0 ;
00848   fFDPOTMC = 1.0 ;
00849   
00850   fNDPOTData = 1.0 ;
00851   fNDPOTMC =  1.0 ;
00852   
00853   int icount = 0;  
00854   int ibm = 0;
00855 
00856   for (Int_t i=0; i<4; i++) {
00857     m0ScaleNC.push_back(1.0);
00858     m0ScaleNCErr.push_back(0.005);
00859     
00860     mScaleNC.push_back(1.0);
00861     mScaleNCErr.push_back(0.005);
00862     
00863   }
00864   
00865   //this is for use in the nd xsection fitting
00866   for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
00867     const pair<NCBeam::Info, Detector::Detector_t> key = it->first;
00868     const NCBeam::Info beamInfo = key.first;
00869     const Detector::Detector_t det = key.second;    
00870     NCBeam* beam = it->second;  
00871     if(det == Detector::kNear){ 
00872       beamsNear.push_back(beam); 
00873       MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "For ND Fit Beam Description " << beamInfo.GetDescription() << " det " <<  beam->GetDetector() <<endl;
00874     }
00875   }
00876 
00877   //filling the vector of BeamMatrix Objects ibm is index for Map
00878   for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
00879     
00880     const pair<NCBeam::Info, Detector::Detector_t> key = it->first;
00881     const NCBeam::Info beamInfo = key.first;
00882     const Detector::Detector_t det = key.second;
00883     
00884     NCBeam* beam = it->second;  
00885 
00886     if(beam->GetBeamType() == BeamType::kL010z185i){    
00887       
00888       string name = (string)beam->GetInfo().GetDescription();
00889       string title = "title_"+name;
00890       string bmname = "beamMatrix_"+name;
00891       
00892       //only 1 BeamMatrix object for ND and FD
00893       if(icount%2 == 0){ 
00894         BMbybeams.push_back(new BeamMatrix(bmname,title,name));
00895         fBMbybeams[beamInfo] = ibm;
00896       }
00897       
00898       
00899       if( beam->GetDetector() == Detector::kNear){ 
00900         
00901         //get all relevant read in histograms
00902         GetEfficiencyHistograms(beam, BMbybeams.at(ibm));
00903         GetBeamMatrix(beam, BMbybeams.at(ibm));
00904         //xsection fit is done here before FillDataMCHistogramsFar
00905         doNDXsectionFit(1, beam, ibm);  
00906         
00907         if(fUseNC) FillDataMCHistogramsNear(beam, NCType::kNC, BMbybeams.at(ibm));
00908         if(fUseCC) FillDataMCHistogramsNear(beam, NCType::kCC, BMbybeams.at(ibm));      
00909 
00910         MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "For Extrapolation Filled Beam with Description " << name << " det  " << det << " ibm " << ibm 
00911                                                      << " BeamsNear.at(ibm)->GetDescription() " << beamsNear.at(ibm)->GetInfo().GetDescription() << endl;
00912 
00913       }
00914       
00915       if( beam->GetDetector() == Detector::kFar){
00916         
00917         BMbybeams[ibm]->GetFDDataHists().at(kNCSig)->Reset("ICE");
00918         BMbybeams[ibm]->GetFDDataHists().at(kCCSig)->Reset("ICE");
00919         if(fUseNC) FillDataMCHistogramsFar(beam, NCType::kNC, BMbybeams.at(ibm));
00920         if(fUseCC) FillDataMCHistogramsFar(beam, NCType::kCC, BMbybeams.at(ibm));
00921         
00922         MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "For Extrapolation Filled Beam with Description " << name << " det  " << det << " ibm " << ibm <<endl ;   
00923         MSG("NCExtrapolationBeamMatrix",Msg::kInfo) 
00924           << "Fitted Parameters: 0-4 " << BMbybeams[ibm]->GetmScaleNC().at(0)  
00925           << " 4-8 " << BMbybeams[ibm]->GetmScaleNC().at(1)  
00926           << " 8-15 " << BMbybeams[ibm]->GetmScaleNC().at(2) 
00927           << " over 15 " << BMbybeams[ibm]->GetmScaleNC().at(3) 
00928           << " size " << BMbybeams[ibm]->GetmScaleNC().size() 
00929           << endl;
00930       }
00931     }
00932     if(icount%2 == 1) ibm++; 
00933     icount++ ;
00934   }
00935 }

void NCExtrapolationBeamMatrix::DoOscillations ( const NC::OscProb::OscPars pars,
BeamMatrix bmtouse 
) [private]

Definition at line 1821 of file NCExtrapolationBeamMatrix.cxx.

References BMbybeams, NC::Utility::CloneFast(), fBMbybeams, fluxCorrectionsFDCCinNC, fNominalBM, fOscBM, BeamMatrix::Get2DHists(), BeamMatrix::Get2DOscHists(), BeamMatrix::GetDataPredHists(), BeamMatrix::GetFitHists(), BeamMatrix::GetmScaleNC(), GetRunTypeFromBM(), NC::OscProb::OscPars::IsEquiv(), NCType::kBaseLineFar, NCType::kCC, kCCBeamNue, kCCBkg, kCCDataRecoEnergySig_FD, kCCDataRecoEnergyWrongSignBkg_FD, kCCDataTrueEnergySelEffCorr_FD, kCCOscNue, kCCSig, kCCTau, kCCTrueEnergy_FD, kCCTrueRecoSig_FD, kCCWrongSignBkg, Msg::kError, BeamType::kL010z185i, NCType::kNC, kNCBeamNue, kNCBkg, kNCDataRecoEnergySig_FD, kNCDataRecoEnergyWrongSignBkg_FD, kNCDataTrueEnergySelEffCorr_FD, kNCOscNue, kNCSig, kNCTau, kNCTrueEnergy_FD, kNCTrueRecoSig_FD, kNCWrongSignBkg, NCType::kNuEToNuE, kNumTrueEnergyBins, NCType::kNuMuToNuE, NCType::kNuMuToNuMu, NCType::kNuMuToNuTau, MSG, normaliseRecoToTrue(), rebinRecoTrueToRecoRecoHist(), t2rCC, t2rCCRebin, t2rNC, t2rNCRebin, NC::OscProb::OscPars::TransitionProbability(), true_bins, and WhichFitParam().

Referenced by ConstructFarSpectrum().

01821                                                                  {
01822 
01823   //means a vector of arrays!
01824   static vector<double*> survivalProbability;
01825   
01826   //means a basic vector of doubles!
01827   vector<double> OscBeamMatrixCC(kNumEnergyBinsFar+1, 0.);  
01828   vector<double> OscBeamMatrixNC(kNumEnergyBinsFar+1, 0.);  
01829   
01830   static bool first = true;
01831   if(first){
01832     first = false;
01833     
01834     fNominalBM.resize(bmtouse->Get2DOscHists().size());
01835     fOscBM.resize(bmtouse->Get2DOscHists().size());
01836     survivalProbability.resize(bmtouse->Get2DOscHists().size());    
01837     
01838     for(unsigned int h = 0; h < bmtouse->Get2DOscHists().size(); ++h){
01839       fNominalBM[h].resize(kNumEnergyBinsFar);
01840       fOscBM[h].resize(kNumEnergyBinsFar);
01841       survivalProbability[h] = new double[kNumTrueEnergyBins];
01842     }
01843   }
01844   
01845   
01846   static const NC::OscProb::OscPars* prevpars = 0;
01847   
01848   // Transition probability does not depend on reco energy - so we
01849   // can precalculate it out here.
01850   // Reuse the calculations from last time round if the oscillations
01851   // parameters haven't changed since then.
01852   if(!prevpars || !pars->IsEquiv(prevpars)){
01853     for(int h = 0; h < int(bmtouse->Get2DOscHists().size()); ++h){
01854       NCType::EOscMode oscType;
01855       NCType::EEventType interaction;
01856       
01857       switch(h){
01858       case kNCSig:
01859       case kCCBkg:
01860       case kNCWrongSignBkg:
01861         interaction = NCType::kNC;
01862         oscType     = NCType::kNuMuToNuMu;
01863         break;
01864       case kCCSig:
01865       case kNCBkg:
01866       case kCCWrongSignBkg:
01867         interaction = NCType::kCC;
01868         oscType     = NCType::kNuMuToNuMu;
01869         break;
01870       case kNCTau:
01871       case kCCTau:
01872         interaction = NCType::kCC;
01873         oscType     = NCType::kNuMuToNuTau;
01874         break;
01875       case kNCBeamNue:
01876       case kCCBeamNue:
01877         interaction = NCType::kCC;
01878         oscType     = NCType::kNuEToNuE;
01879         break;
01880       case kNCOscNue:
01881       case kCCOscNue:
01882         interaction = NCType::kCC;
01883         oscType     = NCType::kNuMuToNuE;
01884         break;
01885       default:
01886         MSG("NCExtrapolationBeamMatrix", Msg::kError) << "Interaction type or flavor is incorrect" << endl;
01887         return;
01888       }
01889 
01890       double* sph = survivalProbability[h];
01891       for(int j = 0; j < kNumTrueEnergyBins; ++j){
01892         sph[j] = pars->TransitionProbability(oscType,
01893                                              interaction,
01894                                              NCType::kBaseLineFar,
01895                                              (true_bins[j+1] + true_bins[j]) /2.0 );
01896       } // end for j
01897     } // end for h
01898   } // end if
01899   
01900   delete prevpars;
01901   // Sadly OscPars::Clone isn't labelled const for technical reasons.
01902   // Need to cast away constness so that we can call it.
01903   prevpars = const_cast<NC::OscProb::OscPars*>(pars)->Clone();
01904   
01905   
01906   //This is the original simple version
01907   /*
01908     for(int i = 0; i < kNumEnergyBinsFar; ++i ){
01909     for(int h = 0; h < int(bmtouse->Get2DOscHists().size()); ++h){
01910     fNominalBM[h][i] = 0.;
01911     fOscBM[h][i] = 0.;
01912     for(int j = 0; j < kNumTrueEnergyBins; ++j){
01913         //double xsecFitParam = mScaleNC.at(WhichFitParam(j+1)) ;
01914         double xsecFitParam = 1.0 ;
01915         double jbineq = (true_bins[j+1] + true_bins[j]) /2.0 ; //(j+1)*0.1 ;// make binwidth
01916         //if(jbineq > 120.0) jbineq = 120.0 ;
01917         //if(h == kNCSig ) MSG("NCExtrapolationBeamMatrix",Msg::kInfo) << "PURITAN BIN j " << j << " jbineq " << jbineq << " mScaleNC.at(WhichFitParam(jbineq)) " << mScaleNC.at(WhichFitParam(jbineq)) << endl;
01918         if(h == kNCSig || h == kNCWrongSignBkg || h == kCCBkg) xsecFitParam = mScaleNC.at(WhichFitParam(jbineq)) ;  
01919         const double tmpd = bmtouse->Get2DOscHists()[h]->GetCellContent( i+1, j+1 );
01920         fNominalBM[h][i] += tmpd; // *xsecFitParam;
01921         fOscBM[h][i] += survivalProbability[j][h]*tmpd*xsecFitParam;;
01922       }//end loop over true bins
01923       }//end loop over histograms
01924     }
01925   */
01926   // And this is the optimized version
01927   const int nx = kNumEnergyBinsFar+2;
01928   for(unsigned int h = 0; h < bmtouse->Get2DOscHists().size(); ++h){
01929     const double* t2ra =  bmtouse->Get2DOscHists()[h]->fArray;
01930     vector<double>& nh = fNominalBM[h];
01931     vector<double>& oh = fOscBM[h];
01932     for(int i = 0; i < kNumEnergyBinsFar; ++i){
01933       double& nhi = nh[i];
01934       double& ohi = oh[i];
01935       nhi = 0;
01936       ohi = 0;
01937       const double* sph = survivalProbability[h];
01938       const int oset = i+1+nx;
01939       for(int j = 0; j < kNumTrueEnergyBins; ++j){
01940         double xsecFitParam = 1.0 ;
01941         double jbineq = (true_bins[j+1] + true_bins[j]) /2.0 ;
01942         //using the nd xsecweights on the Wrong sign bkg and the nominal NC should you wish to study MC only
01943         if(h == kNCSig || h == kNCWrongSignBkg) xsecFitParam = bmtouse->GetmScaleNC().at(WhichFitParam(jbineq)) ;  
01944         //if(h == kNCSig || h == kNCWrongSignBkg || h == kCCBkg) xsecFitParam = bmtouse->GetmScaleNC().at(WhichFitParam(jbineq)) ;  
01945         if(h == kNCBkg) xsecFitParam = fluxCorrectionsFDCCinNC.at(i) ;
01946         const double tmpd = t2ra[nx*j+oset];
01947         nhi += tmpd;
01948         ohi += sph[j]*tmpd*xsecFitParam;
01949       }//end loop over true bins 
01950     } // end for i
01951   }//end loop over histograms
01952 
01953 
01954   //need to truetoreco the BM prediction.
01955   //4 step plan 
01956   //1. multiply oscs into true  
01957   //2. rebin 2d histogram to be farBinning by farBinning
01958   //3. normalise new 2d histogram
01959   //4. true to reco
01960   
01961   vector<double> tempOscCC(kNumEnergyBinsFar+1,0) ;
01962   vector<double> tempOscNC(kNumEnergyBinsFar+1,0) ;
01963 
01964   //need to use nominal for normalisation or Systematics cancel at true to reco stage
01965   NC::RunUtil::ERunType rt = GetRunTypeFromBM(bmtouse) ;
01966   const NCBeam::Info beamInfoLE = NCBeam::Info(BeamType::kL010z185i, rt);
01967   int leNomBeam = fBMbybeams.find(beamInfoLE)->second ;    
01968 
01969   t2rCCRebin->Reset("ICE");
01970   t2rNCRebin->Reset("ICE");
01971 
01972   t2rCC = (TH2D*)CloneFast(bmtouse->Get2DHists()[kCCTrueRecoSig_FD]);
01973   t2rNC = (TH2D*)CloneFast(bmtouse->Get2DHists()[kNCTrueRecoSig_FD]);
01974   
01975   //step 1
01976   for(int i = 1; i < kNumTrueEnergyBins+1; ++i ){//true
01977     for(int j = 1; j < kNumEnergyBinsFar+1; ++j ){//reco
01978       const double spCC = survivalProbability[kCCSig][i-1];
01979       const double spNC = survivalProbability[kNCSig][i-1];
01980       t2rCC->SetBinContent(i,j, spCC * t2rCC->GetBinContent(i,j) );
01981       t2rNC->SetBinContent(i,j, spNC * t2rNC->GetBinContent(i,j) );
01982     }
01983   }
01984   //step 2
01985 
01986   rebinRecoTrueToRecoRecoHist(t2rCC, t2rCCRebin); 
01987   rebinRecoTrueToRecoRecoHist(t2rNC, t2rNCRebin); 
01988 
01989   //step 3
01990   normaliseRecoToTrue(t2rCCRebin,BMbybeams[leNomBeam]->GetMCHists().at( kCCTrueEnergy_FD ) );
01991   normaliseRecoToTrue(t2rNCRebin,BMbybeams[leNomBeam]->GetMCHists().at( kNCTrueEnergy_FD ) );
01992 
01993   //step 4
01994   for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){//true
01995     for(int j = 1; j < kNumEnergyBinsFar+1; ++j ){//reco
01996       //OscBeamMatrixCC[j-1] += ( bmtouse->GetMCHists().at(kCCTrueEnergy_FD)->GetBinContent(i) * t2rCCRebin->GetBinContent(i,j)) ; 
01997       //OscBeamMatrixNC[j-1] += ( bmtouse->GetMCHists().at(kNCTrueEnergy_FD)->GetBinContent(i) * t2rNCRebin->GetBinContent(i,j)) ; 
01998       
01999       OscBeamMatrixCC[j-1] += ( bmtouse->GetDataPredHists().at(kCCDataTrueEnergySelEffCorr_FD)->GetBinContent(i) * t2rCCRebin->GetBinContent(i,j)) ; 
02000       OscBeamMatrixNC[j-1] += ( bmtouse->GetDataPredHists().at(kNCDataTrueEnergySelEffCorr_FD)->GetBinContent(i) * t2rNCRebin->GetBinContent(i,j)) ; 
02001       
02002     }
02003   }
02004 
02005   bmtouse->GetDataPredHists().at(kNCDataRecoEnergySig_FD)->Reset("ICE");
02006   bmtouse->GetDataPredHists().at(kNCDataRecoEnergyWrongSignBkg_FD)->Reset("ICE");
02007   
02008   bmtouse->GetDataPredHists().at(kCCDataRecoEnergySig_FD)->Reset("ICE");
02009   bmtouse->GetDataPredHists().at(kCCDataRecoEnergyWrongSignBkg_FD)->Reset("ICE");
02010 
02011   
02012   for(int i = 0; i < kNumEnergyBinsFar; ++i ){
02013     
02014     // might want the MC as well at some point
02015     
02016     bmtouse->GetDataPredHists().at(kNCDataRecoEnergySig_FD)->SetBinContent(i+1,OscBeamMatrixNC[i] );
02017     bmtouse->GetDataPredHists().at(kNCDataRecoEnergyWrongSignBkg_FD)->SetBinContent(i+1,fOscBM[kNCWrongSignBkg][i] );
02018     
02019     bmtouse->GetDataPredHists().at(kCCDataRecoEnergySig_FD)->SetBinContent(i+1,OscBeamMatrixCC[i]);
02020     bmtouse->GetDataPredHists().at(kCCDataRecoEnergyWrongSignBkg_FD)->SetBinContent(i+1,fOscBM[kCCWrongSignBkg][i]);
02021 
02022 
02023     const double totalMCNCBM_fd = OscBeamMatrixNC[i] //fOscBM[kNCSig][i] //OscBeamMatrixNC[i]      
02024       + fOscBM[kNCWrongSignBkg][i] 
02025       + fOscBM[kNCBkg][i]
02026       + fOscBM[kNCTau][i]
02027       + fOscBM[kNCBeamNue][i]
02028       + fOscBM[kNCOscNue][i];
02029     
02030     const double totalMCCCBM_fd = OscBeamMatrixCC[i] //fOscBM[kCCSig][i] //OscBeamMatrixCC[i]    
02031       + fOscBM[kCCWrongSignBkg][i] 
02032       + fOscBM[kCCBkg][i]
02033       + fOscBM[kCCTau][i]
02034       + fOscBM[kCCBeamNue][i]
02035       + fOscBM[kCCOscNue][i];
02036     
02037     
02038     bmtouse->GetFitHists()[kNCSig]->SetBinContent(i+1,totalMCNCBM_fd);
02039     bmtouse->GetFitHists()[kCCSig]->SetBinContent(i+1,totalMCCCBM_fd);
02040     
02041   }
02042 
02043   delete t2rNC ;
02044   delete t2rCC ;
02045   
02046   OscBeamMatrixNC.clear();
02047   OscBeamMatrixCC.clear();
02048 }

void NCExtrapolationBeamMatrix::DoPurityCorrectionND ( BeamMatrix bmtouse  )  [private]
void NCExtrapolationBeamMatrix::DoRecoToTrueND ( BeamMatrix bmtouse,
int  beamindex = -1 
) [private]

Definition at line 1292 of file NCExtrapolationBeamMatrix.cxx.

References BMforFit, NC::Utility::CloneFast(), BeamMatrix::Get2DHists(), BeamMatrix::GetDataPredHists(), BeamMatrix::GetMCHists(), kCCDataRecoEnergyPurCorr_ND, kCCDataTrueEnergyPurCorr_ND, kCCRecoEnergyPurCorr_ND, kCCRecoTrueSig_ND, kCCTrueEnergyPurCorr_ND, and t2rCCnd.

Referenced by ConstructFarSpectrum(), and correctNDNCFlux().

01292                                                   {
01293 
01294 
01295   vector<double> tempMC(kNumEnergyBinsFar+1,0) ;
01296   vector<double> tempData(kNumEnergyBinsFar+1,0) ;
01297   
01298   t2rCCnd = (TH2D*)CloneFast(bmtouse->Get2DHists()[kCCRecoTrueSig_ND]);
01299   
01300   double bincontent = 0;  
01301   if(beamindex < 0){
01302     for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){//reco
01303       for(int j = 1; j < kNumEnergyBinsFar+1; ++j ){//true
01304         // I think I want to use the ndmc to normalise to get cancellations 
01305         // Normalise reco to true matrix
01306         if( bmtouse->GetMCHists().at(kCCRecoEnergyPurCorr_ND)->GetBinContent(i) > 0){
01307           bincontent = (t2rCCnd->GetBinContent(j,i) / bmtouse->GetMCHists().at(kCCRecoEnergyPurCorr_ND)->GetBinContent(i));
01308         t2rCCnd->SetBinContent(j,i, bincontent);
01309         }
01310       else t2rCCnd->SetBinContent(j,i,0);
01311       }
01312     }
01313   }
01314   else{
01315     for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){//reco
01316       for(int j = 1; j < kNumEnergyBinsFar+1; ++j ){//true
01317         // I think I want to use the ndmc to normalise to get cancellations 
01318         // Normalise reco to true matrix
01319         if( bmtouse->GetMCHists().at(kCCRecoEnergyPurCorr_ND)->GetBinContent(i) > 0){
01320           bincontent = (t2rCCnd->GetBinContent(j,i) /  BMforFit.at(beamindex)->GetMCHists().at(kCCRecoEnergyPurCorr_ND)->GetBinContent(i));
01321           t2rCCnd->SetBinContent(j,i, bincontent);
01322         }
01323         else t2rCCnd->SetBinContent(j,i,0);
01324       }
01325     }
01326   }
01327 
01328   for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){//reco
01329     for(int j = 1; j < kNumEnergyBinsFar+1; ++j ){//true
01330       // correct reco to true
01331       tempMC[j-1] += (bmtouse->GetMCHists().at(kCCRecoEnergyPurCorr_ND)->GetBinContent(i) * t2rCCnd->GetBinContent(j,i) );
01332       tempData[j-1] += (bmtouse->GetDataPredHists().at(kCCDataRecoEnergyPurCorr_ND)->GetBinContent(i) * t2rCCnd->GetBinContent(j,i) );
01333     }
01334   }
01335   
01336   for(int i=1;i<=kNumEnergyBinsFar;i++){
01337     bmtouse->GetMCHists().at(kCCTrueEnergyPurCorr_ND)->SetBinContent(i,tempMC[i-1]);
01338     bmtouse->GetDataPredHists().at(kCCDataTrueEnergyPurCorr_ND)->SetBinContent(i,tempData[i-1]);
01339     
01340   }
01341 
01342   delete t2rCCnd ;
01343   tempMC.clear() ;
01344   tempData.clear() ;
01345 
01346   return;  
01347 }

void NCExtrapolationBeamMatrix::DoUnoscPurityCorrectionFD ( BeamMatrix bmtouse  )  [private]

Definition at line 1790 of file NCExtrapolationBeamMatrix.cxx.

References BeamMatrix::GetDataPredHists(), BeamMatrix::GetMCHists(), kCCDataRecoEnergyPred_FD, kCCDataRecoEnergyPurCorr_FD, kCCPurity_FD, kCCRecoEnergyPred_FD, kCCRecoEnergyPurCorr_FD, kNCDataRecoEnergyPred_FD, kNCDataRecoEnergyPurCorr_FD, kNCPurity_FD, kNCRecoEnergyPred_FD, and kNCRecoEnergyPurCorr_FD.

Referenced by ConstructFarSpectrum().

01790                                               {
01791  
01792   for(int i = 0; i < kNumEnergyBinsFar+1; ++i ){
01793     //CC
01794     if(bmtouse->GetMCHists().at(kCCPurity_FD)->GetBinContent(i) > 0 ){
01795       bmtouse->GetMCHists().at(kCCRecoEnergyPred_FD)->SetBinContent(i, bmtouse->GetMCHists().at(kCCRecoEnergyPurCorr_FD)->GetBinContent(i)
01796                                                                     / bmtouse->GetMCHists().at(kCCPurity_FD)->GetBinContent(i));
01797       bmtouse->GetDataPredHists().at(kCCDataRecoEnergyPred_FD)->SetBinContent(i, bmtouse->GetDataPredHists().at(kCCDataRecoEnergyPurCorr_FD)->GetBinContent(i)
01798                                                                               / bmtouse->GetMCHists().at(kCCPurity_FD)->GetBinContent(i));
01799     }
01800     else{
01801       bmtouse->GetMCHists().at(kCCRecoEnergyPred_FD)->SetBinContent(i, 0);
01802       bmtouse->GetDataPredHists().at(kCCDataRecoEnergyPred_FD)->SetBinContent(i, 0);
01803     }  
01804   
01805     //NC
01806     if(bmtouse->GetMCHists().at(kNCPurity_FD)->GetBinContent(i) > 0 ){
01807       bmtouse->GetMCHists().at(kNCRecoEnergyPred_FD)->SetBinContent(i, bmtouse->GetMCHists().at(kNCRecoEnergyPurCorr_FD)->GetBinContent(i)
01808                                                                     / bmtouse->GetMCHists().at(kNCPurity_FD)->GetBinContent(i));
01809       bmtouse->GetDataPredHists().at(kNCDataRecoEnergyPred_FD)->SetBinContent(i, bmtouse->GetDataPredHists().at(kNCDataRecoEnergyPurCorr_FD)->GetBinContent(i)
01810                                                                               / bmtouse->GetMCHists().at(kNCPurity_FD)->GetBinContent(i));
01811     }
01812     else{
01813       bmtouse->GetMCHists().at(kNCRecoEnergyPred_FD)->SetBinContent(i, 0);
01814       bmtouse->GetDataPredHists().at(kNCDataRecoEnergyPred_FD)->SetBinContent(i, 0);
01815     }  
01816   }
01817   return;  
01818 }

void NCExtrapolationBeamMatrix::DoUnoscTrueToRecoFD ( BeamMatrix bmtouse  )  [private]

Definition at line 1721 of file NCExtrapolationBeamMatrix.cxx.

References BMbybeams, NC::Utility::CloneFast(), fBMbybeams, BeamMatrix::Get2DHists(), BeamMatrix::GetDataPredHists(), NCBeam::Info::GetDescription(), BeamMatrix::GetMCHists(), GetRunTypeFromBM(), kCCDataRecoEnergyPurCorr_FD, kCCDataTrueEnergySelEffCorr_FD, kCCRecoEnergyPurCorr_FD, kCCTrueEnergy_FD, kCCTrueEnergySelEffCorr_FD, kCCTrueRecoSig_FD, Msg::kDebug, BeamType::kL010z185i, kNCDataRecoEnergyPurCorr_FD, kNCDataTrueEnergySelEffCorr_FD, kNCRecoEnergyPurCorr_FD, kNCTrueEnergy_FD, kNCTrueEnergySelEffCorr_FD, kNCTrueRecoSig_FD, MSG, normaliseRecoToTrue(), rebinRecoTrueToRecoRecoHist(), t2rCC, t2rCCRebin, t2rNC, and t2rNCRebin.

Referenced by ConstructFarSpectrum().

01721                                         {
01722   
01723   
01724   vector<double> tempCCMC(kNumEnergyBinsFar+1,0) ;
01725   vector<double> tempNCMC(kNumEnergyBinsFar+1,0) ;
01726   vector<double> tempCCData(kNumEnergyBinsFar+1,0) ;
01727   vector<double> tempNCData(kNumEnergyBinsFar+1,0) ;
01728 
01729   //need to use nominal for normalisation or Systematics cancel at true to reco stage
01730   
01731   NC::RunUtil::ERunType rt = GetRunTypeFromBM(bmtouse) ;
01732   const NCBeam::Info beamInfoLE = NCBeam::Info(BeamType::kL010z185i, rt);
01733   int leNomBeam = fBMbybeams.find(beamInfoLE)->second ;    
01734     
01735   MSG("NCExtrapolationBeamMatrix", Msg::kDebug) << "Normalising TrueToOsc using " << beamInfoLE.GetDescription() << " leNomBeam " << leNomBeam << endl;  
01736 
01737   //should be ok to use these here
01738   t2rCCRebin->Reset("ICE");
01739   t2rNCRebin->Reset("ICE");
01740   
01741   t2rCC = (TH2D*)CloneFast(bmtouse->Get2DHists()[kCCTrueRecoSig_FD]);
01742   t2rNC = (TH2D*)CloneFast(bmtouse->Get2DHists()[kNCTrueRecoSig_FD]);
01743   
01744   //see DoOscillations for the 4 step plan
01745   //step 1
01746   //unosc so no need for this step
01747   //step 2
01748 
01749   rebinRecoTrueToRecoRecoHist(t2rCC, t2rCCRebin); 
01750   rebinRecoTrueToRecoRecoHist(t2rNC, t2rNCRebin); 
01751 
01752   //step 3
01753   normaliseRecoToTrue(t2rCCRebin,BMbybeams[leNomBeam]->GetMCHists().at( kCCTrueEnergy_FD ) );
01754   normaliseRecoToTrue(t2rNCRebin,BMbybeams[leNomBeam]->GetMCHists().at( kNCTrueEnergy_FD ) );
01755 
01756   //step 4
01757   for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){//true
01758     for(int j = 1; j < kNumEnergyBinsFar+1; ++j ){//reco
01759       // correct reco to true CC
01760       tempCCMC[j-1] += (bmtouse->GetMCHists().at(kCCTrueEnergySelEffCorr_FD)->GetBinContent(i) * t2rCCRebin->GetBinContent(i,j) );
01761       tempCCData[j-1] += (bmtouse->GetDataPredHists().at(kCCDataTrueEnergySelEffCorr_FD)->GetBinContent(i) * t2rCCRebin->GetBinContent(i,j));
01762       
01763       // correct reco to true NC
01764       tempNCMC[j-1] += (bmtouse->GetMCHists().at(kNCTrueEnergySelEffCorr_FD)->GetBinContent(i) * t2rNCRebin->GetBinContent(i,j) );
01765       tempNCData[j-1] += (bmtouse->GetDataPredHists().at(kNCDataTrueEnergySelEffCorr_FD)->GetBinContent(i) * t2rNCRebin->GetBinContent(i,j) );
01766     }
01767   }
01768   
01769   for(int i = 0; i < kNumEnergyBinsFar; ++i ){//reco
01770     bmtouse->GetMCHists().at(kCCRecoEnergyPurCorr_FD)->SetBinContent(i+1,tempCCMC[i]);
01771     bmtouse->GetDataPredHists().at(kCCDataRecoEnergyPurCorr_FD)->SetBinContent(i+1,tempCCData[i]);
01772     
01773     // correct reco to true NC
01774     bmtouse->GetMCHists().at(kNCRecoEnergyPurCorr_FD)->SetBinContent(i+1,tempNCMC[i]);
01775     bmtouse->GetDataPredHists().at(kNCDataRecoEnergyPurCorr_FD)->SetBinContent(i+1,tempNCData[i]);
01776   }
01777   
01778   tempCCMC.clear() ;
01779   tempCCData.clear() ;
01780   tempNCMC.clear() ;
01781   tempNCData.clear() ;
01782   
01783   delete t2rCC;
01784   delete t2rNC;
01785   
01786   return;  
01787 }

void NCExtrapolationBeamMatrix::DoXsectionCorrectionFD ( BeamMatrix bmtouse  )  [private]

Definition at line 1607 of file NCExtrapolationBeamMatrix.cxx.

References fFDFidVolMass, fFDPOTData, fFDPOTMC, fFluxPOT, BeamMatrix::GetDataPredHists(), BeamMatrix::GetMCHists(), kCCDataTrueEnergyFlux_FD, kCCDataTrueEnergyMatrix_FD, kCCTrueEnergyFlux_FD, kCCTrueEnergyMatrix_FD, Msg::kInfo, kNCDataTrueEnergyFlux_FD, kNCDataTrueEnergyMatrix_FD, kNCTrueEnergyFlux_FD, kNCTrueEnergyMatrix_FD, and MAXMSG.

Referenced by ConstructFarSpectrum().

01607                                            {
01608 
01609   // Doesn't actually do Xsections as they are folded into the BeamMatrix
01610   // it does do the fiducial mass and pot correction
01611   
01612   MAXMSG("NCExtrapolationBeamMatrix", Msg::kInfo, 1) << "FD Fid Vol " << fFDFidVolMass 
01613                                                      << " FD POT MC " << fFDPOTMC
01614                                                      << " FD POT Data " << fFDPOTData
01615                                                      << endl;  
01616   
01617 
01618   for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
01619     
01620     //CC
01621     bmtouse->GetMCHists().at(kCCTrueEnergyFlux_FD)->SetBinContent(i,bmtouse->GetMCHists().at(kCCTrueEnergyMatrix_FD)->GetBinContent(i)
01622                                                                   *(fFDFidVolMass*fFDPOTMC)/fFluxPOT);
01623     bmtouse->GetDataPredHists().at(kCCDataTrueEnergyFlux_FD)->SetBinContent(i,bmtouse->GetDataPredHists().at(kCCDataTrueEnergyMatrix_FD)->GetBinContent(i)
01624                                                                             *(fFDFidVolMass*fFDPOTData));
01625     
01626     //NC
01627     bmtouse->GetMCHists().at(kNCTrueEnergyFlux_FD)->SetBinContent(i,bmtouse->GetMCHists().at(kNCTrueEnergyMatrix_FD)->GetBinContent(i)
01628                                                                   *(fFDFidVolMass*fFDPOTMC)/fFluxPOT);
01629     bmtouse->GetDataPredHists().at(kNCDataTrueEnergyFlux_FD)->SetBinContent(i,bmtouse->GetDataPredHists().at(kNCDataTrueEnergyMatrix_FD)->GetBinContent(i)
01630                                                                             *(fFDFidVolMass*fFDPOTData));
01631     
01632   }
01633 }

void NCExtrapolationBeamMatrix::DoXsectionCorrectionND ( BeamMatrix bmtouse  )  [private]

Definition at line 1487 of file NCExtrapolationBeamMatrix.cxx.

References fFluxPOT, fNDFidVolMass, fNDPOTData, fNDPOTMC, BeamMatrix::GetDataPredHists(), BeamMatrix::GetMCHists(), kCCDataTrueEnergyEffCorr_ND, kCCDataTrueEnergyFlux_ND, kCCTrueEnergyEffCorr_ND, kCCTrueEnergyFlux_ND, Msg::kInfo, and MAXMSG.

Referenced by ConstructFarSpectrum().

01487                                            {
01488   
01489   // Doesn't actually do Xsections as they are folded into the BeamMatrix
01490   // it does do the fiducial mass and pot correction
01491 
01492   MAXMSG("NCExtrapolationBeamMatrix", Msg::kInfo, 1) << "ND Fid Vol " << fNDFidVolMass 
01493                                                << " ND POT MC " << fNDPOTMC
01494                                                << " ND POT Data " << fNDPOTData
01495                                                << endl;
01496 
01497   
01498   for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
01499     bmtouse->GetMCHists().at(kCCTrueEnergyFlux_ND)->SetBinContent(i,bmtouse->GetMCHists().at(kCCTrueEnergyEffCorr_ND)->GetBinContent(i)*fFluxPOT/(fNDFidVolMass*fNDPOTMC));
01500     bmtouse->GetDataPredHists().at(kCCDataTrueEnergyFlux_ND)->SetBinContent(i,bmtouse->GetDataPredHists().at(kCCDataTrueEnergyEffCorr_ND)->GetBinContent(i)/(fNDFidVolMass*fNDPOTData));
01501   }
01502 }

bool NCExtrapolationBeamMatrix::EnableNearToFarExtrapolation ( bool  enable  )  [virtual]

Enable or disable corrections to FD spectra based on ND spectra.

Parameters:
enable true to perform extrapolation. false to give produce unextrapolated spectra
Returns:
The previous value of this setting

The default value of this setting should be to true

Reimplemented from NCExtrapolation.

Definition at line 1239 of file NCExtrapolationBeamMatrix.cxx.

References fDoExtrapolation.

01240 {
01241   const bool ret = fDoExtrapolation;
01242   fDoExtrapolation = enable;
01243   return ret;
01244 }

void NCExtrapolationBeamMatrix::FillDataMCHistogramsFar ( NCBeam beam,
NCType::EEventType  nccc,
BeamMatrix bmtouse 
) [private]

Definition at line 681 of file NCExtrapolationBeamMatrix.cxx.

References MuELoss::e, FillHistogramSmoothed(), FillHistogramSmoothed2D(), BeamMatrix::Get2DHists(), BeamMatrix::Get2DOscHists(), NCEnergyBin::GetDataInformation(), NCEnergyBin::GetDataVectorSize(), NCBeam::GetEnergyBin(), BeamMatrix::GetFDDataHists(), NCEnergyBin::GetMCBackgroundInformation(), NCEnergyBin::GetMCBackgroundVectorSize(), NCEnergyBin::GetMCBeamNuEInformation(), NCEnergyBin::GetMCBeamNuEVectorSize(), BeamMatrix::GetMCHists(), NCEnergyBin::GetMCInformation(), NCEnergyBin::GetMCNuTauInformation(), NCEnergyBin::GetMCNuTauVectorSize(), NCEnergyBin::GetMCOscNuEInformation(), NCEnergyBin::GetMCOscNuEVectorSize(), NCEnergyBin::GetMCSignalVectorSize(), BeamMatrix::GetmScaleNC(), NCBeam::GetNumberEnergyBins(), NCType::kCC, kCCBeamNue, kCCBkg, kCCOscNue, kCCPurity_FD, kCCRecoBkg_FD, kCCRecoEnergy_FD, kCCRecoEnergyBeamNueBkg_FD, kCCSelRecoEnergy_FD, kCCSig, kCCTau, kCCTrueEnergy_FD, kCCTrueEnergyBkg_FD, kCCTrueEnergyXsecFit_FD, kCCTrueRecoBkg_FD, kCCTrueRecoSig_FD, kCCWrongSignBkg, Detector::kFar, NCType::kNC, kNCBeamNue, kNCBkg, kNCOscNue, kNCPurity_FD, kNCRecoBkg_FD, kNCRecoEnergy_FD, kNCRecoEnergyBeamNueBkg_FD, kNCSelRecoEnergy_FD, kNCSig, kNCTau, kNCTrueEnergy_FD, kNCTrueEnergyBkg_FD, kNCTrueEnergyXsecFit_FD, kNCTrueRecoBkg_FD, kNCTrueRecoSig_FD, kNCWrongSignBkg, kNumEnergyBinsFar, n, NueConvention::nue, and WhichFitParam().

Referenced by DoneFilling().

00684 {
00685   using Detector::kFar;
00686 
00687   //loop over NCEnergy bin objects in the far detector
00688   double trueEnergy    = 0.;
00689   double recoShowerE   = 0.;
00690   double recoMuonE     = 0.;
00691   double weight        = 0.;
00692   int    flavor        = 0;
00693 
00694   int sig         = kNCSig;
00695   int bkg         = kNCBkg;
00696   int tau         = kNCTau;
00697   int beamnue     = kNCBeamNue;
00698   int oscnue      = kNCOscNue;
00699   int wrongsignbkg = kNCWrongSignBkg;
00700 
00701   int selreco_FD = kNCSelRecoEnergy_FD;
00702   int reco_FD = kNCRecoEnergy_FD;
00703   
00704   int trueE = kNCTrueEnergy_FD;
00705   int trueExsec = kNCTrueEnergyXsecFit_FD;
00706   int trueEBkg = kNCTrueEnergyBkg_FD;
00707 
00708   int recototrue = kNCTrueRecoSig_FD;
00709   int recototrueBkg = kNCTrueRecoBkg_FD; 
00710   
00711   int recoBkg = kNCRecoBkg_FD;
00712   int purity = kNCPurity_FD;
00713 
00714   int nue = kNCRecoEnergyBeamNueBkg_FD ;
00715 
00716   if(nccc == NCType::kCC){
00717     sig      = kCCSig;
00718     bkg      = kCCBkg;
00719     tau      = kCCTau;
00720     beamnue  = kCCBeamNue;
00721     oscnue   = kCCOscNue;
00722     wrongsignbkg = kCCWrongSignBkg;
00723   
00724     selreco_FD = kCCSelRecoEnergy_FD;
00725     reco_FD = kCCRecoEnergy_FD;
00726     
00727     trueE = kCCTrueEnergy_FD;
00728     trueExsec = kCCTrueEnergyXsecFit_FD;
00729     trueEBkg = kCCTrueEnergyBkg_FD;
00730 
00731     recototrue = kCCTrueRecoSig_FD;
00732     recototrueBkg = kCCTrueRecoBkg_FD; 
00733 
00734     recoBkg = kCCRecoBkg_FD;
00735     purity = kCCPurity_FD;
00736 
00737     nue = kCCRecoEnergyBeamNueBkg_FD ;
00738 }
00739 
00740    bmtouse->Get2DOscHists()[sig]->Reset( "ICE" );
00741    bmtouse->Get2DOscHists()[bkg]->Reset( "ICE" );
00742    bmtouse->Get2DOscHists()[tau]->Reset( "ICE" );
00743    bmtouse->Get2DOscHists()[beamnue]->Reset( "ICE" );
00744    bmtouse->Get2DOscHists()[oscnue]->Reset( "ICE" );
00745    bmtouse->Get2DOscHists()[wrongsignbkg]->Reset( "ICE" );
00746 
00747    bmtouse->GetMCHists()[selreco_FD]->Reset( "ICE" );
00748    bmtouse->GetMCHists()[reco_FD]->Reset( "ICE" );
00749    
00750    bmtouse->GetMCHists()[trueE]->Reset( "ICE" );
00751    bmtouse->GetMCHists()[trueExsec]->Reset( "ICE" );
00752    bmtouse->GetMCHists()[trueEBkg]->Reset( "ICE" );
00753    
00754    bmtouse->Get2DHists()[recototrue]->Reset( "ICE" );
00755    bmtouse->Get2DHists()[recototrueBkg]->Reset( "ICE" );
00756 
00757    bmtouse->GetMCHists()[recoBkg]->Reset( "ICE" );
00758    bmtouse->GetMCHists()[purity]->Reset( "ICE" );
00759    
00760    bmtouse->GetMCHists()[nue]->Reset( "ICE" );
00761    
00762    const int B = beam->GetNumberEnergyBins(nccc);
00763    for(int b = 0; b < B; ++b){
00764      NCEnergyBin* bin = beam->GetEnergyBin(b, nccc);
00765      
00766      const int mcSize         = bin->GetMCSignalVectorSize();
00767      const int mcSize_bkg     = bin->GetMCBackgroundVectorSize();
00768      const int mcSize_tau     = bin->GetMCNuTauVectorSize();
00769      const int mcSize_beamnue = bin->GetMCBeamNuEVectorSize();
00770      const int mcSize_oscnue  = bin->GetMCOscNuEVectorSize();
00771      
00772      //numu
00773      for(int e = 0; e < mcSize; ++e){
00774        bin->GetMCInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00775        if(nccc == NCType::kNC) recoMuonE = 0;
00776        if(flavor == 14) FillHistogramSmoothed(bmtouse->GetMCHists()[trueE], trueEnergy, weight);           
00777        //The weights for the xsection go into the true to reco histogram and are applied at the true to reco stage. 
00778        //Relevant Backgrounds are adjusted for xsecs in the DoOscillations method
00779        //When fitting for a NC Clean systematic it only shifts the nc selected and wrong sign bkg NOT cc bkg... small effect.             
00780        if(flavor == 14)FillHistogramSmoothed2D( bmtouse->Get2DOscHists()[sig], recoShowerE + recoMuonE, trueEnergy, weight);      
00781        if(flavor != 14)FillHistogramSmoothed2D( bmtouse->Get2DOscHists()[wrongsignbkg], recoShowerE + recoMuonE, trueEnergy, weight);
00782        if(nccc == NCType::kNC) weight *= bmtouse->GetmScaleNC().at(WhichFitParam(trueEnergy));
00783        if(flavor == 14) FillHistogramSmoothed(bmtouse->GetMCHists()[purity], recoShowerE + recoMuonE, weight);    
00784        FillHistogramSmoothed(bmtouse->GetMCHists()[reco_FD], recoShowerE + recoMuonE, weight);
00785        if(flavor == 14) FillHistogramSmoothed2D(bmtouse->Get2DHists()[recototrue], trueEnergy, recoShowerE + recoMuonE, weight);
00786        if(flavor == 14) FillHistogramSmoothed(bmtouse->GetMCHists()[trueExsec], trueEnergy, weight);           
00787      }
00788 
00789      //numu bkg
00790      for(int e = 0; e < mcSize_bkg; ++e){
00791        bin->GetMCBackgroundInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00792        if(nccc == NCType::kNC) recoMuonE = 0;
00793        if(flavor == 14) FillHistogramSmoothed2D(bmtouse->Get2DHists()[recototrueBkg], trueEnergy, recoShowerE + recoMuonE, weight);
00794        if(flavor == 14) FillHistogramSmoothed(bmtouse->GetMCHists()[trueEBkg], trueEnergy, weight);     
00795        if(nccc == NCType::kCC) weight *= bmtouse->GetmScaleNC().at(WhichFitParam(trueEnergy));
00796        FillHistogramSmoothed(bmtouse->GetMCHists()[recoBkg], recoShowerE + recoMuonE, weight);      
00797        //Weighting this here rather than DoOsc method
00798        FillHistogramSmoothed2D( bmtouse->Get2DOscHists()[bkg], recoShowerE + recoMuonE, trueEnergy, weight);
00799      }
00800      
00801      //nutau
00802      for(int e = 0; e < mcSize_tau; ++e){
00803        bin->GetMCNuTauInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00804        if(nccc == NCType::kNC) recoMuonE = 0;
00805        FillHistogramSmoothed2D( bmtouse->Get2DOscHists()[tau], recoShowerE + recoMuonE, trueEnergy, weight);
00806      }
00807      
00808      //beamnue
00809      for(int e = 0; e < mcSize_beamnue; ++e){
00810        bin->GetMCBeamNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00811        if(nccc == NCType::kNC) recoMuonE = 0;
00812        FillHistogramSmoothed(bmtouse->GetMCHists()[nue], recoShowerE + recoMuonE, weight);
00813        FillHistogramSmoothed2D( bmtouse->Get2DOscHists()[beamnue], recoShowerE + recoMuonE, trueEnergy, weight);
00814      }
00815     
00816      //oscnue
00817      for(int e = 0; e < mcSize_oscnue; ++e){
00818        bin->GetMCOscNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00819        if(nccc == NCType::kNC) recoMuonE = 0;
00820        FillHistogramSmoothed2D( bmtouse->Get2DOscHists()[oscnue], recoShowerE + recoMuonE, trueEnergy, weight);
00821      }
00822     
00823      const int N = bin->GetDataVectorSize();
00824      for(int n = 0; n < N; ++n){
00825        double energy, weight;
00826        bin->GetDataInformation(energy, weight, n);
00827        FillHistogramSmoothed(bmtouse->GetFDDataHists()[sig], energy, weight);
00828      } // end for n
00829    } // end for b
00830 
00831    for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
00832      double denom = bmtouse->GetMCHists().at(reco_FD)->GetBinContent(i) + bmtouse->GetMCHists().at(recoBkg)->GetBinContent(i) + bmtouse->GetMCHists().at(nue)->GetBinContent(i) ;
00833      bmtouse->GetMCHists().at(selreco_FD)->SetBinContent(i, denom ) ;    
00834      if(denom > 0)
00835        bmtouse->GetMCHists().at(purity)->SetBinContent(i,bmtouse->GetMCHists().at(purity)->GetBinContent(i) / denom ) ;
00836      else bmtouse->GetMCHists().at(purity)->SetBinContent(i,0);
00837    }
00838 }

void NCExtrapolationBeamMatrix::FillDataMCHistogramsNear ( NCBeam beam,
NCType::EEventType  nccc,
BeamMatrix bmtouse,
bool  systFit = false 
) [private]

Definition at line 510 of file NCExtrapolationBeamMatrix.cxx.

References BMbybeams, MuELoss::e, fBMbybeams, FillHistogramSmoothed(), FillHistogramSmoothed2D(), BeamMatrix::Get2DHists(), NCEnergyBin::GetDataInformation(), BeamMatrix::GetDataPredHists(), NCEnergyBin::GetDataVectorSize(), NCBeam::Info::GetDescription(), NCBeam::GetEnergyBin(), NCBeam::GetInfo(), NCEnergyBin::GetMCBackgroundInformation(), NCEnergyBin::GetMCBackgroundVectorSize(), NCEnergyBin::GetMCBeamNuEInformation(), NCEnergyBin::GetMCBeamNuEVectorSize(), BeamMatrix::GetMCHists(), NCEnergyBin::GetMCInformation(), NCEnergyBin::GetMCNuTauInformation(), NCEnergyBin::GetMCNuTauVectorSize(), NCEnergyBin::GetMCOscNuEInformation(), NCEnergyBin::GetMCOscNuEVectorSize(), NCEnergyBin::GetMCSignalVectorSize(), BeamMatrix::GetmScaleNC(), NCBeam::GetNumberEnergyBins(), NCBeam::GetRunType(), NCType::kCC, kCCDataRecoEnergy_ND, kCCPurity_ND, kCCRecoBkg_ND, kCCRecoEnergy_ND, kCCRecoEnergyBeamNueBkg_ND, kCCRecoEnergyOscNue_ND, kCCRecoEnergyTau_ND, kCCRecoTrueBkg_ND, kCCRecoTrueSig_ND, kCCSelRecoEnergy_ND, kCCTrueEnergy_ND, kCCTrueEnergyBkg_ND, Msg::kDebug, BeamType::kL010z185i, NCType::kNC, kNCDataRecoEnergy_ND, kNCPurity_ND, kNCRecoBkg_ND, kNCRecoEnergy_ND, kNCRecoEnergyBeamNueBkg_ND, kNCRecoEnergyOscNue_ND, kNCRecoEnergyTau_ND, kNCRecoTrueBkg_ND, kNCRecoTrueSig_ND, kNCSelRecoEnergy_ND, kNCTrueEnergy_ND, kNCTrueEnergyBkg_ND, Detector::kNear, kNumEnergyBinsFar, MSG, Munits::second, and WhichFitParam().

Referenced by doNDXsectionFit(), and DoneFilling().

00514 {
00515   using Detector::kNear;
00516 
00517   //loop over NCEnergy bin objects in the near detector
00518   //to get the energies from the data events in each bin
00519   //and fill the histograms
00520 
00521   int selreco_ND = kNCSelRecoEnergy_ND;
00522   int reco_ND = kNCRecoEnergy_ND;
00523   int recoBkg = kNCRecoBkg_ND;
00524   int purity = kNCPurity_ND;
00525   int tau = kNCRecoEnergyTau_ND;
00526   int oscnue =  kNCRecoEnergyOscNue_ND ;
00527   int beamnue = kNCRecoEnergyBeamNueBkg_ND;
00528 
00529   int recototrue = kNCRecoTrueSig_ND;
00530   int recototrueBkg = kNCRecoTrueBkg_ND; 
00531 
00532   int recoData_ND = kNCDataRecoEnergy_ND;
00533   int trueE = kNCTrueEnergy_ND;
00534   int trueEBkg = kNCTrueEnergyBkg_ND;
00535 
00536   if(nccc == NCType::kCC){
00537 
00538     selreco_ND = kCCSelRecoEnergy_ND;
00539     reco_ND = kCCRecoEnergy_ND;
00540     recoBkg = kCCRecoBkg_ND;
00541     purity = kCCPurity_ND;
00542     tau = kCCRecoEnergyTau_ND;
00543     oscnue =  kCCRecoEnergyOscNue_ND ;
00544     beamnue = kCCRecoEnergyBeamNueBkg_ND;
00545 
00546     recototrue =  kCCRecoTrueSig_ND;
00547     recototrueBkg = kCCRecoTrueBkg_ND; 
00548 
00549     recoData_ND = kCCDataRecoEnergy_ND;
00550     trueE =  kCCTrueEnergy_ND;
00551     trueEBkg = kCCTrueEnergyBkg_ND;
00552 
00553   }
00554   
00555   //
00556   //Turns out that reseting is important
00557   //otherwise histograms pile up on each other
00558   //
00559   
00560   bmtouse->GetMCHists()[selreco_ND]->Reset( "ICE" );
00561   bmtouse->GetMCHists()[reco_ND]->Reset( "ICE" );
00562   bmtouse->GetMCHists()[recoBkg]->Reset( "ICE" );
00563   bmtouse->GetMCHists()[purity]->Reset( "ICE" );
00564   bmtouse->GetMCHists()[tau]->Reset( "ICE" );
00565   bmtouse->GetMCHists()[beamnue]->Reset( "ICE" );
00566   bmtouse->GetMCHists()[oscnue]->Reset( "ICE" );
00567   
00568   bmtouse->GetMCHists()[trueE]->Reset( "ICE" );
00569   bmtouse->GetMCHists()[trueEBkg]->Reset( "ICE" );
00570   
00571   bmtouse->Get2DHists()[recototrue]->Reset( "ICE" );
00572   bmtouse->Get2DHists()[recototrueBkg]->Reset( "ICE" );
00573   
00574   bmtouse->GetDataPredHists()[recoData_ND]->Reset( "ICE" );
00575 
00576   for(int b  = 0; b < beam->GetNumberEnergyBins(nccc); ++b){
00577     NCEnergyBin* bin = beam->GetEnergyBin(b, nccc);
00578 
00579     const int dataSize       = bin->GetDataVectorSize();
00580     const int mcSize         = bin->GetMCSignalVectorSize();
00581     const int mcSize_bkg     = bin->GetMCBackgroundVectorSize();
00582     const int mcSize_tau     = bin->GetMCNuTauVectorSize();
00583     const int mcSize_beamnue = bin->GetMCBeamNuEVectorSize();
00584     const int mcSize_oscnue  = bin->GetMCOscNuEVectorSize();
00585         
00586     // The variables that get filled with event information
00587     double trueEnergy    = 0;
00588     double recoShowerE   = 0;
00589     double recoMuonE     = 0;
00590     double energy        = 0;
00591     double weight        = 0;
00592     int    flavor        = 0;
00593 
00594     for(int e = 0; e < dataSize; ++e){
00595       bin->GetDataInformation(energy, weight, e);
00596       FillHistogramSmoothed(bmtouse->GetDataPredHists()[recoData_ND], energy, weight);
00597     }
00598 
00599     //------------------------------
00600     // numu
00601     //no nc weights * here as it doesn't matter. Only need CC spectrum to look right
00602     //in the BeamMatrix used in the actual extrapolation. The ND Fit BeamMatrix obviously
00603     //can't have the weights as they are used to calculate them.
00604     for(int e = 0; e < mcSize; ++e){
00605       bin->GetMCInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00606       if(nccc == NCType::kNC) recoMuonE = 0;
00607       FillHistogramSmoothed( bmtouse->GetMCHists()[reco_ND], recoShowerE + recoMuonE, weight);
00608       if(flavor == 14) FillHistogramSmoothed( bmtouse->GetMCHists()[purity], recoShowerE + recoMuonE, weight);    
00609       if(flavor == 14) FillHistogramSmoothed( bmtouse->GetMCHists()[trueE], trueEnergy, weight);    
00610       if(flavor == 14) FillHistogramSmoothed2D(bmtouse->Get2DHists()[recototrue], trueEnergy, recoShowerE + recoMuonE, weight);
00611     }
00612 
00613   // numubkg
00614     for(int e = 0; e < mcSize_bkg; ++e){
00615       bin->GetMCBackgroundInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00616       if(nccc == NCType::kNC) recoMuonE = 0;
00617       if(nccc == NCType::kCC) weight *= bmtouse->GetmScaleNC().at(WhichFitParam(trueEnergy)) ;      
00618       // need to adjust CCBkg by the nd xsections. The weights are default 1 for NDFit
00619       if(flavor == 14) FillHistogramSmoothed(bmtouse->GetMCHists()[trueEBkg], trueEnergy, weight);    
00620       FillHistogramSmoothed(bmtouse->GetMCHists()[recoBkg], recoShowerE + recoMuonE, weight);             
00621       FillHistogramSmoothed2D(bmtouse->Get2DHists()[recototrueBkg], trueEnergy, recoShowerE + recoMuonE ,weight);
00622     }
00623   
00624     //------------------------------
00625     // nutau
00626     for(int e = 0; e < mcSize_tau; ++e){
00627       bin->GetMCNuTauInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00628       if(nccc == NCType::kNC) recoMuonE = 0;
00629       FillHistogramSmoothed(bmtouse->GetMCHists()[tau], recoShowerE + recoMuonE, weight);
00630     }
00631 
00632     //------------------------------
00633     // beam nue
00634     for(int e = 0; e < mcSize_beamnue; ++e){
00635       bin->GetMCBeamNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00636       if(nccc == NCType::kNC) recoMuonE = 0;
00637       FillHistogramSmoothed(bmtouse->GetMCHists()[beamnue], recoShowerE + recoMuonE, weight);
00638     }
00639 
00640     //------------------------------
00641     // oscillated nue
00642     for(int e = 0; e < mcSize_oscnue; ++e){
00643       bin->GetMCOscNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00644       if(nccc == NCType::kNC) recoMuonE = 0;
00645       FillHistogramSmoothed(bmtouse->GetMCHists()[oscnue], recoShowerE + recoMuonE, weight);
00646     }
00647   } // end for b
00648   
00649   
00650   // copying the data into the systematic beamTypes.
00651   // NDFit has it's own way of doing this only need LE as that is all
00652   // thats extrapolated.   
00653   // works because nominal comes before systematics in beams vector.
00654   if(!systFit && !beam->GetInfo().IsNominal()){
00655 
00656     NC::RunUtil::ERunType runNum = beam->GetRunType();
00657     const NCBeam::Info beamInfo = NCBeam::Info(BeamType::kL010z185i, runNum);
00658     int thisBeam = fBMbybeams.find(beam->GetInfo())->second ;
00659     int nomBeam = fBMbybeams.find(beamInfo)->second ;
00660 
00661     MSG("NCExtrapolationBeamMatrix", Msg::kDebug) << "Filling Syst Beam " << beam->GetInfo().GetDescription() << " thisBeam " << thisBeam  
00662                                                   << " with data from  " << beamInfo.GetDescription() << " nomBeam " << nomBeam <<endl;
00663     for(int i = 1; i < kNumEnergyBinsFar+1; ++i )
00664       bmtouse->GetDataPredHists()[recoData_ND]->SetBinContent(i,BMbybeams[nomBeam]->GetDataPredHists()[recoData_ND]->GetBinContent(i));
00665     
00666   }
00667   
00668   for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
00669     // purity correction histogram
00670     double denom = bmtouse->GetMCHists().at(reco_ND)->GetBinContent(i) + bmtouse->GetMCHists().at(recoBkg)->GetBinContent(i) + bmtouse->GetMCHists().at(beamnue)->GetBinContent(i);
00671     //make NC/cc selected histos
00672     bmtouse->GetMCHists().at(selreco_ND)->SetBinContent(i, denom ) ;
00673     if(denom > 0)
00674       bmtouse->GetMCHists().at(purity)->SetBinContent(i,bmtouse->GetMCHists().at(purity)->GetBinContent(i) / denom ) ;
00675     else bmtouse->GetMCHists().at(purity)->SetBinContent(i,0);
00676   }
00677   
00678 }

void NCExtrapolationBeamMatrix::FillHistogramSmoothed ( TH1D *  h,
double  x,
double  weight 
) [inline, private]

Definition at line 497 of file NCExtrapolationBeamMatrix.cxx.

References HistSmoothHelper().

Referenced by FillDataMCHistogramsFar(), FillDataMCHistogramsNear(), and FillNDHistsForXSectionFit().

00498 {
00499   HistSmoothHelper(h, false, weight, x);
00500 }

void NCExtrapolationBeamMatrix::FillHistogramSmoothed2D ( TH2D *  h,
double  x,
double  y,
double  weight 
) [inline, private]

Definition at line 503 of file NCExtrapolationBeamMatrix.cxx.

References HistSmoothHelper().

Referenced by FillDataMCHistogramsFar(), and FillDataMCHistogramsNear().

00504 {
00505   HistSmoothHelper(h, true, weight, x, y);
00506 }

void NCExtrapolationBeamMatrix::FillNDHistsForXSectionFit ( NCBeam beam_nd,
NCType::EEventType  nccc,
bool  doSyst 
) [private]

Definition at line 2324 of file NCExtrapolationBeamMatrix.cxx.

References dataNDhist, MuELoss::e, FillHistogramSmoothed(), fluxCorrectionsNDNC, fluxCorrectionsNDNCbkg, NCBeam::GetBeamType(), NCEnergyBin::GetDataInformation(), NCEnergyBin::GetDataVectorSize(), NCBeam::GetEnergyBin(), NCEnergyBin::GetMCBackgroundInformation(), NCEnergyBin::GetMCBackgroundVectorSize(), NCEnergyBin::GetMCBeamNuEInformation(), NCEnergyBin::GetMCBeamNuEVectorSize(), NCEnergyBin::GetMCInformation(), NCEnergyBin::GetMCNuTauInformation(), NCEnergyBin::GetMCNuTauVectorSize(), NCEnergyBin::GetMCOscNuEInformation(), NCEnergyBin::GetMCOscNuEVectorSize(), NCEnergyBin::GetMCSignalVectorSize(), NCBeam::GetNumberEnergyBins(), Msg::kInfo, BeamType::kL010z185i, BeamType::kL100z200i, BeamType::kL250z200i, NCType::kNC, kNCDataHE, kNCDataLE, kNCDataME, kNCMCHE, kNCMCLE, kNCMCME, MAXMSG, ND_XsecByParams, ND_XSectionFit, ND_XSectionFitRew, NCExtrapolation::Reset(), WhichFitParam(), and WhichFluxCorrParam().

Referenced by ScaleBeamsToSameNumberOfEvents().

02326                                       {
02327 
02328   int data = -1; 
02329   int mc = -1;
02330   int beamindex = -1;
02331  
02332   if( beam_nd->GetBeamType() == BeamType::kL010z185i){
02333     data = kNCDataLE; 
02334     mc =  kNCMCLE;     
02335     beamindex = 0 ;
02336   }
02337  
02338   if( beam_nd->GetBeamType() == BeamType::kL100z200i){
02339     data = kNCDataME; 
02340     mc =  kNCMCME; 
02341     beamindex = 1;
02342   }
02343   
02344   if( beam_nd->GetBeamType() == BeamType::kL250z200i){
02345     data = kNCDataHE; 
02346     mc =  kNCMCHE;
02347     beamindex = 2;
02348   }
02349   //no cc at the moment. Not needed.
02350 
02351   ND_XSectionFit[data]->Reset("ICE");
02352   ND_XSectionFitRew[data]->Reset("ICE");
02353   for(UInt_t ipar = 0 ; ipar < ND_XsecByParams.at(0).size() ; ipar++)
02354     ND_XsecByParams[beamindex][ipar]->Reset("ICE");
02355   
02356   ND_XSectionFit[mc]->Reset("ICE");
02357 
02358   for(int b  = 0; b < beam_nd->GetNumberEnergyBins(nccc); ++b){
02359     NCEnergyBin* bin = beam_nd->GetEnergyBin(b, nccc);
02360 
02361     const int dataSize       = bin->GetDataVectorSize();
02362     const int mcSize         = bin->GetMCSignalVectorSize();
02363     const int mcSize_bkg     = bin->GetMCBackgroundVectorSize();
02364     const int mcSize_tau     = bin->GetMCNuTauVectorSize();
02365     const int mcSize_beamnue = bin->GetMCBeamNuEVectorSize();
02366     const int mcSize_oscnue  = bin->GetMCOscNuEVectorSize();
02367     
02368     MAXMSG("NCExtrapolationBeamMatrix",Msg::kInfo,3) << "data " << data << " data size " << dataSize  <<" mc  " << mc << " mc size" << mcSize << endl;
02369      
02370 
02371     // The variables that get filled with event information
02372     double trueEnergy    = 0;
02373     double recoShowerE   = 0;
02374     double recoMuonE     = 0;
02375     double energy        = 0;
02376     double weight        = 0;
02377     int flavor = 0;
02378     for(int e = 0; e < dataSize; ++e){
02379       bin->GetDataInformation(energy, weight, e);
02380       FillHistogramSmoothed(ND_XSectionFit[data], energy, weight);
02381       FillHistogramSmoothed(ND_XSectionFitRew[data], energy, weight);
02382     }
02383     //------------------------------
02384     // numu
02385     for(int e = 0; e < mcSize; ++e){ 
02386       bin->GetMCInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
02387       if(nccc == NCType::kNC){ 
02388         recoMuonE = 0;
02389         weight *= fluxCorrectionsNDNC[beamindex][WhichFluxCorrParam(trueEnergy)] ;
02390       }
02391       FillHistogramSmoothed(ND_XsecByParams[beamindex][WhichFitParam(trueEnergy)], recoShowerE + recoMuonE, weight);
02392     }
02393     //bkg
02394     for(int e = 0; e < mcSize_bkg; ++e){
02395       bin->GetMCBackgroundInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
02396       if(nccc == NCType::kNC) {
02397         recoMuonE = 0;
02398         weight *= fluxCorrectionsNDNCbkg[beamindex][WhichFluxCorrParam(trueEnergy)] ;
02399       }
02400       FillHistogramSmoothed(ND_XsecByParams[beamindex][4], recoShowerE + recoMuonE, weight) ;
02401     }
02402 
02403     //------------------------------
02404     // nutau
02405     for(int e = 0; e < mcSize_tau; ++e){
02406       bin->GetMCNuTauInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
02407       if(nccc == NCType::kNC) recoMuonE = 0;
02408       FillHistogramSmoothed(ND_XsecByParams[beamindex][4], recoShowerE + recoMuonE, weight) ;
02409     }
02410 
02411     //------------------------------
02412     // beam nue
02413     for(int e = 0; e < mcSize_beamnue; ++e){
02414       bin->GetMCBeamNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
02415       if(nccc == NCType::kNC)recoMuonE = 0;
02416       FillHistogramSmoothed(ND_XsecByParams[beamindex][4], recoShowerE + recoMuonE, weight) ;
02417     }
02418     
02419     //------------------------------
02420     // oscillated nue
02421     for(int e = 0; e < mcSize_oscnue; ++e){
02422       bin->GetMCOscNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
02423       if(nccc == NCType::kNC) recoMuonE = 0;
02424       FillHistogramSmoothed(ND_XsecByParams[beamindex][4], recoShowerE + recoMuonE, weight) ;
02425     }
02426   } // end for b
02427 
02428   
02429   for(int i = 1 ; i < kNumEnergyBinsFar+1 ; i++){
02430     double binsum = 0.0;
02431     for(UInt_t ipar = 0 ; ipar < ND_XsecByParams.at(0).size() ; ipar++){
02432       binsum += ND_XsecByParams[beamindex][ipar]->GetBinContent(i);
02433     }
02434     ND_XSectionFit[mc]->SetBinContent(i,binsum);
02435   }
02436 
02437   if(doSyst)    
02438     for(int i = 1 ; i < kNumEnergyBinsFar+1 ; i++){
02439       ND_XSectionFit[data]->SetBinContent(i,dataNDhist[beamindex][0]->GetBinContent(i));
02440       ND_XSectionFitRew[data]->SetBinContent(i,dataNDhist[beamindex][0]->GetBinContent(i));
02441     }
02442 }

void NCExtrapolationBeamMatrix::FillResultHistograms ( NCBeam::Info  beamInfo,
NCType::EEventType  nccc,
std::vector< std::vector< double > > &  nominal,
std::vector< std::vector< double > > &  osc 
) [private]

Referenced by WriteResources().

virtual void NCExtrapolationBeamMatrix::FindSpectraForPars ( const NC::OscProb::OscPars oscPars,
const NC::SystPars systPars,
std::vector< NCBeam::Info beamsToUse,
std::vector< TH1 * > &  exps,
std::vector< TH1 * > &  obss 
) [private, virtual]

Override this in the derived class.

Parameters:
oscPars Oscillation parameters
systPars Systematic shifts
beamsToUse Only fill exps and obss with spectra from beams with these indices
[out] exps To be filled with expected spectra
[out] obss To be filled with the corresponding observed spectra

Implements NCExtrapolation.

void NCExtrapolationBeamMatrix::Fit ( int  PrintLevel  )  [private]

Definition at line 2158 of file NCExtrapolationBeamMatrix.cxx.

References corr, Form(), Msg::kInfo, logFile, LogLikelihoodFunc(), m0ScaleNC, m0ScaleNCErr, mScaleNC, mScaleNCErr, mScaleNeg, mScalePos, and MSG.

Referenced by doNDXsectionFit().

02158                    {
02159 
02160   Double_t arglist[10];
02161   Int_t    ierr;
02162 
02163   TMinuit* gMinuit = new TMinuit(4);
02164 
02165   gMinuit->SetFCN(LogLikelihoodFunc); // Likelihood fit
02166 
02167   arglist[0] = PrintLevel;
02168   gMinuit->mnexcm("SET PRInt",arglist,1,ierr);
02169   if ( PrintLevel<-1 ) 
02170     gMinuit->mnexcm("SET NOWarning",arglist,0,ierr);
02171 
02172   arglist[0] = 1; // added a factor 2 to the likelihood function
02173   gMinuit->mnexcm("SET ERRor",arglist,1,ierr);
02174   
02175   arglist[0] = 1; // minimization strategy
02176   gMinuit->mnexcm("SET STR",arglist,1,ierr);
02177   
02178   // Set starting values
02179   Double_t par[4];
02180   for (Int_t i=0; i<4; i++)
02181     par[i] = m0ScaleNC.at(i);
02182 
02183   // set starting error
02184   Double_t step[4];
02185   for (Int_t i=0; i<4; i++)
02186     step[i] = m0ScaleNCErr.at(i);
02187 
02188   // define parameters
02189   for (Int_t i=0; i<4; i++) {
02190     char* name = Form("scaleNC_%i",i);
02191     //gMinuit->DefineParameter(i,name, par[i], step[i],0., 0.); // unbounded
02192     gMinuit->DefineParameter(i,name, par[i], step[i],0., 5.0); // bounded
02193   }
02194 
02195   MSG("NCExtrapolationBeamMatrix",Msg::kInfo) << "Doing BOUNDED FIT 0 to 5 " << endl;
02196   // start restricted fit first
02197  
02198   gMinuit->FixParameter(0);
02199   gMinuit->FixParameter(1);
02200   gMinuit->FixParameter(2);
02201 
02202   // do a scan
02203  //  arglist[0] = 5; // scan parameter no 4+1
02204 //   arglist[1] = 40;
02205 //   arglist[2] = 0.8;
02206 //   arglist[3] = 1.0;
02207 //   gMinuit->mnexcm("SCAn",arglist,4,ierr);
02208 //   TCanvas *cnew = new TCanvas("cnew","cnew");
02209 //   TGraph *gr = (TGraph*)gMinuit->GetPlot();
02210 //   cnew->cd();
02211 //   gr->Draw("ap");
02212   
02213   gMinuit->Migrad();
02214 
02215   gMinuit->Release(2);  
02216   gMinuit->Migrad();
02217 
02218   gMinuit->Release(1);  
02219   gMinuit->Migrad();
02220 
02221   gMinuit->Release(0);  
02222   gMinuit->Migrad();
02223 
02224   // do MINOS error analysis
02225   gMinuit->mnmnos();
02226 
02227   // check convergence (still to do)
02228   TString convergence = gMinuit->fCstatu; 
02229   MSG("NCExtrapolationBeamMatrix",Msg::kInfo) << "Convergence status: "
02230                             << convergence << endl;
02231 
02232   // get result
02233   for (Int_t i=0; i<4; i++)
02234     gMinuit->GetParameter( i, mScaleNC.at(i),  mScaleNCErr.at(i) );
02235 
02236   //Double_t mScaleNeg[4]; // positive and negative errors 
02237                                           // from MINOS
02238   //Double_t mScalePos[4];
02239   //Double_t corr[4];
02240   // get proper asymmetric errors from MINOS fit
02241   for (Int_t i=0; i<4; i++)
02242     gMinuit->mnerrs(i,mScalePos[i], mScaleNeg[i],
02243                     mScaleNCErr[i],
02244                     corr[i] );
02245 
02246   // plot result  
02247   
02248   if ( PrintLevel >-2 ) {
02249 
02250     Double_t LL;
02251     Double_t pars[4];
02252     Int_t    npar = 4;
02253 
02254     for (Int_t i=0; i<4; i++)
02255       pars[i] = mScaleNC.at(i);
02256     //should return what you put in
02257     LogLikelihoodFunc(npar, pars, LL, pars, 0);
02258 
02259 
02260     MSG("NCExtrapolationBeamMatrix",Msg::kInfo) 
02261       << "The final likelihood is: LL=" << setprecision(10) << LL << endl;
02262     MSG("NCExtrapolationBeamMatrix",Msg::kInfo) 
02263       << "Fitted Parameters: 0-4 " << mScaleNC.at(0)  
02264       << " 4-8 " << mScaleNC.at(1)  
02265       << " 8-15 " << mScaleNC.at(2) 
02266       << " over 15 " << mScaleNC.at(3) << endl;
02267 
02268     // write some stuff to a log file
02269     *logFile << "Fit results:" << endl;   
02270     logFile->setf(ios::fixed);
02271     logFile->precision(5);
02272     for (Int_t i=0; i<4; i++) 
02273       *logFile << "Par " << i << ": " << mScaleNC[i] << " +"
02274                << mScalePos[i] << " " << mScaleNeg[i] << " +/- "
02275                << mScaleNCErr.at(i) << ", Corr: " << corr[i] << endl;
02276   
02277     *logFile << "Convergence Status: " << convergence << endl
02278              << "Final Likelihood LL=" << LL << endl;
02279       // << "Written histograms to " << histofile->GetName() << endl;
02280     logFile->unsetf(ios::fixed);
02281     
02282     /*
02283     //not using yet...
02284     // get the number of d.o.f. from chi2 fit
02285     Int_t dof = 0;
02286     if (fUseChi2Function) {
02287       for (UInt_t i = 0; i<fDataFraction.size(); ++i ) {
02288         Double_t chi2 = -1;
02289         Int_t ndf   -1;
02290         Int_t igood = 0;
02291         HistEvisNCd.at(i)->Chi2TestX(HistEvisNCrew.at(i),chi2,ndf,igood, 
02292                                      "UUCHI2");
02293         dof += ndf;
02294         MSG("NCExtrapolationRS",Msg::kDebug) << "Ndf for beam " << i << ": " 
02295                                              << ndf << endl; 
02296         if (igood != 0)
02297           MSG("NCExtrapolationRS",Msg::kWarning) << "Chi2 igood parameter" 
02298                                                  <<" is non-zer: "
02299                                        << igood << endl;
02300       }
02301     }
02302     */
02303 
02304     //This bit of code may be useful if one wants to write out a result tree
02305 
02306     // write same info to fit tree
02307     //    fFitRes.likelihood = LL;
02308     //     for (Int_t i=0; i<4; i++) {
02309     //       fFitRes.parNo = i;
02310     //       fFitRes.inputScale = mInputScale[i];
02311     //       fFitRes.par = (Float_t) mScaleNC[i];
02312     //       fFitRes.err = (Float_t) mScaleNCErr[i];
02313     //       fFitRes.errPos = (Float_t) mScalePos[i];
02314     //       fFitRes.errNeg = (Float_t) mScaleNeg[i];
02315     //       fFitRes.dof = dof;
02316     //       fFitResTree->Fill();    
02317     //     }
02318   }  
02319 
02320   delete gMinuit;
02321 }

void NCExtrapolationBeamMatrix::GetBeamMatrix ( NCBeam beam,
BeamMatrix bmtouse 
) [private]

Definition at line 1505 of file NCExtrapolationBeamMatrix.cxx.

References fBeamMatrixPath, NCBeam::GetRunType(), Msg::kError, Msg::kInfo, NC::RunUtil::kRunAll, NC::RunUtil::kRunI, NC::RunUtil::kRunII, NC::RunUtil::kRunIII, MSG, and BeamMatrix::SetBeamMatrix().

Referenced by doNDXsectionFit(), and DoneFilling().

01505                                                 {
01506 
01507   NC::RunUtil::ERunType runNum = beam->GetRunType();
01508   TString runstring = "" ;
01509   switch(runNum){
01510   case NC::RunUtil::kRunAll:
01511     runstring = "RunI" ; 
01512     MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Using Run Type " << runstring << " Beam Matrix" <<  endl;
01513     break;
01514   case NC::RunUtil::kRunI:
01515     runstring = "RunI" ; 
01516     MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Using Run Type " << runstring << " Beam Matrix" << endl;
01517     break;
01518   case NC::RunUtil::kRunII: 
01519     runstring = "RunII"; 
01520     MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Using Run Type " << runstring << " Beam Matrix" << endl;
01521     break;
01522   case NC::RunUtil::kRunIII: 
01523     runstring = "RunIII"; 
01524     MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Using Run Type " << runstring << " Beam Matrix" << endl;
01525     break;
01526   default:    
01527     runstring = "NORUN";
01528     MSG("NCExtrapolationBeamMatrix", Msg::kError) << "Using Run Type " << runstring << " Beam Matrix" << endl;
01529     }
01530  
01531   MSG("NCExtrapolationBeamMatrix", Msg::kInfo) <<  "Getting Beam Matrix from: " << fBeamMatrixPath << endl;
01532   TFile* beamMatrixFile = new TFile(fBeamMatrixPath,"READ");
01533   if(beamMatrixFile){
01534     bmtouse->SetBeamMatrix((TH2D*)beamMatrixFile->Get("FDVsNDMatrixXSecRW"+runstring));
01535   }
01536   
01537   else{ 
01538     MSG("NCExtrapolationBeamMatrix", Msg::kError) <<  "Can't get BeamMatrix from " << fBeamMatrixPath  << endl;
01539   }
01540 }

void NCExtrapolationBeamMatrix::GetEfficiencyHistograms ( NCBeam beam,
BeamMatrix bmtouse 
) [private]

Definition at line 1350 of file NCExtrapolationBeamMatrix.cxx.

References fExtractionType, fSelEffFilePath, NCBeam::GetBeamType(), NCBeam::GetInfo(), NCBeam::GetRunType(), Msg::kError, kExtractNCFromCCFlux, kExtractNCFromCCFluxND, Msg::kInfo, BeamType::kL010z185i, BeamType::kL100z200i, BeamType::kL250z200i, NCType::kNCCCExtractionANN, NCType::kNCCCExtractionANNFar, NCType::kNCCCExtractionANNNear, NCType::kNCCCExtractionCuts, NCType::kNCCCExtractionCutswkNN, kRecoCCEffFDHist, kRecoCCEffNDHist, kRecoNCEffFDHist, kRecoNCEffNDHist, NC::RunUtil::kRunAll, NC::RunUtil::kRunI, NC::RunUtil::kRunII, NC::RunUtil::kRunIII, kSelCCEffFDHist, kSelCCEffNDHist, kSelCCinNCEffFDHist, kSelCCinNCEffNDHist, kSelNCEffFDHist, kSelNCEffNDHist, MSG, and BeamMatrix::SetEffCorrHist().

Referenced by doNDXsectionFit(), and DoneFilling().

01350                                                           {
01351 
01352   MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Getting Efficiency Histograms from: " << fSelEffFilePath << " Extraction Type " << fExtractionType << endl
01353                                                << "For Run "<< beam->GetRunType() 
01354                                                << " Description " << beam->GetInfo().GetDescription() << endl;  
01355   TFile* selEffFile = new TFile(fSelEffFilePath,"READ");
01356   if(selEffFile){
01357     
01358     NC::RunUtil::ERunType runNum = beam->GetRunType();
01359     TString runstring = "" ;
01360     switch(runNum){
01361     case NC::RunUtil::kRunAll:
01362       runstring = "RunAll" ; 
01363       MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Using Run Type " << runstring ;
01364       break;   
01365     case NC::RunUtil::kRunI:
01366       runstring = "RunI" ; 
01367       MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Using Run Type " << runstring ;
01368       break;
01369     case NC::RunUtil::kRunII: 
01370       runstring = "RunII"; 
01371       MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Using Run Type " << runstring ;
01372       break;
01373     case NC::RunUtil::kRunIII: 
01374       runstring = "RunIII"; 
01375       MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Using Run Type " << runstring ;
01376       break;
01377     default:    
01378       runstring = "NORUN GIVEN";
01379       MSG("NCExtrapolationBeamMatrix", Msg::kError) << "Using Run Type " << runstring ;
01380     }
01381   
01382     TString ndbeamstring = "";
01383     BeamType::BeamType_t bt = beam->GetBeamType();
01384     switch(bt){
01385     case BeamType::kL010z185i: 
01386       ndbeamstring = "LE" ; 
01387       MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << " Using Beam Type " << ndbeamstring  << " for Efficiency Correction" << endl;
01388       break;   
01389     case BeamType::kL100z200i:
01390       ndbeamstring = "ME" ; 
01391       MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << " Using Beam Type " << ndbeamstring << " for Efficiency Correction" << endl;
01392       break;   
01393     case BeamType::kL250z200i:
01394       ndbeamstring = "HE" ; 
01395       MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << " Using Beam Type " << ndbeamstring << " for Efficiency Correction" << endl;
01396       break;   
01397     default:    
01398       runstring = "NO MATCHING BEAM GIVEN";
01399       MSG("NCExtrapolationBeamMatrix", Msg::kError) << " Using Beam Type " << ndbeamstring << " for Efficiency Correction" << endl;
01400     }
01401 
01402 
01403     switch(fExtractionType){
01404     case NCType::kNCCCExtractionCutswkNN:
01405       bmtouse->SetEffCorrHist(kSelNCEffNDHist,(TH1D*)selEffFile->Get("selCorrectionNCCutswkNNND"+runstring+ndbeamstring));
01406       bmtouse->SetEffCorrHist(kSelCCEffNDHist,(TH1D*)selEffFile->Get("selCorrectionCCCutswkNNND"+runstring+ndbeamstring));
01407       bmtouse->SetEffCorrHist(kSelNCEffFDHist,(TH1D*)selEffFile->Get("selCorrectionNCCutswkNNFD"+runstring+ndbeamstring)); 
01408       bmtouse->SetEffCorrHist(kSelCCEffFDHist,(TH1D*)selEffFile->Get("selCorrectionCCCutswkNNFD"+runstring+ndbeamstring));
01409       bmtouse->SetEffCorrHist(kSelCCinNCEffNDHist,(TH1D*)selEffFile->Get("selCorrectionCCinNCCutswkNNND"+runstring+ndbeamstring)); 
01410       bmtouse->SetEffCorrHist(kSelCCinNCEffFDHist,(TH1D*)selEffFile->Get("selCorrectionCCinNCCutswkNNFD"+runstring+ndbeamstring)); 
01411       break;    
01412     case NCType::kNCCCExtractionCuts:
01413       bmtouse->SetEffCorrHist(kSelNCEffNDHist,(TH1D*)selEffFile->Get("selCorrectionNCTOmND"+runstring+ndbeamstring));
01414       bmtouse->SetEffCorrHist(kSelCCEffNDHist,(TH1D*)selEffFile->Get("selCorrectionCCTOmND"+runstring+ndbeamstring));
01415       bmtouse->SetEffCorrHist(kSelNCEffFDHist,(TH1D*)selEffFile->Get("selCorrectionNCTOmFD"+runstring+ndbeamstring)); 
01416       bmtouse->SetEffCorrHist(kSelCCEffFDHist,(TH1D*)selEffFile->Get("selCorrectionCCTOmFD"+runstring+ndbeamstring)); 
01417       bmtouse->SetEffCorrHist(kSelCCinNCEffNDHist,(TH1D*)selEffFile->Get("selCorrectionCCinNCTOmND"+runstring+ndbeamstring)); 
01418       bmtouse->SetEffCorrHist(kSelCCinNCEffFDHist,(TH1D*)selEffFile->Get("selCorrectionCCinNCTOmFD"+runstring+ndbeamstring)); 
01419       break;
01420     case NCType::kNCCCExtractionANN:
01421     case NCType::kNCCCExtractionANNNear:
01422     case NCType::kNCCCExtractionANNFar:
01423       bmtouse->SetEffCorrHist(kSelNCEffNDHist,(TH1D*)selEffFile->Get("selCorrectionNCANNND"+runstring+ndbeamstring));
01424       bmtouse->SetEffCorrHist(kSelCCEffNDHist,(TH1D*)selEffFile->Get("selCorrectionCCANNND"+runstring+ndbeamstring));
01425       bmtouse->SetEffCorrHist(kSelNCEffFDHist,(TH1D*)selEffFile->Get("selCorrectionNCANNFD"+runstring+ndbeamstring)); 
01426       bmtouse->SetEffCorrHist(kSelCCEffFDHist,(TH1D*)selEffFile->Get("selCorrectionCCANNFD"+runstring+ndbeamstring)); 
01427       bmtouse->SetEffCorrHist(kSelCCinNCEffNDHist,(TH1D*)selEffFile->Get("selCorrectionCCinNCANNND"+runstring+ndbeamstring)); 
01428       bmtouse->SetEffCorrHist(kSelCCinNCEffFDHist,(TH1D*)selEffFile->Get("selCorrectionCCinNCANNFD"+runstring+ndbeamstring)); 
01429       break;
01430     default:
01431       bmtouse->SetEffCorrHist(kSelNCEffNDHist,0);
01432       bmtouse->SetEffCorrHist(kSelCCEffNDHist,0);
01433       bmtouse->SetEffCorrHist(kSelNCEffFDHist,0);
01434       bmtouse->SetEffCorrHist(kSelCCEffFDHist,0);
01435       bmtouse->SetEffCorrHist(kSelCCinNCEffNDHist,0); 
01436       bmtouse->SetEffCorrHist(kSelCCinNCEffFDHist,0); 
01437     }    
01438     
01439     bmtouse->SetEffCorrHist(kRecoNCEffNDHist,(TH1D*)selEffFile->Get("recoNCEffCorrHistND"+runstring+ndbeamstring));
01440     bmtouse->SetEffCorrHist(kRecoCCEffNDHist,(TH1D*)selEffFile->Get("recoCCEffCorrHistND"+runstring+ndbeamstring));
01441     bmtouse->SetEffCorrHist(kRecoNCEffFDHist,(TH1D*)selEffFile->Get("recoNCEffCorrHistFD"+runstring+ndbeamstring));
01442     bmtouse->SetEffCorrHist(kRecoCCEffFDHist,(TH1D*)selEffFile->Get("recoCCEffCorrHistFD"+runstring+ndbeamstring));
01443 
01444     bmtouse->SetEffCorrHist(kExtractNCFromCCFluxND,(TH1D*)selEffFile->Get("extractNCFromCCFluxND"+runstring+ndbeamstring));
01445     bmtouse->SetEffCorrHist(kExtractNCFromCCFlux,(TH1D*)selEffFile->Get("extractNCFromCCFlux"+runstring+ndbeamstring));
01446   }
01447   
01448   else{ 
01449     MSG("NCExtrapolationBeamMatrix", Msg::kError) <<  "Can't get Selection Efficiency Histograms from " 
01450                                                   << fSelEffFilePath  <<  " Extraction Type " <<  fExtractionType << endl;
01451   }
01452 }

virtual TString NCExtrapolationBeamMatrix::GetLongName (  )  const [inline, virtual]

This is the name the extrapolation is known under on plots and such.

Implements NCExtrapolation.

Definition at line 46 of file NCExtrapolationBeamMatrix.h.

00046 {return "BeamMatrix";}

void NCExtrapolationBeamMatrix::GetNCFromCCFlux ( BeamMatrix bmtouse  )  [private]

Definition at line 1568 of file NCExtrapolationBeamMatrix.cxx.

References BeamMatrix::GetDataPredHists(), BeamMatrix::GetEffCorrHists(), BeamMatrix::GetMCHists(), kCCDataTrueEnergyMatrix_FD, kCCTrueEnergyMatrix_FD, kExtractNCFromCCFlux, kNCDataTrueEnergyMatrix_FD, and kNCTrueEnergyMatrix_FD.

Referenced by ConstructFarSpectrum().

01568                                     {
01569   
01570   for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
01571 
01572     //not really sure where to put this either here or any of the following steps.
01573     //But the histograms already exist so it may as well be here. 
01574     bmtouse->GetMCHists().at(kNCTrueEnergyMatrix_FD)->SetBinContent(i, bmtouse->GetMCHists().at(kCCTrueEnergyMatrix_FD)->GetBinContent(i) 
01575                                                                     * bmtouse->GetEffCorrHists().at(kExtractNCFromCCFlux)->GetBinContent(i));
01576     bmtouse->GetDataPredHists().at(kNCDataTrueEnergyMatrix_FD)->SetBinContent(i, bmtouse->GetDataPredHists().at(kCCDataTrueEnergyMatrix_FD)->GetBinContent(i) 
01577                                                                               * bmtouse->GetEffCorrHists().at(kExtractNCFromCCFlux)->GetBinContent(i));
01578   }
01579 }

void NCExtrapolationBeamMatrix::GetNCFromCCFluxND ( BeamMatrix bmtouse  )  [private]

Definition at line 2673 of file NCExtrapolationBeamMatrix.cxx.

References BeamMatrix::GetEffCorrHists(), BeamMatrix::GetMCHists(), kCCTrueEnergyEffCorr_ND, kExtractNCFromCCFluxND, and kNCTrueEnergyEffCorrFromCC_ND.

Referenced by correctNDNCFlux().

02673                                       {
02674   for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
02675     //not really sure where to put this either here or any of the following steps.
02676     //But the histograms already exist so it may as well be here.
02677     bmtouse->GetMCHists().at(kNCTrueEnergyEffCorrFromCC_ND)->SetBinContent(i, bmtouse->GetMCHists().at(kCCTrueEnergyEffCorr_ND)->GetBinContent(i) 
02678                                                                                  * bmtouse->GetEffCorrHists().at(kExtractNCFromCCFluxND)->GetBinContent(i));
02679   }
02680 }

std::vector<const TH1* > NCExtrapolationBeamMatrix::GetNearMCSpectra ( std::vector< NCBeam::Info beamsToUse  )  [virtual]

Reimplemented from NCExtrapolation.

NC::RunUtil::ERunType NCExtrapolationBeamMatrix::GetRunTypeFromBM ( BeamMatrix BM  )  [private]

Definition at line 2756 of file NCExtrapolationBeamMatrix.cxx.

References NC::RunUtil::kRunAll, NC::RunUtil::kRunI, NC::RunUtil::kRunII, and NC::RunUtil::kRunIII.

Referenced by DoEfficencyCorrectionFD(), DoOscillations(), and DoUnoscTrueToRecoFD().

02756                                 {
02757 
02758   TString name = BM->GetName();
02759 
02760   //order matters here
02761   if(name.Contains("All Runs")) return NC::RunUtil::kRunAll ;
02762   else if(name.Contains("Run III")) return NC::RunUtil::kRunIII ;
02763   else if(name.Contains("Run II")) return NC::RunUtil::kRunII ;
02764   else if(name.Contains("Run I")) return NC::RunUtil::kRunI ;
02765 
02766   else assert(0 && "Failing to get Run Type"); 
02767 
02768 }

virtual TString NCExtrapolationBeamMatrix::GetShortName (  )  const [inline, virtual]

This is the name used to name things in the output file etc.

Implements NCExtrapolation.

Definition at line 45 of file NCExtrapolationBeamMatrix.h.

00045 {return "BeamMatrix";}

template<class T >
void NCExtrapolationBeamMatrix::HistSmoothHelper ( T *  h,
bool  is2D,
double  weight,
double  x,
double  y = -1 
) [inline, private]

Definition at line 425 of file NCExtrapolationBeamMatrix.cxx.

References fSmoothWidth.

Referenced by FillHistogramSmoothed(), and FillHistogramSmoothed2D().

00426 {
00427   if(fSmoothWidth == 0){
00428     // Shortcut all the below and just fill the histogram
00429     if(is2D)
00430       ((TH2D*)h)->Fill(x, y, weight);
00431     else
00432       h->Fill(x, weight);
00433     return;
00434   }
00435 
00436   // Paranoia over some scary text in the documentation for TAxis::GetBin
00437   h->SetBit(TH1::kCanRebin, false);
00438 
00439   TAxis* xaxis = h->GetXaxis();
00440   const int nbinsx = xaxis->GetNbins();
00441 
00442   // Warning - this number will only work in axis calls, not histogram calls
00443   const int xbin = xaxis->FindBin(x);
00444   const double lo = xaxis->GetBinLowEdge(xbin);
00445   const double hi = xaxis->GetBinUpEdge(xbin);
00446 
00447   const int ybin = is2D ? h->GetYaxis()->FindBin(y) : 0;
00448 
00449   const int histBin = xbin + (nbinsx+2)*ybin;
00450 
00451   const double dist_lo = x-lo;
00452   const double dist_hi = hi-x;
00453 
00454   // Don't smear events into the underflow bin
00455   const bool do_lo = xbin > 1 && dist_lo < fSmoothWidth*x;
00456   // Don't smear events into the overflow bin
00457   const bool do_hi = xbin < nbinsx && dist_hi < fSmoothWidth*x;
00458 
00459   double weight_here = 1;
00460   double weight_lo = 0;
00461   double weight_hi = 0;
00462 
00463   if(do_lo){
00464     const double transfer = .5*(1-dist_lo/(fSmoothWidth*x));
00465     weight_here -= transfer;
00466     weight_lo += transfer;
00467   }
00468   if(do_hi){
00469     const double transfer = .5*(1-dist_hi/(fSmoothWidth*x));
00470     weight_here -= transfer;
00471     weight_hi += transfer;
00472   }
00473 
00474   assert(weight_here >= .499);
00475   assert(weight_lo <= .501);
00476   assert(weight_hi <= .501);
00477   assert(weight_here <= 1);
00478   assert(weight_lo >= 0);
00479   assert(weight_hi >= 0);
00480   assert(fabs(weight_here+weight_lo+weight_hi-1) < .001);
00481 
00482   h->fArray[histBin] += weight*weight_here;
00483 
00484   if(do_lo){
00485     h->fArray[histBin-1] += weight*weight_lo;
00486   }
00487   if(do_hi){
00488     h->fArray[histBin+1] += weight*weight_hi;
00489   }
00490 }

void NCExtrapolationBeamMatrix::LogLikelihoodFunc ( Int_t &  npar,
Double_t *  gin,
Double_t &  f,
Double_t *  pars,
Int_t  iflag 
) [static, private]

Definition at line 2475 of file NCExtrapolationBeamMatrix.cxx.

References beamsFit, bmFit, Msg::kInfo, kNCDataHE, kNCDataLE, kNCDataME, kNCMCHE, kNCMCLE, kNCMCME, MAXMSG, MSG, ND_XsecByParams, ND_XSectionFitRew, and size.

Referenced by Fit().

02480   {  
02481   iflag++; // NOT USED FOR ANYTHING, JUST TO QUIET THE COMPILER
02482   if (0)  cout << gin << npar << endl; 
02483 
02484   f = 0;
02485 
02486 //   //  cout << "Parameter 5: " << pars[4]; 
02487   
02488 //   // cross section reweighting MC using neugen parameters
02489 //   //  sFit->Reweight(ma_qe,ma_res,kno112,kno122,kno113,kno123,false);
02490   
02491 // reset reweighted histograms
02492   for(UInt_t i = 1; i < bmFit->ND_XSectionFitRew.size(); i+=2){
02493     bmFit->ND_XSectionFitRew.at(i)->Reset("ICE");
02494   }
02495   
02496   // now fill Reweighted Histograms
02497   // do something faster than histogram filling
02498   
02499   for(UInt_t i = 0; i < (bmFit->beamsFit).size(); ++i){
02500     //maybe some check on nd only here    
02501     int data = -1 ;
02502     int mc = -1;
02503     //int beamindex = -1;
02504     
02505     if(i == 0){data = kNCDataLE; mc = kNCMCLE; } 
02506     if(i == 1){data = kNCDataME; mc = kNCMCME; } 
02507     if(i == 2){data = kNCDataHE; mc = kNCMCHE; } 
02508 
02509     // loop over the 4 energy ranges in the true NC selected spectrum
02510     for(int ibin = 1 ; ibin < kNumEnergyBinsFar+1 ; ibin++){
02511       double binsum = 0.0;
02512       for(UInt_t ipar = 0 ; ipar <  bmFit->ND_XsecByParams.at(0).size()-1 ; ipar++){
02513         binsum +=  bmFit->ND_XsecByParams[i][ipar]->GetBinContent(ibin) * pars[ipar];
02514       }
02515       // add in the contribution from all the backgrounds
02516       binsum +=  bmFit->ND_XsecByParams[i][bmFit->ND_XsecByParams.at(0).size()-1]->GetBinContent(ibin);
02517       bmFit->ND_XSectionFitRew.at(mc)->SetBinContent(ibin,binsum);
02518     }
02519   
02520     
02521     // now calculate the LL contributions (including under and overflows
02522     
02523     // fitting with Maximum Likelihood will revisit this if I need to
02524     //if (!sFit->fUseChi2Function) { 
02525     MAXMSG("NCExtrapolationBeamMatrix",Msg::kInfo,1) << "Using LL test..." << endl;
02526     
02527     // this is for fitting without the bottom two bins; again revisit if I need to...
02528     //if (sFit->fCutLowBins) 
02529     //NCStartBin +=3; // cut out bins 0-2
02530     // get the index right on these things...
02531     int startbin = 1;
02532     int endbin =  bmFit->ND_XSectionFitRew.at(mc)->GetNbinsX()+1 ;
02533     for(int ibin = startbin ; ibin < endbin; ibin++){
02534       if(bmFit->ND_XSectionFitRew.at(mc)->GetBinContent(ibin) < 0 ){
02535         MSG("NCExtrapolationBeamMatrix",Msg::kInfo) << "bmFit->ND_XSectionFitRew.at(mc)->GetBinContent("<<ibin<<") Gone Negative " 
02536                                                     << bmFit->ND_XSectionFitRew.at(mc)->GetBinContent(ibin) << " setting to 1.0 " << endl;
02537       }
02538       f -= bmFit->ND_XSectionFitRew.at(data)->GetBinContent(ibin)
02539         * log(bmFit->ND_XSectionFitRew.at(mc)->GetBinContent(ibin));
02540     }
02541     f +=  bmFit->ND_XSectionFitRew.at(mc)->Integral(startbin,endbin); 
02542   }
02543 //    else {
02544 // fitting with a Chi2 function
02545 //  MAXMSG("NCExtrapolationRS",Msg::kInfo,1) << "Using Chi2 test..." << endl;
02546 //  f += sFit->GetChiSquare(0,i); // add chi2 for NC histos  
02547 //}
02548 // loop over fBeams.size()
02549   
02550   f *= 2; // include factor 2 in order to keep error definitions 
02551   // as in a chi2 function
02552   // this also means, the UP value for errors is 1 (not 0.5)
02553   //reseting the ones that need resetting
02554   //for(UInt_t i = 1; i < bmFit->ND_XSectionFitRew.size(); i+=2){
02555   //bmFit->ND_XSectionFitRew.at(i)->Reset("ICE");
02556   //}
02557 
02558   return;
02559   }

void NCExtrapolationBeamMatrix::makeFluxCorrectionVector ( BeamMatrix bmtouse,
int  beamindex 
) [private]

Definition at line 2726 of file NCExtrapolationBeamMatrix.cxx.

References BMforFitTypes, fluxCorrectionsNDNC, fluxCorrectionsNDNCbkg, BeamMatrix::GetMCHists(), kCCTrueEnergyPurCorrFluxCorr_ND, Msg::kDebug, kNCTrueEnergy_ND, kNCTrueEnergyBkg_ND, kNCTrueEnergyPurCorrFluxCorr_ND, and MSG.

Referenced by correctNDNCFlux().

02726                                                           {
02727 
02728   MSG("NCExtrapolationBeamMatrix", Msg::kDebug) << "MakeFluxCorrectionsVector for " << bmtouse->GetName() << " bmindex " <<  bmindex 
02729                                                 << " normalised to " << BMforFitTypes[0][bmindex]->GetName() << endl;
02730 
02731   for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){
02732     double val = 0.0 ;
02733     if(bmtouse->GetMCHists().at(kNCTrueEnergy_ND)->GetBinContent(i) > 0 )
02734       val = bmtouse->GetMCHists().at(kNCTrueEnergyPurCorrFluxCorr_ND)->GetBinContent(i) / BMforFitTypes[0][bmindex]->GetMCHists().at(kNCTrueEnergy_ND)->GetBinContent(i);
02735     
02736     fluxCorrectionsNDNC[bmindex][i-1] = val ;
02737     
02738     val = 0.0 ;
02739     
02740     if(bmtouse->GetMCHists().at(kNCTrueEnergyBkg_ND)->GetBinContent(i) > 0 )
02741       val = bmtouse->GetMCHists().at(kCCTrueEnergyPurCorrFluxCorr_ND)->GetBinContent(i) / BMforFitTypes[0][bmindex]->GetMCHists().at(kNCTrueEnergyBkg_ND)->GetBinContent(i);
02742     
02743     fluxCorrectionsNDNCbkg[bmindex][i-1] = val ;
02744   }
02745 }

void NCExtrapolationBeamMatrix::normaliseRecoToTrue ( TH2D *  t2rRebin,
TH1D *  norm 
) [private]

Definition at line 2634 of file NCExtrapolationBeamMatrix.cxx.

Referenced by DoOscillations(), and DoUnoscTrueToRecoFD().

02634                                                {
02635 
02636   for(int i = 1; i < kNumEnergyBinsFar+1; ++i ){//reco
02637     for(int j = 1; j < kNumEnergyBinsFar+1; ++j ){//reco
02638       double trueEnergyBinValue = norm->GetBinContent(i) ;
02639       if( trueEnergyBinValue > 0) t2rRebin->SetBinContent(i,j, t2rRebin->GetBinContent(i,j)/trueEnergyBinValue );
02640       else t2rRebin->SetBinContent(i,j,0);
02641       
02642     }
02643   }
02644 }

void NCExtrapolationBeamMatrix::Prepare ( const Registry r  )  [virtual]

Read whatever values you need out of the registry to initialize yourself. Please remember to chain up to the NCExtrapolation implementation too.

Reimplemented from NCExtrapolation.

Definition at line 384 of file NCExtrapolationBeamMatrix.cxx.

References fBeamMatrixPath, fExtractionType, fSelEffFilePath, fSmoothWidth, Registry::Get(), Msg::kInfo, NCType::kNCCuts, logFile, MSG, and SetFiducialMasses().

00385 {
00386   MSG("NCExtrapolationBeamMatrix", Msg::kInfo) <<  "PREPARING BEAM MATRIX EXTRAPOLATION" << endl;
00387 
00388   // the base class takes care of turning the systematic parameters on/off
00389   // and setting their adjusted values if any
00390 
00391   //might need to change this to be passed in as part of the registry
00392   NCType::ECuts cutsType ;
00393   cutsType = NCType::kNCCuts ;
00394   SetFiducialMasses(cutsType) ;
00395     
00396   NCExtrapolation::Prepare(r);  
00397 
00398   int         tmpi;
00399   double      tmpd;
00400   const char* tmps;
00401 
00402   //  if(r.Get("UE3SqrVal",             tmpd)) fUE3Sqr                     = tmpd;
00403 
00404   if (r.Get("ExtractionType", tmpi)) fExtractionType=tmpi;
00405 
00406   // get location of beam matrix helpers
00407   if (r.Get("BeamMatrixFileLocation", tmps)) fBeamMatrixPath=tmps;
00408    // get location of selection and reconstruction efficiency histograms
00409   if (r.Get("SelectionEfficiencyFileLocation", tmps)) fSelEffFilePath=tmps;
00410    if(r.Get("FarNearSmoothingWidth", tmpd)) fSmoothWidth = tmpd;
00411 
00412   //Logfile
00413   r.Get("FileName",tmps);
00414   string outFile(tmps);
00415   string::size_type start = outFile.rfind("/") + 1;
00416   outFile.erase(start,outFile.size());
00417   outFile += "ExtrapolationBeamMatrixFit.log"; 
00418   logFile = new ofstream(outFile.c_str(),ios::out | ios::app);
00419   
00420 }

void NCExtrapolationBeamMatrix::rebinRecoTrueToRecoRecoHist ( TH2D *  t2r,
TH2D *  t2rRebin 
) [private]

Definition at line 2615 of file NCExtrapolationBeamMatrix.cxx.

References kEnergyBinsFar, kNumTrueEnergyBins, and true_bins.

Referenced by DoOscillations(), and DoUnoscTrueToRecoFD().

02615                                                       {
02616 
02617   for(int y = 1; y < kNumEnergyBinsFar+1; ++y ){
02618     int x = 1;
02619     double trueAdd = 0.0 ;
02620     for(UInt_t curx = 1 ; curx <  kNumTrueEnergyBins+1 ; curx++){
02621       if(true_bins[curx] < kEnergyBinsFar[x]){
02622         trueAdd += t2r->GetBinContent(curx,y); 
02623       }      
02624       else{
02625         trueAdd += t2r->GetBinContent(curx,y); 
02626         t2rRebin->SetBinContent(x,y,trueAdd);
02627         trueAdd = 0.;
02628         x++;    
02629       }
02630     }
02631   }
02632 }

void NCExtrapolationBeamMatrix::ScaleBeamsToSameNumberOfEvents ( int  beamType,
bool  doSyst 
) [private]

Definition at line 2561 of file NCExtrapolationBeamMatrix.cxx.

References beamsFit, BMforFit, CopyNDDataForSystFit(), correctNDNCFlux(), FillNDHistsForXSectionFit(), Msg::kInfo, BeamType::kL010z185i, BeamType::kL100z200i, BeamType::kL250z200i, NCType::kNC, mInputScale, MSG, ND_XSectionFit, and size.

Referenced by doNDXsectionFit().

02561                                                          {
02562 
02563   //need to scale all histograms to the same number of events. For now scale them to the HE MC stats? Its what Tobi did.
02564   //I guess it has the lowest stats when you get the full set.
02565   
02566   //even = 0 odd =1
02567   int other = beamType-1 ;
02568   if(beamType%1) other = beamType+1 ;
02569 
02570   MSG("NCExtrapolationBeamMatrix", Msg::kInfo) <<  "Scaling beams to beamType: " << beamType << " " <<  other <<endl;
02571  
02572   for(UInt_t i = 0; i < (beamsFit).size(); ++i){
02573     int bmind = -1 ;
02574     if(beamsFit.at(i)->GetBeamType() == BeamType::kL010z185i ) bmind = 0 ;
02575     if(beamsFit.at(i)->GetBeamType() == BeamType::kL100z200i ) bmind = 1 ;
02576     if(beamsFit.at(i)->GetBeamType() == BeamType::kL250z200i ) bmind = 2 ;
02577     CopyNDDataForSystFit(BMforFit.at(i),doSyst,bmind);
02578     correctNDNCFlux( BMforFit.at(i),bmind);
02579     FillNDHistsForXSectionFit(beamsFit.at(i), NCType::kNC, doSyst);
02580   }
02581   
02582   //can't reweight them all to the same histogram. loses POT weighting for a start.
02583   //need to reweight each data mc pair to the same. in this case to data exposure
02584 
02585  for(UInt_t i = 0; i <  ND_XSectionFit.size() ; i+=2){
02586 
02587    if(ND_XSectionFit.at(i)->Integral() > 0){
02588      double scale =  ND_XSectionFit.at(beamType)->Integral() / ND_XSectionFit.at(i)->Integral();
02589      mInputScale.push_back(scale);
02590      mInputScale.push_back(scale);
02591    }
02592    else{  
02593      mInputScale.push_back(0);
02594      mInputScale.push_back(0);
02595    }
02596    
02597  }
02598  
02599  MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Scaling Beam Parameters  "
02600                                               <<  "mInputScale 0: " << mInputScale.at(0)
02601                                               <<  " mInputScale 1: " << mInputScale.at(1)
02602                                               <<  " mInputScale 2: " << mInputScale.at(2)
02603                                               <<  " mInputScale 3: " << mInputScale.at(3)
02604                                               <<  " mInputScale 4: " << mInputScale.at(4)
02605                                               <<  " mInputScale 5: " << mInputScale.at(5)
02606                                               <<  " mInputScale 6: " << mInputScale.at(6)
02607                                               <<  " mInputScale 7: " << mInputScale.at(7)
02608                                               <<  " mInputScale 8: " << mInputScale.at(8)
02609                                               <<  " mInputScale 9: " << mInputScale.at(9)
02610                                               <<  " mInputScale 10: " << mInputScale.at(10)
02611                                               <<  " mInputScale 11: " << mInputScale.at(11)
02612                                               << endl ;
02613 }

void NCExtrapolationBeamMatrix::SetFiducialMasses ( NCType::ECuts  cutsType  )  [private]

Definition at line 1581 of file NCExtrapolationBeamMatrix.cxx.

References fFDFidVolMass, fFluxPOT, fNDFidVolMass, NCType::kCCCuts, NCType::kNCCCFidCuts, and NCType::kNCCuts.

Referenced by Prepare().

01581                                        {
01582   
01583   fFluxPOT = 1.0;
01584   
01585   switch(cutsType){
01586   case NCType::kCCCuts:
01587     fNDFidVolMass = -9999.9 ;
01588     fFDFidVolMass = -9999.9 ;
01589     break;
01590   case NCType::kNCCuts:
01591     //fNDFidVolMass = 26874.0; //kg cedar era
01592     fNDFidVolMass = 27022.2; //kg
01593     //fFDFidVolMass = 3726961.; //kg cedar era
01594     fFDFidVolMass = 3947780.; //kg
01595     break;
01596   case NCType::kNCCCFidCuts:
01597     fNDFidVolMass = -9999.9 ;
01598     fFDFidVolMass = -9999.9 ;   
01599     break;
01600   default:
01601     fNDFidVolMass = -9999.9 ;
01602     fFDFidVolMass = -9999.9 ;   
01603   }
01604 }

int NCExtrapolationBeamMatrix::WhichFitParam ( double  trueEnergy  )  [private]

Definition at line 2461 of file NCExtrapolationBeamMatrix.cxx.

Referenced by DoOscillations(), FillDataMCHistogramsFar(), FillDataMCHistogramsNear(), and FillNDHistsForXSectionFit().

02461                                 {
02462   if(trueEnergy <=4.0) return 0;
02463   else if(trueEnergy <= 8.0) return 1;
02464   else if(trueEnergy <= 15.0) return 2;
02465   else if(trueEnergy <= 120.0) return 3;
02466   // will need to put errors in at some point
02467   else return 4;
02468 } 

int NCExtrapolationBeamMatrix::WhichFluxCorrParam ( double  trueEnergy  )  [private]

Definition at line 2747 of file NCExtrapolationBeamMatrix.cxx.

References kEnergyBinsFar.

Referenced by FillNDHistsForXSectionFit().

02747                                      {
02748   
02749   for(int i = 0 ; i <  kNumEnergyBinsFar ; i++){
02750     if(trueEnergy <= kEnergyBinsFar[i+1]) return i;
02751   }
02752   return -1 ;
02753 } 

void NCExtrapolationBeamMatrix::WriteResources ( const NC::OscProb::OscPars trueOscPars  )  [virtual]

gDirectory will point to the output file. Please remember to chain up to the NCExtrapolation implementation.

trueOscPars may be zero (eg obviously in case of fit to real data)

Reimplemented from NCExtrapolation.

Definition at line 1046 of file NCExtrapolationBeamMatrix.cxx.

References BMbybeams, BMforFit, ConstructFarSpectrum(), NCExtrapolation::fBeams, fBMbybeams, NCExtrapolation::fCoordConv, NCExtrapolation::fDoSystematicInterpolation, FillResultHistograms(), NCExtrapolation::fInterpolator, fNCalls, fNominalBM, Form(), fOscBM, NCExtrapolation::fUseCC, NCExtrapolation::fUseNC, NCBeam::GetBeamType(), NCExtrapolation::GetBestFitOscPars(), NCExtrapolation::GetBestFitSysts(), NCBeam::Info::GetDescription(), NCBeam::GetDetector(), NCBeam::GetInfo(), NC::ISpectrumInterpolator::GetInterpolatedSpectra(), handler(), IgnoreErr(), it, NCType::kCC, kCCSig, Detector::kFar, Msg::kInfo, BeamType::kL010z185i, NCType::kNC, kNCSig, MSG, NC::Utility::MultiplyFast(), ND_XSectionFit, ND_XSectionFitRew, Registry::Set(), Registry::UnLockKeys(), Registry::UnLockValues(), and NC::CoordinateConverter::VectorFromSystPars().

01047 {
01048   MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "nCalls: " << fNCalls << endl;
01049 
01050   NC::SystPars systPars = GetBestFitSysts();
01051   const NC::OscProb::OscPars* oscPars = GetBestFitOscPars();
01052   
01053   for(fBeams_t::iterator it=fBeams.begin(); it!=fBeams.end(); ++it){
01054     NCBeam::Info bmInfo=it->first.first;
01055     // it->first.second is the detector. We'll just do the FD for now, so ignore it
01056     NCBeam* beam=it->second;
01057     
01058     int bmequiv = fBMbybeams.find(bmInfo)->second;
01059 
01060     //only le beams are filled only need to do one of them.
01061     if(beam->GetBeamType() !=  BeamType::kL010z185i) continue ;
01062     if(beam->GetDetector() !=  Detector::kFar) continue ;
01063     
01064     MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Write Resources for " << bmInfo.GetDescription() << " bmequiv " << bmequiv  << endl;
01065     
01066     if(fDoSystematicInterpolation){
01067       const vector<double> shift = fCoordConv.VectorFromSystPars(systPars);
01068       vector<TH1*> interp = fInterpolator->GetInterpolatedSpectra(shift);
01069       
01070       ConstructFarSpectrum(oscPars, BMbybeams[bmequiv]);    
01071       
01072       MultiplyFast(BMbybeams[bmequiv]->GetFitHists()[kNCSig], interp[0]);
01073       MultiplyFast(BMbybeams[bmequiv]->GetFitHists()[kCCSig], interp[2]);
01074       
01075     }
01076     
01077     else ConstructFarSpectrum(oscPars, BMbybeams[bmequiv]);
01078     
01079     
01080     if(fUseNC) FillResultHistograms(bmInfo, NCType::kNC, fNominalBM, fOscBM);
01081     if(fUseCC) FillResultHistograms(bmInfo, NCType::kCC, fNominalBM, fOscBM);  
01082     
01083     
01084     // get a pointer to the current directory
01085     // this is one of the output files
01086     
01087     TDirectory* saveDir = gDirectory;
01088     
01089     MSG("NCExtrapolationBeamMatrix", Msg::kInfo) << "Writing out for " <<  beam->GetInfo().GetDescription()  << endl;
01090     //writing out histograms. Any others wanted can be easily added ME or HE beam for example
01091     for(unsigned int i = 0; i < 2; ++i) BMbybeams[bmequiv]->GetFDDataHists()[i]->Write();
01092     for(unsigned int i = 0; i < BMbybeams[bmequiv]->GetMCHists().size(); ++i) BMbybeams[bmequiv]->GetMCHists().at(i)->Write();
01093     for(unsigned int i = 0; i < BMbybeams[bmequiv]->GetDataPredHists().size(); ++i) BMbybeams[bmequiv]->GetDataPredHists().at(i)->Write();
01094     for(unsigned int i = 0; i < BMbybeams[bmequiv]->Get2DHists().size(); ++i)  BMbybeams[bmequiv]->Get2DHists().at(i)->Write();
01095     for(unsigned int i = 0; i < BMbybeams[bmequiv]->GetFitHists().size();  ++i)  BMbybeams[bmequiv]->GetFitHists().at(i)->Write();
01096     for(unsigned int i = 0; i < ND_XSectionFitRew.size();  ++i)  ND_XSectionFitRew[i]->Write();
01097     for(unsigned int i = 0; i < ND_XSectionFit.size();  ++i)  ND_XSectionFit[i]->Write();
01098     for(unsigned int i = 0; i < BMbybeams[bmequiv]->GetEffCorrHists().size();  ++i)  BMbybeams[bmequiv]->GetEffCorrHists().at(i)->Write();
01099     
01100     BMbybeams[bmequiv]->GetBeamMatrix()->Write();
01101    
01102     //writing the result of the x-section fit into the file
01103     TString parname ;
01104     TString parout ;
01105     
01106     //stolen from NCExtrapolation
01107     ErrorHandlerFunc_t handler = GetErrorHandler();
01108     SetErrorHandler(IgnoreErr);
01109     if(!gDirectory->cd("xSecFit")) gDirectory->mkdir("xSecFit","xSecFit")->cd();
01110     SetErrorHandler(handler);
01111 
01112     Registry* xsecFit = new Registry;
01113     TString bmname = BMbybeams[bmequiv]->GetName();
01114     TString xsecname = "xSecFit_registry"+bmname;
01115     xsecFit->SetName(xsecname);
01116     xsecFit->UnLockKeys();
01117     xsecFit->UnLockValues();
01118     for (Int_t i=0; i<4; i++){ 
01119       parname = (TString) Form("Par %i: ",i);
01120       //parout = (TString) Form("%e +%e %e +/- %e, Corr: %e",mScaleNC[i],mScalePos[i],mScaleNeg[i],mScaleNCErr[i],corr[i]);
01121       parout = (TString) Form("%e +%e %e +/- %e, Corr: %e",BMbybeams[bmequiv]->GetmScaleNC().at(i),0.,0.,0.,0.);
01122       xsecFit->Set(parname,parout);
01123     }
01124     xsecFit->Write();
01125     saveDir->cd();
01126     delete xsecFit;
01127 
01128   }
01129   
01130   for(unsigned int i = 0; i < 40 ; ++i) BMforFit.at(2)->GetMCHists().at(i)->Write();
01131   for(unsigned int i = 0; i < 40 ; ++i) BMforFit.at(1)->GetMCHists().at(i)->Write();
01132   for(unsigned int i = 0; i < 40 ; ++i) BMforFit.at(0)->GetMCHists().at(i)->Write();
01133   for(unsigned int i = 0; i < 11 ; ++i) BMforFit.at(2)->GetDataPredHists().at(i)->Write();
01134   for(unsigned int i = 0; i < 11 ; ++i) BMforFit.at(1)->GetDataPredHists().at(i)->Write();
01135   for(unsigned int i = 0; i < 11 ; ++i) BMforFit.at(0)->GetDataPredHists().at(i)->Write();
01136   for(unsigned int i = 0; i < BMforFit.at(2)->GetEffCorrHists().size();  ++i)  BMforFit[2]->GetEffCorrHists().at(i)->Write();
01137   for(unsigned int i = 0; i < BMforFit.at(1)->GetEffCorrHists().size();  ++i)  BMforFit[1]->GetEffCorrHists().at(i)->Write();
01138   for(unsigned int i = 0; i < BMforFit.at(0)->GetEffCorrHists().size();  ++i)  BMforFit[0]->GetEffCorrHists().at(i)->Write();
01139 
01140   //for(unsigned int i = 0; i < 5 ;  ++i){
01141   //ND_XsecByParams[0][i]->Write();
01142   //ND_XsecByParams[1][i]->Write();
01143   //ND_XsecByParams[2][i]->Write();
01144   //}
01145 
01146   delete oscPars;
01147   NCExtrapolation::WriteResources(trueOscPars);    
01148 
01149 }


Member Data Documentation

std::vector< NCBeam* > NCExtrapolationBeamMatrix::beamsFit [private]
std::vector< NCBeam* > NCExtrapolationBeamMatrix::beamsNear [private]

Definition at line 186 of file NCExtrapolationBeamMatrix.h.

Referenced by doNDXsectionFit(), and DoneFilling().

std::vector<std::vector< BeamMatrix* > > NCExtrapolationBeamMatrix::BMforFitTypes [private]

Definition at line 196 of file NCExtrapolationBeamMatrix.h.

Referenced by doNDXsectionFit(), and makeFluxCorrectionVector().

double NCExtrapolationBeamMatrix::corr[4] [private]

Definition at line 215 of file NCExtrapolationBeamMatrix.h.

Referenced by doNDXsectionFit(), and Fit().

std::vector<std::vector<TH1D* > > NCExtrapolationBeamMatrix::dataNDhist [private]

Definition at line 107 of file NCExtrapolationBeamMatrix.h.

Definition at line 105 of file NCExtrapolationBeamMatrix.h.

Definition at line 106 of file NCExtrapolationBeamMatrix.h.

Definition at line 110 of file NCExtrapolationBeamMatrix.h.

Definition at line 108 of file NCExtrapolationBeamMatrix.h.

Definition at line 109 of file NCExtrapolationBeamMatrix.h.

Definition at line 121 of file NCExtrapolationBeamMatrix.h.

Definition at line 120 of file NCExtrapolationBeamMatrix.h.

Referenced by GetBeamMatrix(), and Prepare().

Definition at line 219 of file NCExtrapolationBeamMatrix.h.

Referenced by GetEfficiencyHistograms(), and Prepare().

Definition at line 223 of file NCExtrapolationBeamMatrix.h.

Definition at line 54 of file NCExtrapolationBeamMatrix.h.

Referenced by DoXsectionCorrectionFD(), and SetFiducialMasses().

Definition at line 55 of file NCExtrapolationBeamMatrix.h.

Referenced by DoneFilling(), and DoXsectionCorrectionFD().

Definition at line 56 of file NCExtrapolationBeamMatrix.h.

Referenced by DoneFilling(), and DoXsectionCorrectionFD().

std::vector<Double_t > NCExtrapolationBeamMatrix::fluxCorrectionsFDCCinNC [private]
std::vector<std::vector<Double_t > > NCExtrapolationBeamMatrix::fluxCorrectionsNDNC [private]
std::vector<std::vector<Double_t > > NCExtrapolationBeamMatrix::fluxCorrectionsNDNCbkg [private]

Counts number of calls to the method.

Definition at line 139 of file NCExtrapolationBeamMatrix.h.

Referenced by NCExtrapolationBeamMatrix(), and WriteResources().

Definition at line 54 of file NCExtrapolationBeamMatrix.h.

Referenced by DoXsectionCorrectionND(), and SetFiducialMasses().

Definition at line 55 of file NCExtrapolationBeamMatrix.h.

Referenced by DoneFilling(), and DoXsectionCorrectionND().

Definition at line 56 of file NCExtrapolationBeamMatrix.h.

Referenced by DoneFilling(), and DoXsectionCorrectionND().

std::vector<std::vector<double> > NCExtrapolationBeamMatrix::fNominal [private]

Definition at line 141 of file NCExtrapolationBeamMatrix.h.

std::vector<std::vector<double> > NCExtrapolationBeamMatrix::fNominalBM [private]

Definition at line 144 of file NCExtrapolationBeamMatrix.h.

Referenced by DoOscillations(), and WriteResources().

std::vector<std::vector<double> > NCExtrapolationBeamMatrix::fOsc [private]

Definition at line 142 of file NCExtrapolationBeamMatrix.h.

std::vector<std::vector<double> > NCExtrapolationBeamMatrix::fOscBM [private]

Definition at line 145 of file NCExtrapolationBeamMatrix.h.

Referenced by DoOscillations(), and WriteResources().

Definition at line 96 of file NCExtrapolationBeamMatrix.h.

Referenced by GetEfficiencyHistograms(), and Prepare().

Definition at line 147 of file NCExtrapolationBeamMatrix.h.

Definition at line 152 of file NCExtrapolationBeamMatrix.h.

std::ofstream* NCExtrapolationBeamMatrix::logFile [private]

Definition at line 217 of file NCExtrapolationBeamMatrix.h.

Referenced by Fit(), and Prepare().

std::vector<Double_t> NCExtrapolationBeamMatrix::m0ScaleNC [private]

Definition at line 203 of file NCExtrapolationBeamMatrix.h.

Referenced by doNDXsectionFit(), DoneFilling(), and Fit().

std::vector<Double_t> NCExtrapolationBeamMatrix::m0ScaleNCErr [private]

Definition at line 204 of file NCExtrapolationBeamMatrix.h.

Referenced by doNDXsectionFit(), DoneFilling(), and Fit().

std::vector<Double_t> NCExtrapolationBeamMatrix::mInputScale [private]

Definition at line 199 of file NCExtrapolationBeamMatrix.h.

Referenced by doNDXsectionFit(), and ScaleBeamsToSameNumberOfEvents().

std::vector<Double_t> NCExtrapolationBeamMatrix::mScaleNC [private]

Definition at line 200 of file NCExtrapolationBeamMatrix.h.

Referenced by doNDXsectionFit(), DoneFilling(), and Fit().

std::vector<Double_t> NCExtrapolationBeamMatrix::mScaleNCErr [private]

Definition at line 201 of file NCExtrapolationBeamMatrix.h.

Referenced by doNDXsectionFit(), DoneFilling(), and Fit().

Definition at line 212 of file NCExtrapolationBeamMatrix.h.

Referenced by doNDXsectionFit(), and Fit().

Definition at line 214 of file NCExtrapolationBeamMatrix.h.

Referenced by doNDXsectionFit(), and Fit().

std::vector<std::vector<TH1D* > > NCExtrapolationBeamMatrix::ND_XsecByParams [private]
std::vector< TH1D* > NCExtrapolationBeamMatrix::ND_XSectionFit [private]
std::vector< TH1D* > NCExtrapolationBeamMatrix::ND_XSectionFitRew [private]

Definition at line 112 of file NCExtrapolationBeamMatrix.h.

Referenced by DoRecoToTrueND(), and NCExtrapolationBeamMatrix().

Definition at line 113 of file NCExtrapolationBeamMatrix.h.

Referenced by NCExtrapolationBeamMatrix().

Definition at line 221 of file NCExtrapolationBeamMatrix.h.

std::vector<double> NCExtrapolationBeamMatrix::true_bins [private]

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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1