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

NCExtrapolationPID Class Reference

Extrapolation binning events by PID value. More...

#include <NCExtrapolationPID.h>

Inheritance diagram for NCExtrapolationPID:

NCExtrapolation List of all members.

Public Types

typedef std::map< NCBeam::Info,
PIDSpectrum * > 		SpecMap

Public Member Functions

 NCExtrapolationPID ()
virtual ~NCExtrapolationPID ()
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 std::vector< const
TH1 * > 		GetNearMCSpectra (std::vector< NCBeam::Info > beamsToUse)
virtual void CleanupSpectra (std::vector< TH1 * > exps, std::vector< TH1 * > obss)
 Called after FindSpectraForPars() to delete necessary spectra.
virtual void AddEvent (NCEventInfo eventInfo, bool useMCAsData, NCType::EFileType fileType, NCBeam::Info beamInfo)
void DoneFilling ()
virtual 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 void Prepare (const Registry &r)
virtual void Reset (bool nearData, bool farData, bool nearMC, bool farMC)
 Causes the extrapolation to forget everything it knows about the requested spectra.
virtual bool EnableNearToFarExtrapolation (bool enable)
 Enable or disable corrections to FD spectra based on ND spectra.

Static Public Member Functions

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

Protected Member Functions

std::vector< TH1 * > GetExpHists (std::vector< NCBeam::Info > beamsToUse, const NC::OscProb::OscPars *coords, std::vector< TH1 * > ndSpectra=std::vector< TH1 * >()) const
 Get the expected histograms in the format required for NCUtils' log likelihood method.
std::vector< TH1 * > GetObsHists (std::vector< NCBeam::Info > beamsToUse) const
 Get the observed histograms in the format required for NCUtils' log likelihood method.
void WriteSpectra (const NC::OscProb::OscPars *coords, const NC::SystPars &s, std::string suffix="")
 Write the spectra, evaluated at the parameters coords, prettily to gDirectory.
void DoneFillingOneBeam (NCBeam::Info beamInfo)
void WriteMCSpectra (const char *filename)
 Save the unoscillated MC spectra to filename so they can be.
void ReadMCSpectra (const char *filename)
 Read the MC spectra from filename, replacing the existing spectra.

Protected Attributes

SpecMap fSpectra
 The 2D spectra, indexed by beam.
bool fUseCCEnergy
bool fUseSelectionEnergy
int fTrueBinFactor
 How much finer to make the true energy binning than the reco.
int fNumPIDBins
 How many bins in PID to use.
double fPIDmin
 The minimum value of PID to use.
double fPIDmax
 The maximum value of PID to use.
std::string fMCInFilename
 The filename to read the MC spectra from.
std::string fMCOutFilename
 The filename to write the MC spectra to.
bool fDoFarNearCorrection
 Whether to do the far/near correction.
bool fEmulateFarNear
 Whether to emulate the Far/Near method by having just two bins.

Private Member Functions

void PrintMap (POTMap &m, const char *name)
 Print out a POT map for debugging purposes.
void PrintSpectraMap ()
void CopyBinContents (TH1D *hFrom, TH1F *hTo) const
 Copy the bin contents of hFrom into hTo. Used by WriteResources.

Detailed Description

Extrapolation binning events by PID value.

Definition at line 23 of file NCExtrapolationPID.h.

Member Typedef Documentation

typedef std::map<NCBeam::Info, PIDSpectrum*> NCExtrapolationPID::SpecMap
 

Definition at line 27 of file NCExtrapolationPID.h.

Constructor & Destructor Documentation

NCExtrapolationPID::NCExtrapolationPID (   ) 
 

Definition at line 34 of file NCExtrapolationPID.cxx.

00035 :  fUseCCEnergy(true), fUseSelectionEnergy(true), fTrueBinFactor(5), fNumPIDBins(20),
00036   fPIDmin(0.2), fPIDmax(1.1),
00037   fMCInFilename(""), fMCOutFilename(""),
00038   fDoFarNearCorrection(true), fEmulateFarNear(false)
00039 {
00040 }

NCExtrapolationPID::~NCExtrapolationPID (   )  [virtual]
 

Definition at line 44 of file NCExtrapolationPID.cxx.

References fSpectra.

00045 {
00046   for (SpecMap::iterator it=fSpectra.begin(); it!=fSpectra.end(); ++it)
00047     delete it->second;
00048 }

Member Function Documentation

void NCExtrapolationPID::AddEvent (  NCEventInfo  eventInfo,
bool  useMCAsData,
NCType::EFileType  fileType,
NCBeam::Info  beamInfo
)  [virtual]
 

Reimplemented from NCExtrapolation.

Definition at line 53 of file NCExtrapolationPID.cxx.

References NCExtrapolation::AddEvent(), NCEventInfo::analysis, ANtpHeaderInfo::dataType, det, ANtpHeaderInfo::detector, Detector::Detector_t, fEmulateFarNear, fMCInFilename, fNumPIDBins, fPIDmax, fPIDmin, fSpectra, fTrueBinFactor, fUseCCEnergy, NCBeam::Info::GetDescription(), NCEventInfo::header, ANtpRecoInfo::inFiducialVolume, ANtpAnalysisInfo::isCC, ANtpAnalysisInfo::isNC, ANtpRecoInfo::isSimpleCutsClean, ANtpRecoInfo::nuEnergy, ANtpRecoInfo::nuEnergyCC, ANtpRecoInfo::nuEnergyNC, NCEventInfo::reco, ANtpAnalysisInfo::separationParameterCut, NCEventInfo::truth, and ANtpRecoInfo::weight.

00057 {
00058   // This call fills the NCBeam objects with the event. We don't need it
00059   // for the fit, but it does make some pretty plots at the end.
00060   NCExtrapolation::AddEvent(info, useMCAsData, fileType, beamInfo);
00061 
00062   // Some events don't get into the NC or the CC spectrum - exclude
00063   // them if we're emulating the FarNear method
00064   if (fEmulateFarNear &&
00065       !(info.analysis->isNC || info.analysis->isCC)) return;
00066 
00067   // Add the necessary spectrum for this beam if it's not already there
00068   if (fSpectra.find(beamInfo)==fSpectra.end()) {
00069     // If this happened when we got our MC spectra from disk,
00070     // something is probably wrong
00071     assert(fMCInFilename=="");
00072 
00073     string name = string("PID_") + beamInfo.GetDescription().Data();
00074     fSpectra[beamInfo]=new PIDSpectrum(name, "title",
00075                                        fNumPIDBins, fPIDmin, fPIDmax, // PID bins
00076                                        // fNumPIDBins will be -1 if emulating FarNear
00077                                        fTrueBinFactor, beamInfo,
00078                                        info.analysis->separationParameterCut);
00079   }
00080 
00081   // If we loaded the spectra from file, we don't need to add any MC
00082   // events.  But add fake data events in the far detector - this
00083   // isn't the main use case, but it is useful for testing
00084   if (fMCInFilename!="" &&
00085       (info.header->detector==Detector::kNear ||
00086        (info.header->dataType==SimFlag::kMC && !useMCAsData))) return;
00087 
00088   if(!fUseNC && info.analysis->isNC) return;
00089   if(!fUseCC && info.analysis->isCC) return;
00090 
00091   if (!info.reco->inFiducialVolume) return;
00092 
00093   Detector::Detector_t det=Detector::Detector_t(info.header->detector);
00094   if (det==Detector::kNear && !info.reco->isSimpleCutsClean) return;
00095 
00096   double energy;
00097   if (fUseSelectionEnergy) {
00098     energy=info.reco->nuEnergy;
00099   } else {
00100     energy=fUseCCEnergy ? info.reco->nuEnergyCC : info.reco->nuEnergyNC;
00101   }
00102 
00103   // Don't use events with energy<0: this happens for NC energy for
00104   // events with no shower. Also obey energy threshold option.
00105   if (energy<0 || energy<fEnergyThreshold) return;
00106 
00107   if (info.header->dataType==SimFlag::kData || useMCAsData) {
00108     // For fake data, reco->weight contains all the scaling we need (including POT)
00109     fSpectra[beamInfo]->FillData(det,
00110                                 info.truth, info.reco, info.analysis, energy,
00111                                 info.reco->weight);
00112   } else {
00113     fSpectra[beamInfo]->FillMC(det,
00114                               info.truth, info.reco, info.analysis, energy,
00115                               info.reco->weight);
00116   }
00117 
00118 }

void NCExtrapolationPID::CleanupSpectra (  std::vector< TH1 * >  exps,
std::vector< TH1 * >  obss
)  [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.

Definition at line 197 of file NCExtrapolationPID.cxx.

00198 {
00199   // TODO - lots of time goes into deleting histograms.
00200   for (unsigned int n = 0; n < exp.size(); ++n) delete exp[n];
00201 }

void NCExtrapolationPID::CopyBinContents (  TH1D *  hFrom,
TH1F *  hTo
)  const [private]
 

Copy the bin contents of hFrom into hTo. Used by WriteResources.

Definition at line 638 of file NCExtrapolationPID.cxx.

00639 {
00640   for(int i = 0; i < kNumEnergyBinsFar; ++i)
00641     hTo->SetBinContent(i+1, hFrom->GetBinContent(i+1));
00642 }

const Registry & NCExtrapolationPID::DefaultConfig (   )  [static]
 

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

Reimplemented from NCExtrapolation.

Definition at line 478 of file NCExtrapolationPID.cxx.

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

Referenced by Prepare().

00479 {
00480   static Registry r;
00481 
00482   r.UnLockValues();
00483 
00484   // Energy measure to use. Can be "CC", "NC" or "Selection" (which
00485   // uses NC or CC based on the event type given by the selection)
00486   r.Set("PIDFitEnergyType",   "CC");
00487 
00488   // Number of bins in PID to use
00489   r.Set("PIDFitPIDNbins",          20);
00490   // Minimum, maximum value of PID
00491   r.Set("PIDFitPIDmin",          0.2);
00492   r.Set("PIDFitPIDmax",          1.1);
00493 
00494   // How much finer to make the true energy binning than the reco
00495   // energy binning in the PIDSpectrum objects
00496   r.Set("PIDFitTrueBinFactor",   5);
00497 
00498   // The file to read the MC spectra from/write them to
00499   r.Set("PIDFitWriteMCToFile", "");
00500   r.Set("PIDFitReadMCFromFile", "");
00501 
00502   // Whether to do the far/near correction
00503   r.Set("PIDFitDoFarNearCorrection", true);
00504 
00505   // Emulate the FarNear method - overrides the number of bins setting
00506   // and the energy type if set to true
00507   r.Set("PIDFitEmulateFarNear", false);
00508 
00509   r.LockValues();
00510   return r;
00511 
00512 }

void NCExtrapolationPID::DoneFilling (   )  [virtual]
 

Call this when done with filling the MC and data: normalizes the true-to-reco maps, along with anything else that needs doing. This calls DoneFillingOneBeam() for each beam

Reimplemented from NCExtrapolation.

Definition at line 453 of file NCExtrapolationPID.cxx.

References DoneFillingOneBeam(), fMCInFilename, and fSpectra.

00454 {
00455   cout << "NCExtrapolationPID::DoneFilling()" << endl;
00456   // Only need to do this if we didn't load the spectra from file
00457   if(fMCInFilename == ""){
00458     for (SpecMap::iterator sit=fSpectra.begin(); sit!=fSpectra.end(); ++sit) {
00459       DoneFillingOneBeam(sit->first);
00460     }
00461   }
00462 }

void NCExtrapolationPID::DoneFillingOneBeam (  NCBeam::Info  beamInfo  )  [protected]
 

See also:
DoneFilling()

Definition at line 466 of file NCExtrapolationPID.cxx.

References fSpectra, and NCBeam::Info::GetDescription().

Referenced by DoneFilling().

00467 {
00468   cout << "NCExtrapolationPID::DoneFillingOneBeam("
00469        << beamInfo.GetDescription() << ")" << endl;
00470 
00471   cout << "Normalizing true-to-recos" << endl;
00472   fSpectra[beamInfo]->NormalizeTrueToRecos();
00473 
00474 }

bool NCExtrapolationPID::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 588 of file NCExtrapolationPID.cxx.

References fDoFarNearCorrection.

00589 {
00590   const bool ret = fDoFarNearCorrection;
00591   fDoFarNearCorrection = enable;
00592   return ret;
00593 }

void NCExtrapolationPID::FindSpectraForPars (  const NC::OscProb::OscPars *  oscPars,
const NC::SystPars &  systPars,
std::vector< NCBeam::Info >  beamsToUse,
std::vector< TH1 * > &  exps,
std::vector< TH1 * > &  obss
)  [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.

Definition at line 122 of file NCExtrapolationPID.cxx.

References NC::Utility::CloneFast(), GetExpHists(), NC::ISpectrumInterpolator::GetInterpolatedSpectra(), GetNearMCSpectra(), GetObsHists(), NCExtrapolation::InitializeInterpolator(), NC::Utility::MultiplyFast(), s(), and NC::CoordinateConverter::VectorFromSystPars().

00127 {
00128   using namespace NC::Utility;
00129 
00130   if(fDoSystematicInterpolation) InitializeInterpolator(pars);
00131 
00132   if(!fInterpolator){
00133     exps = GetExpHists(beamsToUse, pars);
00134     obss = GetObsHists(beamsToUse);
00135     return;
00136   }
00137 
00138 
00139   const vector<double> shift = fCoordConv.VectorFromSystPars(s);
00140 
00141   vector<TH1*> interp = fInterpolator->GetInterpolatedSpectra(shift);
00142 
00143   vector<TH1*> nd_ratios;
00144   assert(interp.size()%2 == 0); // Every FD has a corresponding ND
00145   // The ND ratios are all the odd numbered ones
00146   for(unsigned int n = 1; n < interp.size(); n += 2){
00147     nd_ratios.push_back(interp[n]);
00148   }
00149 
00150   vector<const TH1*> nd_noms = GetNearMCSpectra(beamsToUse);
00151 
00152   assert(nd_ratios.size() == nd_noms.size());
00153 
00154   vector<TH1*> nd_spectra;
00155   for(unsigned int n = 0; n < nd_noms.size(); ++n){
00156     nd_spectra.push_back(CloneFast(nd_noms[n]));
00157     MultiplyFast(nd_spectra[n], nd_ratios[n]);
00158   }
00159   assert(nd_spectra.size() == nd_noms.size());
00160 
00161   // Pass nd_spectra to GetExpHists so the FD spectra take account of ND shifts
00162   exps = GetExpHists(beamsToUse, pars, nd_spectra);
00163   // Get observed spectra - no interpolation needs to be done on these
00164   obss = GetObsHists(beamsToUse);
00165 
00166   assert(interp.size() == 2*exps.size());
00167   // The FD ratios are the even numbered elements, so their indices correspond
00168   // to twice those in 'exps'
00169   for(unsigned int n = 0; n < exps.size(); ++n){
00170     MultiplyFast(exps[n], interp[2*n]);
00171   }
00172 
00173   // TODO - deleting spectra takes silly amounts of time
00174   for(unsigned int n = 0; n < interp.size(); ++n) delete interp[n];
00175   for(unsigned int n = 0; n < nd_spectra.size(); ++n) delete nd_spectra[n];
00176 }

vector< TH1 * > NCExtrapolationPID::GetExpHists (  std::vector< NCBeam::Info >  beamsToUse,
const NC::OscProb::OscPars *  coords,
std::vector< TH1 * >  ndSpectra = std::vector< TH1 * >()
)  const [protected]
 

Get the expected histograms in the format required for NCUtils' log likelihood method.

Definition at line 328 of file NCExtrapolationPID.cxx.

References fDoFarNearCorrection, and fSpectra.

Referenced by FindSpectraForPars().

00331 {
00332   vector<TH1*> expHists;
00333   for (int i=0; i<(int)beamsToUse.size(); ++i) {
00334     NCBeam::Info beamInfo=beamsToUse[i];
00335     SpecMap::const_iterator spectrumForBeamInfo=fSpectra.find(beamInfo);
00336     assert(spectrumForBeamInfo!=fSpectra.end());
00337     assert(spectrumForBeamInfo->second);
00338     expHists.push_back(spectrumForBeamInfo->second->
00339                        GetPredicted(pars, 
00340                                     "", 
00341                                     fDoFarNearCorrection,
00342                                     ndSpectra.empty() ? 0 : ndSpectra[i]));
00343   }
00344 
00345   return expHists;
00346 }

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

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

Implements NCExtrapolation.

Definition at line 60 of file NCExtrapolationPID.h.

00060 { return "2d PID fit"; }

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

Reimplemented from NCExtrapolation.

Definition at line 181 of file NCExtrapolationPID.cxx.

References fSpectra.

Referenced by FindSpectraForPars().

00182 {
00183   vector<const TH1*> expHists;
00184   for (int i=0; i<(int)beamsToUse.size(); ++i) {
00185     NCBeam::Info beamInfo=beamsToUse[i];
00186     SpecMap::const_iterator spectrumForBeamInfo=fSpectra.find(beamInfo);
00187     assert(spectrumForBeamInfo!=fSpectra.end());
00188     assert(spectrumForBeamInfo->second);
00189     expHists.push_back(spectrumForBeamInfo->second->GetNearMC());
00190   }
00191 
00192   return expHists;
00193 }

vector< TH1 * > NCExtrapolationPID::GetObsHists (  std::vector< NCBeam::Info >  beamsToUse  )  const [protected]
 

Get the observed histograms in the format required for NCUtils' log likelihood method.

Definition at line 350 of file NCExtrapolationPID.cxx.

References fSpectra.

Referenced by FindSpectraForPars().

00352 {
00353   vector<TH1*> obsHists;
00354   for (int i=0; i<(int)beamsToUse.size(); ++i) {
00355     NCBeam::Info beamInfo=beamsToUse[i];
00356     obsHists.push_back(fSpectra.find(beamInfo)->second->GetDataHist());
00357   }
00358 
00359   return obsHists;
00360 }

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

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

Implements NCExtrapolation.

Definition at line 59 of file NCExtrapolationPID.h.

00059 { return "PID"; }

void NCExtrapolationPID::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 516 of file NCExtrapolationPID.cxx.

References DefaultConfig(), fDoFarNearCorrection, fEmulateFarNear, fMCInFilename, fMCOutFilename, fNumPIDBins, fPIDmax, fPIDmin, fSpectra, fTrueBinFactor, fUseCCEnergy, fUseSelectionEnergy, Registry::Get(), Registry::GetInt(), NCExtrapolation::Prepare(), and ReadMCSpectra().

00517 {
00518   NCExtrapolation::Prepare(r);
00519 
00520   cout << "NCExtrapolationPID::Prepare" << endl;
00521   const char* tmps;
00522   int tmpi;
00523   double tmpd;
00524 
00525   if (r.Get("PIDFitEnergyType", tmps)) {
00526     string energyType(tmps);
00527     assert(energyType=="NC" || energyType=="CC" || energyType=="Selection");
00528     fUseSelectionEnergy = (energyType=="Selection");
00529     fUseCCEnergy = (energyType=="CC");
00530   }
00531 
00532   if (r.Get("PIDFitPIDNbins", tmpi)) fNumPIDBins=tmpi;
00533 
00534   // Minimum, maximum value of PID
00535   if (r.Get("PIDFitPIDmin", tmpd)) fPIDmin=tmpd;
00536   if (r.Get("PIDFitPIDmax", tmpd)) fPIDmax=tmpd;
00537 
00538   cout << "fNumPIDBins=" << fNumPIDBins
00539        << " PID range=[" << fPIDmin << ", " << fPIDmax << "]" << endl;
00540 
00541   if (r.Get("PIDFitTrueBinFactor", tmpi))       fTrueBinFactor=tmpi;
00542   if (r.Get("PIDFitDoFarNearCorrection", tmpi)) fDoFarNearCorrection=tmpi;
00543 
00544   if (r.Get("PIDFitEmulateFarNear", tmpi)) fEmulateFarNear=tmpi;
00545 
00546   if (fEmulateFarNear) {
00547     // If we're emulating the FarNear method, we just have two PID
00548     // bins, but we use fNumPIDBins=-1 to signal to PIDSpectrum that it
00549     // should put the separation point where the cut value is
00550     fNumPIDBins=-1;
00551     fUseSelectionEnergy=true;
00552 
00553     if (r.Get("PIDFitPIDNbins", tmpi) &&
00554         tmpi!=DefaultConfig().GetInt("PIDFitPIDNbins"))
00555       cout << "WARNING: PIDFitPIDNbins was set, but overridden by PIDFitEmulateFarNear"
00556            << endl;
00557   }
00558 
00559   // Files to read the MC spectra from/write them to
00560   if (r.Get("PIDFitReadMCFromFile", tmps)) fMCInFilename=tmps;
00561   if (r.Get("PIDFitWriteMCToFile", tmps))  fMCOutFilename=tmps;
00562 
00563   if (fMCInFilename!="") {
00564     ReadMCSpectra(fMCInFilename.c_str());
00565     // reset the observed spectra - if we got the PIDSpectrum's from file, we'll be
00566     // filling them with mock data presumably anyway
00567     for (SpecMap::iterator sit=fSpectra.begin(); sit!=fSpectra.end(); ++sit) {
00568       TH2F* h=sit->second->GetDataHist();
00569       h->Reset();
00570     }
00571   }
00572 
00573 }

void NCExtrapolationPID::PrintMap (  POTMap &  m,
const char *  name
)  [private]
 

Print out a POT map for debugging purposes.

Definition at line 304 of file NCExtrapolationPID.cxx.

00305 {
00306   cout << "Map " << name << endl;
00307   for (POTMap::const_iterator it=m.begin();
00308        it!=m.end(); ++it) {
00309     cout << it->first.GetDescription() << "\t" << it->second << endl;
00310   }
00311 }

void NCExtrapolationPID::PrintSpectraMap (   )  [private]
 

Definition at line 315 of file NCExtrapolationPID.cxx.

References fSpectra.

00316 {
00317   cout << "fSpectra:" << endl;
00318   for (SpecMap::const_iterator it=fSpectra.begin();
00319        it!=fSpectra.end(); ++it) {
00320     cout << it->first.GetDescription() << "\t" << hex << it->second << dec
00321          << "\t" << it->second->GetBeamInfo().GetTitle() << endl;
00322   }
00323 }

void NCExtrapolationPID::ReadMCSpectra (  const char *  filename  )  [protected]
 

Read the MC spectra from filename, replacing the existing spectra.

Definition at line 607 of file NCExtrapolationPID.cxx.

References fSpectra, Registry::Get(), PIDSpectrum::GetBeamInfo(), and PIDSpectrum::GetNPIDBins().

Referenced by Prepare().

00608 {
00609 
00610   cout << "Reading MC spectra from file " << filename << endl;
00611   TFile f(filename);
00612 
00613   TList *l = f.GetListOfKeys();
00614   TIter it(l);
00615   TObject *obj;
00616   TKey* k;
00617   while (( k=(TKey*)it() )) {
00618     obj=f.Get(k->GetName());
00619     cout << "Trying object " << obj->GetName() << endl;
00620     if (! (obj->IsA() == PIDSpectrum::Class()) ) {
00621       cerr << "Object with name " << obj->GetName() << " not a PIDSpectrum. Bailing" << endl;
00622       exit(1);
00623     }
00624     cout << "...object is a PIDSpectrum. Adding to list" << endl;
00625     PIDSpectrum* p=(PIDSpectrum*)obj;
00626 
00627     // fNumPIDBins is mangled when using fEmulateFarNear, so need this branch
00628     if (fEmulateFarNear) assert(p->GetNPIDBins()==2);
00629     else                 assert(p->GetNPIDBins()==fNumPIDBins);
00630 
00631     fSpectra[p->GetBeamInfo()]=p;
00632   }
00633 
00634 }

void NCExtrapolationPID::Reset (  bool  nearData,
bool  farData,
bool  nearMC,
bool  farMC
)  [virtual]
 

Causes the extrapolation to forget everything it knows about the requested spectra.

This implementation clears the relevant NCBeam objects

Derived classes please override if you handle spectra yourself

Reimplemented from NCExtrapolation.

Definition at line 576 of file NCExtrapolationPID.cxx.

References fSpectra.

00578 {
00579   for (SpecMap::iterator sit=fSpectra.begin(); sit!=fSpectra.end(); ++sit) {
00580     if(nearData) sit->second->ResetNearData();
00581     if(farData) sit->second->ResetData();
00582     if(nearMC) sit->second->ResetNearMC();
00583     if(farMC) sit->second->ResetMC();
00584   }
00585 }

void NCExtrapolationPID::WriteMCSpectra (  const char *  filename  )  [protected]
 

Save the unoscillated MC spectra to filename so they can be.

Definition at line 597 of file NCExtrapolationPID.cxx.

References fSpectra.

Referenced by WriteResources().

00598 {
00599   TFile f(filename, "RECREATE");
00600   for (SpecMap::iterator sit=fSpectra.begin(); sit!=fSpectra.end(); ++sit) {
00601     sit->second->Write();
00602   }
00603 }

void NCExtrapolationPID::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 206 of file NCExtrapolationPID.cxx.

References COPYBINS, NCType::EEventType, fMCOutFilename, fSpectra, NCExtrapolation::GetBestFitOscPars(), NCExtrapolation::GetBestFitSysts(), PIDSpectrum::GetDataHist(), PIDSpectrum::GetPseudoNCCCComponent(), PIDSpectrum::GetPseudoNCCCSpectrum(), NCBeam::MarkHistogramsFilled(), WriteMCSpectra(), NCExtrapolation::WriteResources(), and WriteSpectra().

00207 {
00208 
00209   NC::OscProb::OscPars* bestFit=GetBestFitOscPars();
00210   NC::OscProb::OscPars* oscParsNoOsc = new NC::OscProb::NoOscillations;
00211 
00212   // Hack to fill NCBeam's histograms with something semi-reasonable
00213   // that includes the F/N ratio, for ease of comparison with other
00214   // methods
00215   
00216   // Loop over all of the beams
00217   for(fBeams_t::iterator it=fBeams.begin(); it!=fBeams.end(); ++it){
00218     NCBeam::Info beamInfo=it->first.first;
00219     // it->first.second is the detector. We'll just do the FD for now, so ignore it
00220     NCBeam* beam=it->second;
00221     assert(fSpectra.find(beamInfo)!=fSpectra.end());
00222     PIDSpectrum* spectrum=fSpectra[beamInfo];
00223     
00224     for(int inccc=NCType::kNC; inccc<=NCType::kCC; ++inccc){
00225       // Gah, I still don't understand int <-> enum implicit casting rules
00226       NCType::EEventType nccc=NCType::EEventType(inccc);
00227       TH1D* pseudoNoOsc=spectrum->GetPseudoNCCCSpectrum(nccc, PIDSpectrum::kAll, 
00228                                                         oscParsNoOsc);
00229       TH1D* pseudoFit=spectrum->GetPseudoNCCCSpectrum(nccc, PIDSpectrum::kAll, 
00230                                                       bestFit);
00231       TH1D* pseudoCC=spectrum->GetPseudoNCCCSpectrum(nccc, PIDSpectrum::kNumu, 
00232                                                      bestFit);
00233       TH1D* pseudoBeamNue=spectrum->GetPseudoNCCCSpectrum(nccc, PIDSpectrum::kBeamNue, 
00234                                                           bestFit);
00235       TH1D* pseudoOscNue=spectrum->GetPseudoNCCCSpectrum(nccc, PIDSpectrum::kOscNue, 
00236                                                          bestFit);
00237       TH1D* pseudoNutau=spectrum->GetPseudoNCCCSpectrum(nccc, PIDSpectrum::kNutau, 
00238                                                         bestFit);
00239       
00240       TH1D* pseudoData=spectrum->GetPseudoNCCCComponent(spectrum->GetDataHist(), nccc);
00241       
00242       // The NCBeam plots want background for both spectra (ie, CC in
00243       // the NC spectrum, and NC in the CC spectrum), but we don't
00244       // have the NC component alone, so we have to form it here
00245       TH1D* pseudoBG;
00246       if(nccc==NCType::kNC) pseudoBG=pseudoCC;
00247       else{
00248         // Form the NC BG in the CC spectrum as the difference between
00249         // the total and the various CC parts
00250         pseudoBG=(TH1D*)pseudoFit->Clone("pseudoBG");
00251         pseudoBG->SetDirectory(0);
00252         
00253         pseudoBG->Add(pseudoCC, -1);
00254         // These are going to be tiny and not matter, but may as well
00255         // subtract them anyway
00256         pseudoBG->Add(pseudoBeamNue, -1);
00257         pseudoBG->Add(pseudoOscNue, -1);
00258         pseudoBG->Add(pseudoNutau, -1);
00259       }
00260       
00261 #define COPYBINS(from, to) CopyBinContents(from,                        \
00262                                            beam->Get##to##Histogram(nccc)) 
00263       
00264       COPYBINS(pseudoData,    Data);
00265       COPYBINS(pseudoNoOsc,   MC);
00266       COPYBINS(pseudoFit,     MCFit);
00267       COPYBINS(pseudoBG,      MCBackground);
00268       COPYBINS(pseudoBeamNue, MCFitNuEToNuE);
00269       COPYBINS(pseudoOscNue,  MCFitNuMuToNuE);
00270       COPYBINS(pseudoNutau,   MCFitNuMuToNuTau);
00271 
00272 #undef COPYBINS
00273 
00274     }      
00275     // We filled the histograms manually, and don't need NCBeam's automatic stuff
00276     // to happen. So inform it of that fact.
00277     beam->MarkHistogramsFilled();
00278   }
00279 
00280   delete oscParsNoOsc;
00281 
00282   NCExtrapolation::WriteResources(trueOscPars);
00283 
00284   // Done with the NCBeam stuff, so now we can get on with making my plots
00285     
00286   // Store my plots in a debug directory
00287   gDirectory->mkdir("debug", "debug")->cd();
00288     
00289   if(bestFit){
00290     NC::SystPars systs=GetBestFitSysts();
00291     
00292     WriteSpectra(bestFit, systs, "bestFit");
00293     
00294     if(trueOscPars) WriteSpectra(trueOscPars, systs, "truePars");
00295     
00296   }
00297     
00298   // If requested, save the PIDSpectrum objects out to file
00299   if (fMCOutFilename!="") WriteMCSpectra(fMCOutFilename.c_str());
00300 }

void NCExtrapolationPID::WriteSpectra (  const NC::OscProb::OscPars *  coords,
const NC::SystPars &  s,
std::string  suffix = ""
)  [protected]
 

Write the spectra, evaluated at the parameters coords, prettily to gDirectory.

Definition at line 364 of file NCExtrapolationPID.cxx.

References PIDSpectrum::DrawEnergySlices(), PIDSpectrum::DrawNearSpectra(), PIDSpectrum::DrawPIDSlices(), PIDSpectrum::DrawSpectrum(), PIDSpectrum::EventTypeList(), fSpectra, PIDSpectrum::GetBeamInfo(), NCBeam::Info::GetDescription(), NC::ISpectrumInterpolator::GetInterpolatedSpectra(), PIDSpectrum::GetOnePredicted(), NCBeam::Info::GetRunType(), NCBeam::Info::IsNominal(), MSG, s(), and NC::CoordinateConverter::VectorFromSystPars().

Referenced by WriteResources().

00367 {
00368   if(suffix!="") gDirectory->mkdir(suffix.c_str(), suffix.c_str())->cd();
00369 
00370   vector<TH1*> fd_ratios;
00371   vector<TH1*> nd_ratios;
00372 
00373   if (fInterpolator) {
00374     const vector<double> shift = fCoordConv.VectorFromSystPars(s);
00375     vector<TH1*> ratios = fInterpolator->GetInterpolatedSpectra(shift);
00376 
00377     assert(ratios.size()%2 == 0); // near and far
00378 
00379     for(unsigned int n = 0; n < ratios.size(); n += 2)
00380       nd_ratios.push_back(ratios[n]);
00381     for(unsigned int n = 1; n < ratios.size(); n += 2)
00382       fd_ratios.push_back(ratios[n]);
00383   }
00384 
00385   int specIdx = 0;
00386 
00387   for (SpecMap::iterator sit=fSpectra.begin(); sit!=fSpectra.end(); ++sit) {
00388     PIDSpectrum* spec=sit->second;
00389 
00390     TString beamName=spec->GetBeamInfo().GetDescription();
00391     assert(beamName!="");
00392     gDirectory->mkdir(beamName, beamName)->cd();
00393 
00394     TH1* nd_ratio = 0;
00395     TH1* fd_ratio = 0;
00396 
00397     // Need an interpolator obviously. Interpolating the shifted spectra makes
00398     // no sense. Interpolator doesn't return a kRunAll ratio which
00399     // unfortunately means that the output RunAll spectrum won't be
00400     // interpolated at all.
00401     const bool canInterpolate =
00402       fInterpolator &&
00403       spec->GetBeamInfo().IsNominal() &&
00404       spec->GetBeamInfo().GetRunType() != NC::RunUtil::kRunAll;
00405 
00406     if(canInterpolate){
00407       nd_ratio=nd_ratios[specIdx];
00408       fd_ratio=fd_ratios[specIdx];
00409 
00410       // These should match up, if not we've done the indexing wrong
00411       MSG("NCExtrapolationPID", Msg::kDebug)
00412         << "Wanted " << beamName << " found "
00413         << nd_ratio->GetName() << " and "
00414         << fd_ratio->GetName() << endl;
00415 
00416       ++specIdx;
00417     }
00418 
00419     vector<TCanvas*> v;
00420     v.push_back( spec->DrawSpectrum(coords, fd_ratio, suffix) );
00421 
00422     vector<TCanvas*> vNear=spec->DrawNearSpectra(nd_ratio, suffix);
00423     vector<TCanvas*> vE=spec->DrawEnergySlices(coords, fd_ratio, suffix);
00424     vector<TCanvas*> vPID=spec->DrawPIDSlices(coords, fd_ratio, suffix);
00425 
00426     v.insert(v.end(), vE.begin(),    vE.end());
00427     v.insert(v.end(), vPID.begin(),  vPID.end());
00428     v.insert(v.end(), vNear.begin(), vNear.end());
00429 
00430     for (vector<TCanvas*>::iterator it=v.begin();
00431          it!=v.end(); ++it) {
00432       (*it)->Write();
00433     }
00434 
00435     vector<PIDSpectrum::evtType> evtTypes=spec->EventTypeList();
00436     for (int j=0; j<(int)evtTypes.size(); ++j) {
00437       spec->GetOnePredicted(evtTypes[j], coords, suffix.c_str())->Write();
00438     }
00439 
00440     gDirectory->cd("..");
00441   }
00442 
00443 
00444   gDirectory->cd("..");
00445   // Also get the histograms of the tau component
00446   // (for the "can we see taus?" study)
00447   //for (SpecMap::iterator sit=fSpectra.begin(); sit!=fSpectra.end(); ++sit)
00448   //  sit->second->GetOnePredicted(PIDSpectrum::kNutau, coords, suffix.c_str())->Write();
00449 }

Member Data Documentation

bool NCExtrapolationPID::fDoFarNearCorrection [protected]
 

Whether to do the far/near correction.

Definition at line 137 of file NCExtrapolationPID.h.

Referenced by EnableNearToFarExtrapolation(), GetExpHists(), and Prepare().

bool NCExtrapolationPID::fEmulateFarNear [protected]
 

Whether to emulate the Far/Near method by having just two bins.

Definition at line 143 of file NCExtrapolationPID.h.

Referenced by AddEvent(), and Prepare().

std::string NCExtrapolationPID::fMCInFilename [protected]
 

The filename to read the MC spectra from.

Definition at line 131 of file NCExtrapolationPID.h.

Referenced by AddEvent(), DoneFilling(), and Prepare().

std::string NCExtrapolationPID::fMCOutFilename [protected]
 

The filename to write the MC spectra to.

Definition at line 134 of file NCExtrapolationPID.h.

Referenced by Prepare(), and WriteResources().

int NCExtrapolationPID::fNumPIDBins [protected]
 

How many bins in PID to use.

Definition at line 122 of file NCExtrapolationPID.h.

Referenced by AddEvent(), and Prepare().

double NCExtrapolationPID::fPIDmax [protected]
 

The maximum value of PID to use.

Definition at line 128 of file NCExtrapolationPID.h.

Referenced by AddEvent(), and Prepare().

double NCExtrapolationPID::fPIDmin [protected]
 

The minimum value of PID to use.

Definition at line 125 of file NCExtrapolationPID.h.

Referenced by AddEvent(), and Prepare().

SpecMap NCExtrapolationPID::fSpectra [protected]
 

The 2D spectra, indexed by beam.

Definition at line 102 of file NCExtrapolationPID.h.

Referenced by AddEvent(), DoneFilling(), DoneFillingOneBeam(), GetExpHists(), GetNearMCSpectra(), GetObsHists(), Prepare(), PrintSpectraMap(), ReadMCSpectra(), Reset(), WriteMCSpectra(), WriteResources(), WriteSpectra(), and ~NCExtrapolationPID().

int NCExtrapolationPID::fTrueBinFactor [protected]
 

How much finer to make the true energy binning than the reco.

Definition at line 119 of file NCExtrapolationPID.h.

Referenced by AddEvent(), and Prepare().

bool NCExtrapolationPID::fUseCCEnergy [protected]
 

Whether to use the CC energy estimator. If false, use NC energy estimator. useSelectionEnergy overrides this

Definition at line 108 of file NCExtrapolationPID.h.

Referenced by AddEvent(), and Prepare().

bool NCExtrapolationPID::fUseSelectionEnergy [protected]
 

Whether to use the event energy as defined by NCUtils - ie CC energy if event is selected CC and NC energy if event is selected NC. Overrides fUseCCEnergy

Definition at line 115 of file NCExtrapolationPID.h.

Referenced by Prepare().

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