NCExtrapolationFarNear Class Reference

Extrapolation using the Far over Near method. More...

#include <NCExtrapolationFarNear.h>

Inheritance diagram for NCExtrapolationFarNear:
NCExtrapolation

List of all members.

Public Member Functions

 NCExtrapolationFarNear ()
virtual ~NCExtrapolationFarNear ()
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.
virtual 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 Member Functions

void ConstructFarSpectrum (const NC::OscProb::OscPars *pars)
void FillDataMCHistogramsNear (NCBeam *beam, NCType::EEventType nccc, NC::SystPars systPars)
void FillDataMCHistogramsFar (NCBeam *beam, NCType::EEventType nccc, NC::SystPars systPars)
void FillResultHistograms (NCBeam *beam, NCType::EEventType nccc, std::vector< std::vector< double > > &nominal, std::vector< std::vector< double > > &osc)
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)
 This implementation ignores beamsToUse.
virtual void CleanupSpectra (std::vector< TH1 * > exps, std::vector< TH1 * > obss)
 Called after FindSpectraForPars() to delete necessary spectra.
void GetNearMCSpectra (NCBeam *beam, NC::SystPars systs, TH1 **nc, TH1 **cc)

Private Attributes

int fNCalls
 Counts number of calls to the method.
std::vector< std::vector
< double > > 
fNominal
std::vector< std::vector
< double > > 
fOsc
double fTrue_bin_width
double fUE3Sqr
std::vector< TH2D * > fNom_true_vs_reco
 Histograms for reco fitting work.
std::vector< TH1D * > fFN
 Far over Near histogram.
std::vector< TH1D * > fFD_fit
std::vector< TH1D * > fND_data
std::vector< TH1D * > fND_MC
TH1D * fFD_data [2]
double fSmoothWidth
bool fDoExtrapolation

Detailed Description

Extrapolation using the Far over Near method.

Definition at line 31 of file NCExtrapolationFarNear.h.


Constructor & Destructor Documentation

NCExtrapolationFarNear::NCExtrapolationFarNear (  ) 

Definition at line 51 of file NCExtrapolationFarNear.cxx.

References kEnergyBinsFar, kNumEnergyBinsFar, and kNumTrueEnergyBins.

00052 {
00053   fNCalls = 0;
00054   fDoExtrapolation = true;
00055 
00056   int numBins = kNumEnergyBinsFar;
00057   double bins[kNumEnergyBinsFar+1];
00058 
00059   //use NCType to fill array of bins as it has the right binning already
00060   for( int i = 0; i < numBins+1; ++i ){
00061     bins[i] = kEnergyBinsFar[i];
00062   }
00063 
00064   // Create a true energy binning scheme for the
00065   // F/N matrix
00066 
00067   fTrue_bin_width = bins[numBins]/double(kNumTrueEnergyBins);
00068   vector<double> true_bins(kNumTrueEnergyBins+1, 0);
00069 
00070   for( int i = 0; i < kNumTrueEnergyBins+1; ++i )
00071     true_bins[i] = i*fTrue_bin_width;
00072 
00073 
00074   //--------------------------------------------------
00075   // Far Detector Histograms
00076 
00077   fNom_true_vs_reco.push_back(new TH2D ( "NCNomtruevsreco", "NC True vs. Reco (Nom)",
00078                                          numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
00079   fNom_true_vs_reco.push_back(new TH2D ( "CCNomtruevsreco", "CC True vs. Reco (Nom)",
00080                                          numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
00081   fNom_true_vs_reco.push_back(new TH2D ( "NCNomtruevsreco_bkg", "NC True vs. Reco (Nom)bkg",
00082                                          numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
00083   fNom_true_vs_reco.push_back(new TH2D ( "CCNomtruevsreco_bkg", "CC True vs. Reco (Nom)bkg",
00084                                          numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
00085   fNom_true_vs_reco.push_back(new TH2D ( "NCNomtruevsreco_tau", "NC True vs. Reco (Nom)tau",
00086                                          numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
00087   fNom_true_vs_reco.push_back(new TH2D ( "CCNomtruevsreco_tau", "CC True vs. Reco (Nom)tau",
00088                                          numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
00089   fNom_true_vs_reco.push_back(new TH2D ( "NCNomtruevsreco_beamnue", "NC True vs. Reco (Nom)beamnue",
00090                                          numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
00091   fNom_true_vs_reco.push_back(new TH2D ( "CCNomtruevsreco_beamnue", "CC True vs. Reco (Nom)beamnue",
00092                                          numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
00093   fNom_true_vs_reco.push_back(new TH2D ( "NCNomtruevsreco_oscnue", "NC True vs. Reco (Nom)oscnue",
00094                                          numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
00095   fNom_true_vs_reco.push_back(new TH2D ( "CCNomtruevsreco_oscnue", "CC True vs. Reco (Nom)oscnue",
00096                                          numBins, bins, kNumTrueEnergyBins, &true_bins.at(0) ));
00097 
00098   fFN.push_back( new TH1D( "FNNC" , "FN NC", numBins, bins ) );
00099   fFN.push_back( new TH1D( "FNCC" , "FN CC", numBins, bins ) );
00100 
00101   fFD_fit.push_back( new TH1D( "FD_fit_nc" , "FD_fit NC", numBins, bins ) );
00102   fFD_fit.push_back( new TH1D( "FD_fit_cc" , "FD_fit CC", numBins, bins ) );
00103 
00104   //--------------------------------------------------
00105   // Near Detector Histograms
00106 
00107   fND_data.push_back( new TH1D( "ND_data_NC" , "ND_data_NC", numBins, bins ) );
00108   fND_data.push_back( new TH1D( "ND_data_CC" , "ND_data_CC", numBins, bins ) );
00109 
00110   //------------------------------
00111   // muon neutrinos
00112 
00113   fND_MC.push_back( new TH1D( "ND_MC_NC" , "ND_MC_NC", numBins, bins ) );
00114   fND_MC.push_back( new TH1D( "ND_MC_CC" , "ND_MC_CC", numBins, bins ) );
00115   fND_MC.push_back( new TH1D( "ND_MC_NC_ccbkg" , "ND_MC_NC-ccbkg", numBins, bins ) );
00116   fND_MC.push_back( new TH1D( "ND_MC_CC_ncbkg" , "ND_MC_CC-ncbkg", numBins, bins ) );
00117 
00118   fND_MC.push_back( new TH1D( "ND_MC_NC_tau" , "ND_MC_NC-tau", numBins, bins ) );
00119   fND_MC.push_back( new TH1D( "ND_MC_CC_tau" , "ND_MC_CC-tau", numBins, bins ) );
00120 
00121   fND_MC.push_back( new TH1D( "ND_MC_NC_beamnue" , "ND_MC_NC-beamnue", numBins, bins ) );
00122   fND_MC.push_back( new TH1D( "ND_MC_CC_beamnue" , "ND_MC_CC-beamnue", numBins, bins ) );
00123   fND_MC.push_back( new TH1D( "ND_MC_NC_oscnue" , "ND_MC_NC-oscnue", numBins, bins ) );
00124   fND_MC.push_back( new TH1D( "ND_MC_CC_oscnue" , "ND_MC_CC-oscnue", numBins, bins ) );
00125 
00126   fFD_data[kNCSig] = new TH1D("FD_data_NC", "FD data NC", kNumEnergyBinsFar, kEnergyBinsFar);
00127   fFD_data[kCCSig] = new TH1D("FD_data_CC", "FD data CC", kNumEnergyBinsFar, kEnergyBinsFar);
00128 
00129   fSmoothWidth = 0; // Don't do any smoothing
00130 }

NCExtrapolationFarNear::~NCExtrapolationFarNear (  )  [virtual]

Definition at line 136 of file NCExtrapolationFarNear.cxx.

00136 {}


Member Function Documentation

virtual void NCExtrapolationFarNear::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 NCExtrapolationFarNear::ConstructFarSpectrum ( const NC::OscProb::OscPars pars  )  [private]

Definition at line 463 of file NCExtrapolationFarNear.cxx.

References fDoExtrapolation, fFD_fit, fFN, fND_data, fND_MC, fNom_true_vs_reco, fNominal, fOsc, fTrue_bin_width, NCExtrapolation::fUseCC, NCExtrapolation::fUseNC, NC::OscProb::OscPars::IsEquiv(), NCType::kBaseLineFar, NCType::kCC, kCCBeamNue, kCCBkg, kCCOscNue, kCCSig, kCCTau, Msg::kError, Msg::kInfo, NCType::kNC, kNCBeamNue, kNCBkg, kNCOscNue, kNCSig, kNCTau, NCType::kNuEToNuE, kNumEnergyBinsFar, kNumTrueEnergyBins, NCType::kNuMuToNuE, NCType::kNuMuToNuMu, NCType::kNuMuToNuTau, MSG, n, and NC::OscProb::OscPars::TransitionProbability().

Referenced by WriteResources().

00464 {
00465   static vector<double*> survivalProbability;
00466 
00467   static bool first = true;
00468   if(first){
00469     first = false;
00470 
00471     fNominal.resize(fNom_true_vs_reco.size());
00472     fOsc.resize(fNom_true_vs_reco.size());
00473 
00474     for(unsigned int h = 0; h < fNom_true_vs_reco.size(); ++h){
00475       fNominal[h].resize(kNumEnergyBinsFar);
00476       fOsc[h].resize(kNumEnergyBinsFar);
00477     }
00478 
00479     survivalProbability.resize(fNom_true_vs_reco.size());
00480     for(unsigned int n = 0; n < fNom_true_vs_reco.size(); ++n)
00481       survivalProbability[n] = new double[kNumTrueEnergyBins];
00482   }
00483 
00484   static const NC::OscProb::OscPars* prevpars = 0;
00485 
00486   // Transition probability does not depend on reco energy - so we
00487   // can precalculate it out here.
00488   // Reuse the calculations from last time round if the oscillations
00489   // parameters haven't changed since then.
00490   if(!prevpars || !pars->IsEquiv(prevpars)){
00491     for(int h = 0; h < int(fNom_true_vs_reco.size()); ++h){
00492       NCType::EOscMode oscType;
00493       NCType::EEventType interaction;
00494 
00495       switch(h){
00496       case kNCSig:
00497       case kCCBkg:
00498         interaction = NCType::kNC;
00499         oscType     = NCType::kNuMuToNuMu;
00500         break;
00501       case kCCSig:
00502       case kNCBkg:
00503         interaction = NCType::kCC;
00504         oscType     = NCType::kNuMuToNuMu;
00505         break;
00506       case kNCTau:
00507       case kCCTau:
00508         interaction = NCType::kCC;
00509         oscType     = NCType::kNuMuToNuTau;
00510         break;
00511       case kNCBeamNue:
00512       case kCCBeamNue:
00513         interaction = NCType::kCC;
00514         oscType     = NCType::kNuEToNuE;
00515         break;
00516       case kNCOscNue:
00517       case kCCOscNue:
00518         interaction = NCType::kCC;
00519         oscType     = NCType::kNuMuToNuE;
00520         break;
00521       default:
00522         MSG("NCExtrapolationFarNear", Msg::kError) << "Interaction type or flavor is incorrect" << endl;
00523         return;
00524       }
00525 
00526       double* sph = survivalProbability[h];
00527       for(int j = 0; j < kNumTrueEnergyBins; ++j){
00528         sph[j] = pars->TransitionProbability(oscType,
00529                                              interaction,
00530                                              NCType::kBaseLineFar,
00531                                              (j+0.5)*fTrue_bin_width);
00532       } // end for j
00533     } // end for h
00534   } // end if
00535 
00536 
00537   delete prevpars;
00538   // Sadly OscPars::Clone isn't labelled const for technical reasons.
00539   // Need to cast away constness so that we can call it.
00540   prevpars = const_cast<NC::OscProb::OscPars*>(pars)->Clone();
00541 
00542 
00543   const int nx = kNumEnergyBinsFar+2;
00544   for(unsigned int h = 0; h < fNom_true_vs_reco.size(); ++h){
00545     const double* t2ra = fNom_true_vs_reco[h]->fArray;
00546     vector<double>& nh = fNominal[h];
00547     vector<double>& oh = fOsc[h];
00548     for(int i = 0; i < kNumEnergyBinsFar; ++i){
00549       double& nhi = nh[i];
00550       double& ohi = oh[i];
00551       nhi = 0;
00552       ohi = 0;
00553       const double* sph = survivalProbability[h];
00554       const int oset = i+1+nx;
00555       for(int j = 0; j < kNumTrueEnergyBins; ++j){
00556         //      const double tmpd = t2r->GetCellContent(i+1, j+1);
00557         const double tmpd = t2ra[nx*j+oset];
00558         nhi += tmpd;
00559         ohi += sph[j]*tmpd;
00560       }//end loop over true bins
00561     } // end for i
00562   }//end loop over histograms
00563 
00564   for(int i = 0; i < kNumEnergyBinsFar; ++i){
00565     const double totalMCNC_fd = fOsc[kNCSig][i]
00566       + fOsc[kNCBkg][i]
00567       + fOsc[kNCTau][i]
00568       + fOsc[kNCBeamNue][i]
00569       + fOsc[kNCOscNue][i];
00570 
00571     const double totalMCCC_fd = fOsc[kCCSig][i]
00572       + fOsc[kCCBkg][i]
00573       + fOsc[kCCTau][i]
00574       + fOsc[kCCBeamNue][i]
00575       + fOsc[kCCOscNue][i];
00576 
00577     // TODO - implement this bit in terms of GetNearMCSpectra. Needs some
00578     // work to keep the cacheing correct.
00579     const double totalMCNC_nd = fND_MC.at(kNCSig)->GetBinContent(i+1)
00580       + fND_MC.at(kNCBkg)->GetBinContent(i+1)      // numu
00581       + fND_MC.at(kNCTau)->GetBinContent(i+1)      // nutau
00582       + fND_MC.at(kNCBeamNue)->GetBinContent(i+1)  // beamnue
00583       + fND_MC.at(kNCOscNue)->GetBinContent(i+1);  // oscnue
00584 
00585     const double totalMCCC_nd = fND_MC.at(kCCSig)->GetBinContent(i+1)
00586       + fND_MC.at(kCCBkg)->GetBinContent(i+1)      // numu
00587       + fND_MC.at(kCCTau)->GetBinContent(i+1)      // nutau
00588       + fND_MC.at(kCCBeamNue)->GetBinContent(i+1)  // beamnue
00589       + fND_MC.at(kCCOscNue)->GetBinContent(i+1);  // oscnue
00590 
00591     if(totalMCNC_nd > 0){
00592       fFN[kNCSig]->SetBinContent(i+1, totalMCNC_fd/totalMCNC_nd);
00593       // Set this spectrum up without F/N correction. Will be overwritten if
00594       // fDoExtrapolation. If not then we need it done this way.
00595       fFD_fit[kNCSig]->SetBinContent(i+1, totalMCNC_fd);
00596     }
00597     else if(totalMCNC_nd == 0 && fUseNC){
00598       MSG("NCExtrapolationFarNear", Msg::kInfo) << "MC ND NC HISTO "
00599                                                << "HAS ZERO BIN:"
00600                                                << i+1 << ". PLEASE FIX IT"
00601                                                << endl;
00602     }
00603 
00604     if(totalMCCC_nd > 0){
00605       fFN[kCCSig]->SetBinContent(i+1, totalMCCC_fd/totalMCCC_nd);
00606       // Set this spectrum up without F/N correction. Will be overwritten if
00607       // fDoExtrapolation. If not then we need it done this way.
00608       fFD_fit[kCCSig]->SetBinContent(i+1, totalMCCC_fd);
00609     }
00610     else if(totalMCCC_nd == 0. && fUseCC){
00611       MSG("NCExtrapolationFarNear", Msg::kInfo) << "MC ND CC HISTO "
00612                                                << "HAS ZERO BIN:"
00613                                                << i+1 << ". PLEASE FIX IT"
00614                                                << endl;
00615     }
00616   }//end loop over reco energy bins
00617 
00618   if(fDoExtrapolation){
00619     fFD_fit[kNCSig]->Multiply(fND_data.at(kNCSig), fFN.at(kNCSig));
00620     fFD_fit[kCCSig]->Multiply(fND_data.at(kCCSig), fFN.at(kCCSig));
00621   }
00622 }

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

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

Reimplemented from NCExtrapolation.

Definition at line 139 of file NCExtrapolationFarNear.cxx.

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

00140 {
00141   static Registry r;
00142 
00143   r.UnLockValues();
00144 
00145   r.Set("FarNearSmoothingWidth", 0.05);
00146 
00147   r.LockValues();
00148   return r;
00149 }

bool NCExtrapolationFarNear::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 975 of file NCExtrapolationFarNear.cxx.

References fDoExtrapolation.

00976 {
00977   const bool ret = fDoExtrapolation;
00978   fDoExtrapolation = enable;
00979   return ret;
00980 }

void NCExtrapolationFarNear::FillDataMCHistogramsFar ( NCBeam beam,
NCType::EEventType  nccc,
NC::SystPars  systPars 
) [private]

Definition at line 368 of file NCExtrapolationFarNear.cxx.

References NC::SystPars::CCBkgScale(), MuELoss::e, fFD_data, FillHistogramSmoothed2D(), fNom_true_vs_reco, 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(), NCType::kCC, kCCBeamNue, kCCBkg, kCCOscNue, kCCSig, kCCTau, kFar, Detector::kFar, NCType::kNC, kNCBeamNue, kNCBkg, kNCOscNue, kNCSig, kNCTau, n, NC::SystPars::NCBkgScale(), NC::SystPars::NormScale(), NC::SystPars::ShowerScale(), and NC::SystPars::TrackScale().

Referenced by WriteResources().

00371 {
00372   using Detector::kFar;
00373 
00374   //loop over NCEnergy bin objects in the near detector
00375   //to get the energies from the data events in each bin
00376   //and fill the fFD_data histogram
00377   double trueEnergy    = 0.;
00378   double recoShowerE   = 0.;
00379   double recoMuonE     = 0.;
00380   double weight        = 0.;
00381   int    flavor        = 0;
00382 
00383   int sig         = kNCSig;
00384   int bkg         = kNCBkg;
00385   int tau         = kNCTau;
00386   int beamnue     = kNCBeamNue;
00387   int oscnue      = kNCOscNue;
00388   double bkgShift = systPars.CCBkgScale();
00389 
00390   if(nccc == NCType::kCC){
00391     sig      = kCCSig;
00392     bkg      = kCCBkg;
00393     tau      = kCCTau;
00394     beamnue  = kCCBeamNue;
00395     oscnue   = kCCOscNue;
00396     bkgShift = systPars.NCBkgScale();
00397   }
00398 
00399   fNom_true_vs_reco[sig]->Reset( "ICE" );
00400   fNom_true_vs_reco[bkg]->Reset( "ICE" );
00401   fNom_true_vs_reco[tau]->Reset( "ICE" );
00402   fNom_true_vs_reco[beamnue]->Reset( "ICE" );
00403   fNom_true_vs_reco[oscnue]->Reset( "ICE" );
00404 
00405 
00406   const int B = beam->GetNumberEnergyBins(nccc);
00407   for(int b = 0; b < B; ++b){
00408     const NCEnergyBin* bin = beam->GetEnergyBin(b, nccc);
00409 
00410     const int mcSize         = bin->GetMCSignalVectorSize();
00411     const int mcSize_bkg     = bin->GetMCBackgroundVectorSize();
00412     const int mcSize_tau     = bin->GetMCNuTauVectorSize();
00413     const int mcSize_beamnue = bin->GetMCBeamNuEVectorSize();
00414     const int mcSize_oscnue  = bin->GetMCOscNuEVectorSize();
00415 
00416     //numu
00417     for(int e = 0; e < mcSize; ++e){
00418       bin->GetMCInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00419       if(nccc == NCType::kNC) recoMuonE = 0;
00420       FillHistogramSmoothed2D(fNom_true_vs_reco[sig], recoShowerE*systPars.ShowerScale(kFar) + recoMuonE*systPars.TrackScale(), trueEnergy, weight*systPars.NormScale());
00421     }
00422 
00423     for(int e = 0; e < mcSize_bkg; ++e){
00424       bin->GetMCBackgroundInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00425       if(nccc == NCType::kNC) recoMuonE = 0;
00426       FillHistogramSmoothed2D(fNom_true_vs_reco[bkg], recoShowerE*systPars.ShowerScale(kFar) + recoMuonE*systPars.TrackScale(), trueEnergy, weight*systPars.NormScale()*bkgShift);
00427     }
00428 
00429     //nutau
00430     for(int e = 0; e < mcSize_tau; ++e){
00431       bin->GetMCNuTauInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00432       if(nccc == NCType::kNC) recoMuonE = 0;
00433       FillHistogramSmoothed2D(fNom_true_vs_reco[tau], recoShowerE*systPars.ShowerScale(kFar) + recoMuonE*systPars.TrackScale(), trueEnergy, weight*systPars.NormScale());
00434     }
00435 
00436     //beamnue
00437     for(int e = 0; e < mcSize_beamnue; ++e){
00438       bin->GetMCBeamNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00439       if(nccc == NCType::kNC) recoMuonE = 0;
00440       FillHistogramSmoothed2D(fNom_true_vs_reco[beamnue], recoShowerE*systPars.ShowerScale(kFar) + recoMuonE*systPars.TrackScale(), trueEnergy, weight*systPars.NormScale());
00441     }
00442 
00443     //oscnue
00444     for(int e = 0; e < mcSize_oscnue; ++e){
00445       bin->GetMCOscNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00446       if(nccc == NCType::kNC) recoMuonE = 0;
00447       FillHistogramSmoothed2D(fNom_true_vs_reco[oscnue], recoShowerE*systPars.ShowerScale(kFar) + recoMuonE*systPars.TrackScale(), trueEnergy, weight*systPars.NormScale());
00448     }
00449 
00450     const int N = bin->GetDataVectorSize();
00451     for(int n = 0; n < N; ++n){
00452       double energy, weight;
00453       bin->GetDataInformation(energy, weight, n);
00454       fFD_data[sig]->Fill(energy, weight);
00455     } // end for n
00456   } // end for b
00457 }

void NCExtrapolationFarNear::FillDataMCHistogramsNear ( NCBeam beam,
NCType::EEventType  nccc,
NC::SystPars  systPars 
) [private]

Definition at line 260 of file NCExtrapolationFarNear.cxx.

References NC::SystPars::CCBkgScale(), MuELoss::e, FillHistogramSmoothed(), fND_data, fND_MC, 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(), NCType::kCC, kCCBeamNue, kCCBkg, kCCOscNue, kCCSig, kCCTau, NCType::kNC, kNCBeamNue, kNCBkg, kNCOscNue, kNCSig, kNCTau, kNear, Detector::kNear, NC::SystPars::NCBkgScale(), NC::SystPars::NDCleaningScale(), NC::SystPars::ShowerScale(), and NC::SystPars::TrackScale().

Referenced by GetNearMCSpectra(), and WriteResources().

00263 {
00264   using Detector::kNear;
00265 
00266   //loop over NCEnergy bin objects in the near detector
00267   //to get the energies from the data events in each bin
00268   //and fill the fND_data histogram
00269 
00270 
00271   // Alters the systematics 'weights' to reflect the current
00272   // value returned by the fitter.
00273 
00274   int sig         = kNCSig;
00275   int bkg         = kNCBkg;
00276   int tau         = kNCTau;
00277   int beamnue     = kNCBeamNue;
00278   int oscnue      = kNCOscNue;
00279   double bkgShift = systPars.CCBkgScale();
00280   double cleanShift = 1.;
00281 
00282   if(nccc == NCType::kCC){
00283     sig      = kCCSig;
00284     bkg      = kCCBkg;
00285     tau      = kCCTau;
00286     beamnue  = kCCBeamNue;
00287     oscnue   = kCCOscNue;
00288     bkgShift = systPars.NCBkgScale();
00289   }
00290 
00291   fND_data[sig]->Reset( "ICE" );
00292   fND_MC[sig]->Reset( "ICE" );
00293   fND_MC[bkg]->Reset( "ICE" );
00294   fND_MC[tau]->Reset( "ICE" );
00295   fND_MC[beamnue]->Reset( "ICE" );
00296   fND_MC[oscnue]->Reset( "ICE" );
00297 
00298   for(int b  = 0; b < beam->GetNumberEnergyBins(nccc); ++b){
00299     NCEnergyBin* bin = beam->GetEnergyBin(b, nccc);
00300 
00301     const int dataSize       = bin->GetDataVectorSize();
00302     const int mcSize         = bin->GetMCSignalVectorSize();
00303     const int mcSize_bkg     = bin->GetMCBackgroundVectorSize();
00304     const int mcSize_tau     = bin->GetMCNuTauVectorSize();
00305     const int mcSize_beamnue = bin->GetMCBeamNuEVectorSize();
00306     const int mcSize_oscnue  = bin->GetMCOscNuEVectorSize();
00307 
00308     // The variables that get filled with event information
00309     double trueEnergy    = 0;
00310     double recoShowerE   = 0;
00311     double recoMuonE     = 0;
00312     double energy        = 0;
00313     double weight        = 0;
00314     int    flavor        = 0;
00315 
00316     for(int e = 0; e < dataSize; ++e){
00317       bin->GetDataInformation(energy, weight, e);
00318       fND_data[sig]->Fill(energy, weight);
00319     }
00320 
00321     //------------------------------
00322     // numu
00323     for(int e = 0; e < mcSize; ++e){
00324       bin->GetMCInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00325       if(nccc == NCType::kNC) recoMuonE = 0;
00326       cleanShift = systPars.NDCleaningScale(nccc, recoShowerE);
00327       FillHistogramSmoothed(fND_MC[sig], recoShowerE*systPars.ShowerScale(kNear) + recoMuonE*systPars.TrackScale(), weight*cleanShift);
00328     }
00329 
00330     for(int e = 0; e < mcSize_bkg; ++e){
00331       bin->GetMCBackgroundInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00332       if(nccc == NCType::kNC) recoMuonE = 0;
00333       cleanShift = systPars.NDCleaningScale(nccc, recoShowerE);
00334       FillHistogramSmoothed(fND_MC[bkg], recoShowerE*systPars.ShowerScale(kNear) + recoMuonE*systPars.TrackScale(), weight*bkgShift*cleanShift);
00335     }
00336 
00337     //------------------------------
00338     // nutau
00339     for(int e = 0; e < mcSize_tau; ++e){
00340       bin->GetMCNuTauInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00341       if(nccc == NCType::kNC) recoMuonE = 0;
00342       cleanShift = systPars.NDCleaningScale(nccc, recoShowerE);
00343       FillHistogramSmoothed(fND_MC[tau], recoShowerE*systPars.ShowerScale(kNear) + recoMuonE*systPars.TrackScale(), weight*cleanShift);
00344     }
00345 
00346     //------------------------------
00347     // beam nue
00348     for(int e = 0; e < mcSize_beamnue; ++e){
00349       bin->GetMCBeamNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00350       if(nccc == NCType::kNC) recoMuonE = 0;
00351       cleanShift = systPars.NDCleaningScale(nccc, recoShowerE);
00352       FillHistogramSmoothed(fND_MC[beamnue], recoShowerE*systPars.ShowerScale(kNear) + recoMuonE*systPars.TrackScale(), weight*cleanShift);
00353     }
00354 
00355     //------------------------------
00356     // oscillated nue
00357     for(int e = 0; e < mcSize_oscnue; ++e){
00358       bin->GetMCOscNuEInformation(trueEnergy, recoShowerE, recoMuonE, weight, flavor, e);
00359       if(nccc == NCType::kNC) recoMuonE = 0;
00360       cleanShift = systPars.NDCleaningScale(nccc, recoShowerE);
00361       FillHistogramSmoothed(fND_MC[oscnue], recoShowerE*systPars.ShowerScale(kNear) + recoMuonE*systPars.TrackScale(), weight*cleanShift);
00362     }
00363 
00364   } // end for b
00365 }

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

Definition at line 247 of file NCExtrapolationFarNear.cxx.

References HistSmoothHelper().

Referenced by FillDataMCHistogramsNear().

00248 {
00249   HistSmoothHelper(h, false, weight, x);
00250 }

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

Definition at line 253 of file NCExtrapolationFarNear.cxx.

References HistSmoothHelper().

Referenced by FillDataMCHistogramsFar().

00254 {
00255   HistSmoothHelper(h, true, weight, x, y);
00256 }

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

Referenced by WriteResources().

virtual void NCExtrapolationFarNear::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]

This implementation ignores beamsToUse.

Implements NCExtrapolation.

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

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

Implements NCExtrapolation.

Definition at line 45 of file NCExtrapolationFarNear.h.

00045 {return "Far/Near";}

void NCExtrapolationFarNear::GetNearMCSpectra ( NCBeam beam,
NC::SystPars  systs,
TH1 **  nc,
TH1 **  cc 
) [private]

Definition at line 986 of file NCExtrapolationFarNear.cxx.

References FillDataMCHistogramsNear(), fND_MC, NCType::kCC, kCCBeamNue, kCCBkg, kCCOscNue, kCCTau, kEnergyBinsFar, NCType::kNC, kNCBeamNue, kNCBkg, kNCOscNue, and kNCTau.

00988 {
00989   if(nc){
00990     FillDataMCHistogramsNear(beam, NCType::kNC, systs);
00991 
00992     *nc = new TH1D("", "", kNumEnergyBinsFar, kEnergyBinsFar);
00993 
00994     for(int i = 1; i <= kNumEnergyBinsFar; ++i){
00995       const double totalMCNC_nd = fND_MC.at(kNCSig)->GetBinContent(i)
00996         + fND_MC.at(kNCBkg)->GetBinContent(i)      // numu
00997         + fND_MC.at(kNCTau)->GetBinContent(i)      // nutau
00998         + fND_MC.at(kNCBeamNue)->GetBinContent(i)  // beamnue
00999         + fND_MC.at(kNCOscNue)->GetBinContent(i);  // oscnue
01000 
01001       (*nc)->SetBinContent(i, totalMCNC_nd);
01002     } // end for i
01003   }
01004   if(cc){
01005     FillDataMCHistogramsNear(beam, NCType::kCC, systs);
01006 
01007     *cc = new TH1D("", "", kNumEnergyBinsFar, kEnergyBinsFar);
01008 
01009     for(int i = 1; i <= kNumEnergyBinsFar; ++i){
01010       const double totalMCCC_nd = fND_MC.at(kCCSig)->GetBinContent(i)
01011         + fND_MC.at(kCCBkg)->GetBinContent(i)      // numu
01012         + fND_MC.at(kCCTau)->GetBinContent(i)      // nutau
01013         + fND_MC.at(kCCBeamNue)->GetBinContent(i)  // beamnue
01014         + fND_MC.at(kCCOscNue)->GetBinContent(i);  // oscnue
01015 
01016       (*cc)->SetBinContent(i, totalMCCC_nd);
01017     } // end for i
01018   }
01019 }

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

Reimplemented from NCExtrapolation.

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

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

Implements NCExtrapolation.

Definition at line 44 of file NCExtrapolationFarNear.h.

00044 {return "FarNear";}

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

Definition at line 175 of file NCExtrapolationFarNear.cxx.

References fSmoothWidth.

Referenced by FillHistogramSmoothed(), and FillHistogramSmoothed2D().

00176 {
00177   if(fSmoothWidth == 0){
00178     // Shortcut all the below and just fill the histogram
00179     if(is2D)
00180       ((TH2D*)h)->Fill(x, y, weight);
00181     else
00182       h->Fill(x, weight);
00183     return;
00184   }
00185 
00186   // Paranoia over some scary text in the documentation for TAxis::GetBin
00187   h->SetBit(TH1::kCanRebin, false);
00188 
00189   TAxis* xaxis = h->GetXaxis();
00190   const int nbinsx = xaxis->GetNbins();
00191 
00192   // Warning - this number will only work in axis calls, not histogram calls
00193   const int xbin = xaxis->FindBin(x);
00194   const double lo = xaxis->GetBinLowEdge(xbin);
00195   const double hi = xaxis->GetBinUpEdge(xbin);
00196 
00197   const int ybin = is2D ? h->GetYaxis()->FindBin(y) : 0;
00198 
00199   const int histBin = xbin + (nbinsx+2)*ybin;
00200 
00201   const double dist_lo = x-lo;
00202   const double dist_hi = hi-x;
00203 
00204   // Don't smear events into the underflow bin
00205   const bool do_lo = xbin > 1 && dist_lo < fSmoothWidth*x;
00206   // Don't smear events into the overflow bin
00207   const bool do_hi = xbin < nbinsx && dist_hi < fSmoothWidth*x;
00208 
00209   double weight_here = 1;
00210   double weight_lo = 0;
00211   double weight_hi = 0;
00212 
00213   if(do_lo){
00214     const double transfer = .5*(1-dist_lo/(fSmoothWidth*x));
00215     weight_here -= transfer;
00216     weight_lo += transfer;
00217   }
00218   if(do_hi){
00219     const double transfer = .5*(1-dist_hi/(fSmoothWidth*x));
00220     weight_here -= transfer;
00221     weight_hi += transfer;
00222   }
00223 
00224   assert(weight_here >= .499);
00225   assert(weight_lo <= .501);
00226   assert(weight_hi <= .501);
00227   assert(weight_here <= 1);
00228   assert(weight_lo >= 0);
00229   assert(weight_hi >= 0);
00230   assert(fabs(weight_here+weight_lo+weight_hi-1) < .001);
00231 
00232   h->fArray[histBin] += weight*weight_here;
00233 
00234   if(do_lo){
00235     h->fArray[histBin-1] += weight*weight_lo;
00236   }
00237   if(do_hi){
00238     h->fArray[histBin+1] += weight*weight_hi;
00239   }
00240 }

void NCExtrapolationFarNear::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 156 of file NCExtrapolationFarNear.cxx.

References fSmoothWidth, Registry::Get(), Msg::kInfo, and MSG.

00157 {
00158   MSG("NCExtrapolationFarNear", Msg::kInfo) <<  "PREPARING F/N EXTRAPOLATION" << endl;
00159 
00160   // the base class takes care of turning the systematic parameters on/off
00161   // and setting their adjusted values if any
00162 
00163   NCExtrapolation::Prepare(r);
00164 
00165   double      tmpd;
00166 
00167   //  if(r.Get("UE3SqrVal",             tmpd)) fUE3Sqr                     = tmpd;
00168 
00169   if(r.Get("FarNearSmoothingWidth", tmpd)) fSmoothWidth = tmpd;
00170 }

void NCExtrapolationFarNear::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 861 of file NCExtrapolationFarNear.cxx.

References ConstructFarSpectrum(), fFD_data, fFD_fit, fFN, FillDataMCHistogramsFar(), FillDataMCHistogramsNear(), FillResultHistograms(), fNCalls, fND_data, fND_MC, fNom_true_vs_reco, fNominal, fOsc, NCExtrapolation::fUseCC, NCExtrapolation::fUseNC, NCExtrapolation::GetBeam(), NCExtrapolation::GetBestFitOscPars(), NCExtrapolation::GetBestFitSysts(), NCType::kCC, Detector::kFar, Msg::kInfo, BeamType::kL010z185i, NCType::kNC, Detector::kNear, NC::RunUtil::kRunAll, and MSG.

00862 {
00863   MSG("NCExtrapolationFarNear", Msg::kInfo) << "nCalls: " << fNCalls << endl;
00864 
00865   NCBeam* beam = GetBeam(Detector::kFar,
00866                          NCBeam::Info(BeamType::kL010z185i,
00867                                       NC::RunUtil::kRunAll));
00868   NCBeam* beam_nd = GetBeam(Detector::kNear,
00869                             NCBeam::Info(BeamType::kL010z185i,
00870                                          NC::RunUtil::kRunAll));
00871 
00872 
00873   NC::SystPars systPars = GetBestFitSysts();
00874   const NC::OscProb::OscPars* oscPars = GetBestFitOscPars();
00875 
00876   if(fUseNC) FillDataMCHistogramsNear(beam_nd, NCType::kNC, systPars);
00877   if(fUseCC) FillDataMCHistogramsNear(beam_nd, NCType::kCC, systPars);
00878 
00879   fFD_data[kNCSig]->Reset("ICE");
00880   fFD_data[kCCSig]->Reset("ICE");
00881 
00882   if(fUseNC) FillDataMCHistogramsFar(beam, NCType::kNC, systPars);
00883   if(fUseCC) FillDataMCHistogramsFar(beam, NCType::kCC, systPars);
00884 
00885   //resets the histograms - integral contents error
00886   fFD_fit.at(kNCSig)->Reset("ICE");
00887   fFD_fit.at(kCCSig)->Reset("ICE");
00888 
00889   if(oscPars) ConstructFarSpectrum(oscPars);
00890 
00891   delete oscPars;
00892 
00893   // We need the near detector spectrum to be set up so we can extrapolate it
00894   NCExtrapolation::FillResultHistograms(beam_nd);
00895 
00896   // these calls fill the NCBeam histograms using our calculation, which
00897   // allows for Far/Near ratio to be used. NB that 'beam' is only kRunAll
00898   // which means that the RunI plots in NCBeam come out without extrapolation.
00899   if(fUseNC) FillResultHistograms(beam, NCType::kNC, fNominal, fOsc);
00900   if(fUseCC) FillResultHistograms(beam, NCType::kCC, fNominal, fOsc);
00901 
00902   // Now that the beams should be ready, call up to NCExtrapolation
00903 
00904   NCExtrapolation::WriteResources(trueOscPars);
00905 
00906   // get a pointer to the current directory
00907   // this is one of the output files
00908 
00909   TDirectory* saveDir = gDirectory;
00910 
00911   for(unsigned int i = 0; i < fFD_fit.size();  ++i)  fFD_fit[i]->Write();
00912   for(unsigned int i = 0; i < fFN.size();      ++i)      fFN[i]->Write();
00913   for(unsigned int i = 0; i < fND_data.size(); ++i) fND_data[i]->Write();
00914   for(unsigned int i = 0; i < fND_MC.size();   ++i)   fND_MC[i]->Write();
00915   for(unsigned int i = 0; i < fNom_true_vs_reco.size();++i)
00916     fNom_true_vs_reco[i]->Write();
00917 
00918   saveDir->cd();
00919 }


Member Data Documentation

Definition at line 111 of file NCExtrapolationFarNear.h.

Referenced by ConstructFarSpectrum(), and EnableNearToFarExtrapolation().

Definition at line 107 of file NCExtrapolationFarNear.h.

Referenced by FillDataMCHistogramsFar(), and WriteResources().

std::vector<TH1D*> NCExtrapolationFarNear::fFD_fit [private]

Definition at line 103 of file NCExtrapolationFarNear.h.

Referenced by ConstructFarSpectrum(), and WriteResources().

std::vector<TH1D*> NCExtrapolationFarNear::fFN [private]

Far over Near histogram.

Definition at line 101 of file NCExtrapolationFarNear.h.

Referenced by ConstructFarSpectrum(), and WriteResources().

Counts number of calls to the method.

Definition at line 86 of file NCExtrapolationFarNear.h.

Referenced by WriteResources().

std::vector<TH1D*> NCExtrapolationFarNear::fND_data [private]
std::vector<TH1D*> NCExtrapolationFarNear::fND_MC [private]
std::vector<TH2D*> NCExtrapolationFarNear::fNom_true_vs_reco [private]

Histograms for reco fitting work.

Definition at line 98 of file NCExtrapolationFarNear.h.

Referenced by ConstructFarSpectrum(), FillDataMCHistogramsFar(), and WriteResources().

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

Definition at line 88 of file NCExtrapolationFarNear.h.

Referenced by ConstructFarSpectrum(), and WriteResources().

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

Definition at line 89 of file NCExtrapolationFarNear.h.

Referenced by ConstructFarSpectrum(), and WriteResources().

Definition at line 109 of file NCExtrapolationFarNear.h.

Referenced by HistSmoothHelper(), and Prepare().

Definition at line 91 of file NCExtrapolationFarNear.h.

Referenced by ConstructFarSpectrum().

Definition at line 95 of file NCExtrapolationFarNear.h.


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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1