PIDSpectrum Class Reference

Provides a 2d energy vs PID spectrum for use by NCExtrapolationPID . Holds separate 2d histograms for the NC, CC, beam nue, osc nue and nutau components, along with one for the data. More...

#include <PIDSpectrum.h>

List of all members.

Public Types
enum	evtType {   kNC = 0, kNumu = 1, kBeamNue = 2, kOscNue = 3,   kNutau = 4, kAll = 5 }
Public Member Functions
	PIDSpectrum ()
	Default constructor: Needed for persisting to disk.
	PIDSpectrum (std::string _name, std::string _title, int _nPIDBins, double _PIDmin, double _PIDmax, int _trueBinFactor, NCBeam::Info _beamInfo, double _cutValue=-9999)
virtual	~PIDSpectrum ()
std::vector< PIDSpectrum::evtType >	EventTypeList () const
	The list of event types available.
TH2F *	GetPredicted (const NC::OscProb::OscPars *coords, const char *histname, bool doFN, TH1 *ndSpectrum=0) const
	Get the predicted total histogram for the given coords.
const TH2F *	GetNearMC () const
	Get the near detector Monte Carlo spectrum.
TH2F *	GetOnePredicted (evtType t, const NC::OscProb::OscPars *coords, const char *histname=0) const
	Get the predicted histogram for one event type for the given coords.
void	AddOnePredicted (TH2F *, evtType t, const NC::OscProb::OscPars *coords) const
	Add the predicted histogram for one event type for the given coords into another.
TH1D *	GetPseudoNCCCSpectrum (NCType::EEventType nccc, evtType t, const NC::OscProb::OscPars *coords) const
	Return the "NC" or "CC" spectrum for true event type t at coords.
TH1D *	GetPseudoNCCCComponent (TH2F *h, NCType::EEventType nccc) const
	Return the 1D part of h that is type nccc in the manner of GetPseudoNCCCSpectrum.
TH2F *	GetDataHist () const
	Get the data histogram for this Spectrum.
virtual void	FillMC (Detector::Detector_t det, const ANtpTruthInfoBeam *t, const ANtpRecoInfo *r, const ANtpAnalysisInfo *ana, double E, double scale)
	Add an event to the appropriate MC histogram.
virtual void	FillData (Detector::Detector_t det, const ANtpTruthInfoBeam *t, const ANtpRecoInfo *r, const ANtpAnalysisInfo *ana, double E, double scale=1)
	Add an event to the data histogram.
virtual void	NormalizeTrueToRecos ()
std::map< PIDSpectrum::evtType, TH3F * >	GetTrueToRecoMap () const
	Get the true-to-reco maps for each event type.
std::map< PIDSpectrum::evtType, TH2F * >	GetTrueEMap () const
	Get the true energy spectra for each event type.
NCBeam::Info	GetBeamInfo ()
	Get the beam index for this spectrum.
TCanvas *	DrawSpectrum (const NC::OscProb::OscPars *coords, TH1 *systratios, std::string suffix="") const
std::vector< TCanvas * >	DrawEnergySlices (const NC::OscProb::OscPars *coords, TH1 *systratios, std::string suffix) const
	Draw slices of the energy spectrum for a given set of coords.
std::vector< TCanvas * >	DrawPIDSlices (const NC::OscProb::OscPars *coords, TH1 *systratios, std::string suffix) const
	Draw slices of the PID spectrum for a given set of coords.
std::vector< TCanvas * >	DrawNearSpectra (TH1 *systratios, std::string suffix) const
	Draw the near spectra onto a canvas.
int	GetNPIDBins () const
	Return the number of bins in PID.
void	ResetData ()
	Clears the data histogram.
void	ResetMC ()
	Clears all the Monte Carlo histograms.
void	ResetNearData ()
	Clears the ND data histogram.
void	ResetNearMC ()
	Clears the ND Monte Carlo histogram.
Protected Member Functions
std::string	EvtTypeAsString (evtType t) const
	Return event type as string.
evtType	GetEvtType (const ANtpTruthInfoBeam *t) const
	Get the event type for a given truth.
double	GetEventWeight (evtType t, double trueE, const NC::OscProb::OscPars *coords) const
	Return the oscillation weight for event of type t with true energy trueE at oscillation parameters coords.
virtual TH1F *	GetOscWeightHist (evtType t, const NC::OscProb::OscPars *coords, int nBins, Double_t *binEdges) const
	Get the histogram of oscillation weights in true energy for the given event type.
double	PIDWithinRange (double origPID) const
	Bring origPID within the range [PIDmin, PIDmax ] so that nothing goes into the overflow bins.
Protected Attributes
std::map< evtType, TH2F * >	fCmpTrue
	The true energy components.
std::map< evtType, TH3F * >	fTrueToReco
	The true-to-reco matrices.
TH2F *	fData
	The data spectrum.
TH2F *	fNearMC
	The near MC.
TH2F *	fNearData
	The near data.
std::vector< evtType >	fEvtTypes
	The available event types so we can loop over them.
int	fNumPIDBins
double	fPIDmin
double	fPIDmax
int	fNumBinEdgesDontUseThis
	The number of bin edges.
Double_t *	fPIDBins
int	fTrueBinFactor
	How much finer to make the true energy binning than the reco energy binning.
NCBeam::Info	fBeamInfo
	The beam info for this spectrum.
std::vector< int >	fDrawColors
	Colours for the components of the selected spectrum.
double	fCutValue
	The NC/CC cut value. Only used for the NCBeam plot filling,.
Private Member Functions
Double_t *	MultiplyBins (int factor, int nBinsOrig, const Double_t *origBins) const
	Make a binning based on origBins, but factor times finer.
TH1D *	ProjectHist (TH2 *h, std::string axis, std::string histname, int lobin=0, int hibin=-1) const
	Draw the energy or pid projection of histogram h.
std::vector< TCanvas * >	DrawSlices (std::string variable, const NC::OscProb::OscPars *coords, TH1 *systratios, std::string suffix) const
	Draw the energy or pid slices onto a canvas and return it.
void	FarNearCorrect (TH2F *fdPred, const TH2F *ndData, const TH1 *ndMC) const
	Do the Far/Near correction, ie, multiply fdPred by ndData over ndMC.
	ClassDef (PIDSpectrum, 9)

---

Detailed Description

Provides a 2d energy vs PID spectrum for use by NCExtrapolationPID . Holds separate 2d histograms for the NC, CC, beam nue, osc nue and nutau components, along with one for the data.

The class also performs a simple far/near extrapolation using near detector data and MC

Definition at line 40 of file PIDSpectrum.h.

---

Member Enumeration Documentation

enum PIDSpectrum::evtType

Enumeration values: kNC  kNumu  kBeamNue  kOscNue  kNutau  kAll  Definition at line 62 of file PIDSpectrum.h. Referenced by FillMC(), and GetEvtType(). { kNC=0, kNumu=1, kBeamNue=2, kOscNue=3, kNutau=4, kAll=5 };

---

Constructor & Destructor Documentation

PIDSpectrum::PIDSpectrum (   )  [inline]

Default constructor: Needed for persisting to disk. Definition at line 44 of file PIDSpectrum.h. References fData, fNearData, fNearMC, and fPIDBins. : fData(0), fNearMC(0), fNearData(0), fPIDBins(0) { };

PIDSpectrum::PIDSpectrum (  std::string  _name, std::string  _title, int  _nPIDBins, double  _PIDmin, double  _PIDmax, int  _trueBinFactor, NCBeam::Info  _beamInfo, double  _cutValue = -9999 )

Normal constructor. If _nPIDBins is -1, emulate the FarNear method by having two bins separated at the cut parameter. The cut value is given by the _cutValue argument, and is only used when emulating the FarNear method Definition at line 21 of file PIDSpectrum.cxx. References Nav::GetName(), and kNumEnergyBinsFar. : TNamed(_name, _title), fNumPIDBins(_nPIDBins), fPIDmin(_PIDmin), fPIDmax(_PIDmax), fTrueBinFactor(_trueBinFactor), fBeamInfo(_beamInfo), fCutValue(_cutValue) { fEvtTypes.push_back(kNC); fEvtTypes.push_back(kNumu); fEvtTypes.push_back(kBeamNue); fEvtTypes.push_back(kOscNue); fEvtTypes.push_back(kNutau); fDrawColors.resize(fEvtTypes.size()); fDrawColors[kNC]=kRed; fDrawColors[kNumu]=kBlue; fDrawColors[kBeamNue]=kGreen+2; fDrawColors[kOscNue]=kOrange-3; fDrawColors[kNutau]=kMagenta+2; // Set up the binning for true energy axes Double_t* trueEBins=MultiplyBins(fTrueBinFactor, kNumEnergyBinsFar, kEnergyBinsFar); int nTrueEBins=kNumEnergyBinsFar*fTrueBinFactor; // Gah, have to construct this manually, because ROOT won't let me // mix-and-match the auto and manual binning for TH3s if (fNumPIDBins==-1) { // We're emulating the FarNear method - two bins, separated at the cut value fPIDBins=new Double_t[3]; fPIDBins[0]=fPIDmin; fPIDBins[2]=fPIDmax; fPIDBins[1]=fCutValue; // HACK: we probably shouldn't change this from what was passed // in, but it makes the rest of the code easier fNumPIDBins=2; } else { fPIDBins=new Double_t[fNumPIDBins+1]; for (int i=0; i<fNumPIDBins+1; ++i) fPIDBins[i]=fPIDmin+float(i)*(fPIDmax-fPIDmin)/fNumPIDBins; } // This is the one place where we *should* use this variable fNumBinEdgesDontUseThis=fNumPIDBins+1; for (vector<evtType>::iterator it=fEvtTypes.begin(); it!=fEvtTypes.end(); ++it) { string tmp("hTrue"); tmp+=GetName(); tmp+=EvtTypeAsString(*it); //cout << "Creating true energy hist" << endl; fCmpTrue[*it]=new TH2F(tmp.c_str(), "True energy spectrum;True Energy (GeV);PID", nTrueEBins, trueEBins, fNumPIDBins, fPIDBins); tmp="t2r_"; tmp+=GetName(); tmp+=EvtTypeAsString(*it); //cout << "Creating true to reco hist" << endl; fTrueToReco[*it]=new TH3F(tmp.c_str(), "True-to-reco;Reco Energy;PID;True energy", kNumEnergyBinsFar, kEnergyBinsFar, fNumPIDBins, fPIDBins, nTrueEBins, trueEBins); } delete[] trueEBins; string tmp("hData"); tmp+=GetName(); //cout << "Creating data hist" << endl; fData=new TH2F(tmp.c_str(), "Far Data;Reco Energy;PID", kNumEnergyBinsFar, kEnergyBinsFar, fNumPIDBins, fPIDBins); tmp="hNearMC"; tmp+=GetName(); fNearMC=new TH2F(tmp.c_str(), "Near MC;Reco Energy;PID", kNumEnergyBinsFar, kEnergyBinsFar, fNumPIDBins, fPIDBins); tmp="hNearData"; tmp+=GetName(); fNearData=new TH2F(tmp.c_str(), "Near Data;Reco Energy;PID", kNumEnergyBinsFar, kEnergyBinsFar, fNumPIDBins, fPIDBins); }

virtual PIDSpectrum::~PIDSpectrum (   )  [inline, virtual]

Definition at line 58 of file PIDSpectrum.h. { if (fPIDBins) delete[] fPIDBins; }

---

Member Function Documentation

void PIDSpectrum::AddOnePredicted (  TH2F *  , evtType  t, const NC::OscProb::OscPars *  coords )  const

Add the predicted histogram for one event type for the given coords into another. Parameters: onePred  Histogram to add to. Must be kNumEnergyBinsFar by fNumPIDBins t  The event type to evaluate for coords  The oscillation/systematic parameters to evaluate at Definition at line 180 of file PIDSpectrum.cxx. References fCmpTrue, fNumPIDBins, fTrueBinFactor, fTrueToReco, GetOscWeightHist(), kNumEnergyBinsFar, and MultiplyBins(). Referenced by GetOnePredicted(), and GetPredicted(). { // Set up the binning for true energy axes Double_t* trueEBins=MultiplyBins(fTrueBinFactor, kNumEnergyBinsFar, kEnergyBinsFar); const int nTrueEBins=kNumEnergyBinsFar*fTrueBinFactor; TH1F* oscProb=GetOscWeightHist(t, coords, nTrueEBins, trueEBins); delete[] trueEBins; /* // This is some debugging code that might turn out to be useful again one day if (oscProb->Integral()==0 && (t==kNumu || t==kNC)) { cout << "oscprob empty for type " << EvtTypeAsString(t) << endl; cout << "Coords: "; for (int i=0; i<(int)coords.size(); ++i) cout << coords[i] << " "; cout << endl; } static bool doneAlready=false; if (t==kNC && !doneAlready) { oscProb->Write(); doneAlready=true; } */ // Convolve it all together // Need to use find() because map::operator[] isn't const TH3F* t2r=fTrueToReco.find(t)->second; TH2F* trueCmp=fCmpTrue.find(t)->second; /* if ((t==kNC || t==kNumu) && hit->second->Integral()==0) { cout << "t2r empty for type " << EvtTypeAsString(t) << endl; cout << "Coords: "; for (int i=0; i<(int)coords.size(); ++i) cout << coords[i] << " "; cout << endl; } if ((t==kNC || t==kNumu) && hit2->second->Integral()==0) { cout << "fCmpTrue empty for type " << EvtTypeAsString(t) << endl; cout << "Coords: "; for (int i=0; i<(int)coords.size(); ++i) cout << coords[i] << " "; cout << endl; } */ // TH2 and TH3 convert to 1D bin numbers when asked for a bin by 2D // number. It turns out that this is rather slow (probably since // it's all implemented by the TH1 base class). Since I know I won't // be accessing past the end of the range, and I know what dimension // the histogram has, I can calculate this manually for (hopefully) // a speedup // GetBinContent does bounds-checking, which we don't need, // since the loop limits ensure that we only access real // bins. Instead, use the internal array fArray, which for // some reason is public, to access the bin values without // incurring the time penalty from bounds-checking // // Similarly SetBinContent can be replaced by a direct write to fArray // The disadvantage here is that bin errors etc won't be set, but we // don't need them here // The arrangement of t2r in memory in fArray is my trueE, PID, recoE // most significant first. By putting the loops in this same order // we consider every point in the order they are layed out in memory. // // onePred has the same arrangement, except it doesn't have a trueE direction // // It is important we loop through every single bin (including underflow // and overflow) or the counting of the index will go wrong. // Begin at the start of the true-to-reco array float* pt2r = t2r->fArray; for(int iTrueE = 0; iTrueE <= nTrueEBins+1; ++iTrueE){ const float oscWeight = oscProb->fArray[iTrueE]; // For each new true energy go back to the start of the output array float* pOnePred = onePred->fArray; for(int iPID=0; iPID <= fNumPIDBins+1; ++iPID){ // The trueCmp histogram has its axes the wrong way round, so have // to do this calculation instead of being able to use the same trick // as for t2r and onePred const int trueE_PIDBin = iTrueE+(nTrueEBins+2)*iPID; const float pidWeight = trueCmp->fArray[trueE_PIDBin]; const float oscPidWeight = oscWeight * pidWeight; for(int iRecoE = 0; iRecoE <= kNumEnergyBinsFar+1; ++iRecoE){ const float fillValue = (*pt2r)*oscPidWeight; *pOnePred += fillValue; // Advance to the next entry in each array ++pt2r; ++pOnePred; } // end for iRecoE } // end for iPID } // end for iTrueE /* if ((t==kNC || t==kNumu) && onePred->Integral()<1e-6) { cout << "predicted spectrum empty for type " << EvtTypeAsString(t) << endl; cout << "Coords: "; for (int i=0; i<(int)coords.size(); ++i) cout << coords[i] << " "; cout << endl; } */ // delete oscProb; }

PIDSpectrum::ClassDef (  PIDSpectrum  , 9  )  [private]

vector< TCanvas * > PIDSpectrum::DrawEnergySlices (  const NC::OscProb::OscPars *  coords, TH1 *  systratios, std::string  suffix )  const

Draw slices of the energy spectrum for a given set of coords. systratios is the histogram of ratios from the systematics interpolation, which are applied only to the total spectrum. This is slightly wrong, since the components change with systematics too, but in a non-trivial way (eg, CC background should only affect the CC component, not the others). This way they're all equally wrong, instead of being wrong in different ways. Definition at line 542 of file PIDSpectrum.cxx. References DrawSlices(). Referenced by NCExtrapolationPID::WriteSpectra(). { return DrawSlices("energy", coords, systratios, suffix); }

vector< TCanvas * > PIDSpectrum::DrawNearSpectra (  TH1 *  systratios, std::string  suffix )  const

Draw the near spectra onto a canvas. Definition at line 702 of file PIDSpectrum.cxx. References NC::Utility::CloneFast(), fBeamInfo, fNearData, fNearMC, NCBeam::Info::GetDescription(), and ProjectHist(). Referenced by NCExtrapolationPID::WriteSpectra(). { vector<TCanvas*> ret; TString canvName("nearSpectra "+fBeamInfo.GetDescription()); canvName+=suffix.c_str(); TCanvas* c = new TCanvas(canvName, "Near detector spectra"); c->Divide(2,2); c->cd(1); TH2F* nearMCShifted=(TH2F*)NC::Utility::CloneFast(fNearMC); nearMCShifted->SetName("nearMCShifted"); if (systratios) nearMCShifted->Multiply(systratios); nearMCShifted->Draw("colz"); c->cd(2); fNearData->Draw("colz"); c->cd(3); TH2F* hRatio=(TH2F*)NC::Utility::CloneFast(fNearData); hRatio->SetName("hRatio"); hRatio->Divide(nearMCShifted); hRatio->Draw("colz"); ret.push_back(c); nearMCShifted->SetLineColor(kRed); fNearData->SetLineColor(kBlack); // 1D spectra TCanvas* cE=new TCanvas(canvName+"energy_1D", "Near detector 1D spectra"); cE->cd(); ProjectHist(nearMCShifted, "energy", "nearMC")->DrawCopy(); ProjectHist(fNearData, "energy", "nearData")->DrawCopy("same e"); ret.push_back(cE); TCanvas* cPID=new TCanvas(canvName+"pid_1D", "Near detector 1D spectra"); cPID->cd(); ProjectHist(nearMCShifted, "pid", "nearMC")->DrawCopy(); ProjectHist(fNearData, "pid", "nearData")->DrawCopy("same e"); ret.push_back(cPID); return ret; }

vector< TCanvas * > PIDSpectrum::DrawPIDSlices (  const NC::OscProb::OscPars *  coords, TH1 *  systratios, std::string  suffix )  const

Draw slices of the PID spectrum for a given set of coords. See DrawEnergySlices for a comment on systratios Definition at line 551 of file PIDSpectrum.cxx. References DrawSlices(). Referenced by NCExtrapolationPID::WriteSpectra(). { return DrawSlices("pid", coords, systratios, suffix); }

vector< TCanvas * > PIDSpectrum::DrawSlices (  std::string  variable, const NC::OscProb::OscPars *  coords, TH1 *  systratios, std::string  suffix )  const [private]

Draw the energy or pid slices onto a canvas and return it. Definition at line 573 of file PIDSpectrum.cxx. References EvtTypeAsString(), fBeamInfo, fData, fDrawColors, fEvtTypes, fNumPIDBins, Form(), NCBeam::Info::GetDescription(), GetOnePredicted(), GetPredicted(), and ProjectHist(). Referenced by DrawEnergySlices(), and DrawPIDSlices(). { // Are we drawing energy slices? (If not, we're drawing PID slices) assert(variable == "energy" || variable == "pid"); bool energySlice=(variable=="energy"); vector<TCanvas*> v; // Canvas for the energy-integrated-over-PID plot (or vice versa) // ROOT is really not happy having more than one canvas with the same name // Keep a counter so we can keep the canvases from separate calls to this // function distinct. static int canvN=0; TCanvas* cTotal = new TCanvas((variable+"Total_"+suffix+fBeamInfo.GetDescription()+ Form("_%d", canvN))); // Canvas for the actual slices TCanvas* cSlices = new TCanvas((variable+"Slices_"+suffix+fBeamInfo.GetDescription()+ Form("_%d", canvN))); ++canvN; v.push_back(cTotal); v.push_back(cSlices); cSlices->cd(); int nPlots = energySlice ? fNumPIDBins : kNumEnergyBinsFar; int nDiv = (int)TMath::Ceil(TMath::Sqrt(nPlots)); // Divide with more rows than columns if possible, so we get wider graphs if ((nDiv-1)*nDiv>=nPlots) cSlices->Divide(nDiv-1, nDiv); else if ((nDiv+1)*(nDiv-1)>=nPlots) cSlices->Divide(nDiv-1, nDiv+1); else cSlices->Divide(nDiv, nDiv); TH2F* hTotal=GetPredicted(coords, "total", false); if (systratios) hTotal->Multiply(systratios); // In principle we might be being called with parameters that don't // correspond to the best fit, but this'll do for now hTotal->SetTitle("Best fit"); NC::OscProb::NoOscillations oscParsNoOsc; TH2F* hNoOsc=GetPredicted(&oscParsNoOsc, "noosc", false); hNoOsc->SetTitle("No Oscillations"); // Hold the components in a map map<evtType, TH2F*> cmpHists; for (vector<evtType>::const_iterator it=fEvtTypes.begin(); it!=fEvtTypes.end(); ++it) { // Drawing everything except NC if (*it==kNC) continue; cmpHists[*it]=GetOnePredicted(*it, coords, (variable+"Slices_"+suffix).c_str()); cmpHists[*it]->SetLineColor(fDrawColors[*it]); } hTotal->SetLineColor(kRed); hNoOsc->SetLineColor(kRed); hNoOsc->SetLineStyle(2); fData->SetLineColor(kBlack); // Stupid ROOT and its stupid global nonsense. gStyle->SetOptStat(0); // Draw the spectrum integrated over the other variable cTotal->cd(); TH1F* noOsc1D=(TH1F*)ProjectHist(hNoOsc, variable, "NoOsc")->DrawCopy(); noOsc1D->SetLineStyle(2); if (energySlice) noOsc1D->GetXaxis()->SetRangeUser(0,19.9); ProjectHist(hTotal, variable, "Total")->DrawCopy("same"); for (vector<evtType>::const_iterator it=fEvtTypes.begin(); it!=fEvtTypes.end(); ++it) { // Drawing everything except NC if (*it==kNC) continue; ProjectHist(cmpHists[*it], variable, string("slice")+EvtTypeAsString(*it))->DrawCopy("same"); } ProjectHist(fData, variable, "Data")->DrawCopy("same e"); ((TPad*)gPad)->BuildLegend()->SetFillColor(kWhite); if (!energySlice && noOsc1D->GetEntries()) ((TPad*)gPad)->SetLogy(); // Draw the individual slices for (int i=1; i<nPlots+1; ++i) { cSlices->cd(i); TH1F* noOsc1DSlice=(TH1F*)ProjectHist(hNoOsc, variable, "NoOsc", i, i)->DrawCopy(); noOsc1DSlice->SetLineStyle(2); if (energySlice) noOsc1DSlice->GetXaxis()->SetRangeUser(0,19.9); ProjectHist(hTotal, variable, "Total", i, i)->DrawCopy("same"); for(unsigned int n = 0; n < fEvtTypes.size(); ++n){ // Drawing everything except NC if (fEvtTypes[n]==kNC) continue; ProjectHist(cmpHists[fEvtTypes[n]], variable, EvtTypeAsString(fEvtTypes[n]), i, i)->DrawCopy("same"); } ProjectHist(fData, variable, "Data", i, i)->DrawCopy("same e"); if (!energySlice && noOsc1DSlice->GetEntries()) ((TPad*)gPad)->SetLogy(); } return v; }

TCanvas * PIDSpectrum::DrawSpectrum (  const NC::OscProb::OscPars *  coords, TH1 *  systratios, std::string  suffix = "" )  const

Draw the 2d spectrum for the given coords . The returned canvas contains the predicted spectrum at the given point, the data, and the ratio. Parameters: coords  Oscillate the spectrum with these parameters systratios  The ratios from the systematics interpolator to apply to the MC suffix  Suffix to add to the canvas name, to distinguish it from any others that may be drawn Definition at line 501 of file PIDSpectrum.cxx. References NC::Utility::CloneFast(), fCmpTrue, fData, Nav::GetName(), GetPredicted(), and kNC. Referenced by NCExtrapolationPID::WriteSpectra(). { TH2F* hPred=GetPredicted(coords, "total", false); if (systratios) hPred->Multiply(systratios); string canName("can_"); canName+=GetName(); canName+=suffix; TCanvas *c = new TCanvas(canName.c_str(), GetTitle()); c->Divide(2,2); c->cd(1); hPred->DrawCopy("colz"); TLatex* la=new TLatex; la->SetNDC(); la->DrawLatex(0.2, 0.95, "Predicted"); c->cd(2); fData->DrawCopy("colz"); la->DrawLatex(0.2, 0.95, "Data"); c->cd(3); TH2F* hRatio=(TH2F*)NC::Utility::CloneFast(fData); hRatio->SetName("ratio"); hRatio->Divide(hPred); hRatio->Draw("colz"); la->DrawLatex(0.2, 0.95, "Ratio"); c->cd(4); fCmpTrue.find(kNC)->second->DrawCopy("colz"); la->DrawLatex(0.2, 0.95, "True NC true E"); // delete hPred; return c; }

std::vector<PIDSpectrum::evtType> PIDSpectrum::EventTypeList (   )  const [inline]

The list of event types available. Definition at line 65 of file PIDSpectrum.h. Referenced by NCExtrapolationPID::WriteSpectra(). { return fEvtTypes; }

string PIDSpectrum::EvtTypeAsString (  evtType  t  )  const [protected]

Return event type as string. Definition at line 479 of file PIDSpectrum.cxx. References kBeamNue, kNC, kNumu, kNutau, and kOscNue. Referenced by DrawSlices(), and GetOnePredicted(). { switch (t) { case kNC: return "NC"; case kNumu: return "NuMu"; case kBeamNue: return "BeamNuE"; case kOscNue: return "OscNuE"; case kNutau: return "NuTau"; default: cerr << "EvtTypeAsString got unknown evtType " << t << endl; return "Unknown"; } }

void PIDSpectrum::FarNearCorrect (  TH2F *  fdPred, const TH2F *  ndData, const TH1 *  ndMC )  const [private]

Do the Far/Near correction, ie, multiply fdPred by ndData over ndMC. Definition at line 743 of file PIDSpectrum.cxx. Referenced by GetPredicted(), and GetPseudoNCCCSpectrum(). { // TODO: Replace the data and MC spectra with just the ratio, which // we'll only need to calculate once // Downcast ndMC - ugh. This should be OK because that histogram came // from us initially. TH2F* ndMC = (TH2F*)ndMC1D; // Use the now-standard trick of accessing bin contents directly to // avoid unnecessary bounds-checking for (int i=0; i<fdPred->fN; ++i) { // Don't divide by zero. This should really be left to crash, but // I'm too lazy at the moment if (ndMC->fArray[i]) fdPred->fArray[i]*=ndData->fArray[i]/ndMC->fArray[i]; } }

void PIDSpectrum::FillData (  Detector::Detector_t  det, const ANtpTruthInfoBeam *  t, const ANtpRecoInfo *  r, const ANtpAnalysisInfo *  ana, double  E, double  scale = 1 )  [virtual]

Add an event to the data histogram. Definition at line 410 of file PIDSpectrum.cxx. References det, fData, fNearData, kNumEnergyBinsFar, min, PIDWithinRange(), and ANtpAnalysisInfo::separationParameter. { // Cap energy at 120 GeV // The -0.01 is to make sure the event goes in the 20-120 GeV // bin and not the overflow bin above E=std::min(E,kEnergyBinsFar[kNumEnergyBinsFar]-0.01); double pidVal=PIDWithinRange(ana->separationParameter); if (det==Detector::kFar) fData->Fill(E, pidVal, scale); else if (det==Detector::kNear) fNearData->Fill(E, pidVal, scale); else { cout << "Unknown detector " << det << ". Bailing" << endl; exit(1); } }

void PIDSpectrum::FillMC (  Detector::Detector_t  det, const ANtpTruthInfoBeam *  t, const ANtpRecoInfo *  r, const ANtpAnalysisInfo *  ana, double  E, double  scale )  [virtual]

Add an event to the appropriate MC histogram. Definition at line 382 of file PIDSpectrum.cxx. References det, evtType, fCmpTrue, fNearMC, fTrueToReco, GetEvtType(), kNumEnergyBinsFar, min, ANtpTruthInfo::nuEnergy, PIDWithinRange(), and ANtpAnalysisInfo::separationParameter. { evtType evt(GetEvtType(t)); double pidVal=PIDWithinRange(ana->separationParameter); // Cap energy at 120 GeV // The -0.01 is to make sure the event goes in the 20-120 GeV // bin and not the overflow bin above E=std::min(E,kEnergyBinsFar[kNumEnergyBinsFar]-0.01); if (det==Detector::kFar) { fCmpTrue[evt]->Fill(t->nuEnergy, pidVal, scale); fTrueToReco[evt]->Fill(E, pidVal, t->nuEnergy, scale); } else if (det==Detector::kNear) { fNearMC->Fill(E, pidVal, scale); } else { cout << "Unknown detector " << det << ". Bailing" << endl; exit(1); } }

NCBeam::Info PIDSpectrum::GetBeamInfo (   )  [inline]

Get the beam index for this spectrum. Definition at line 154 of file PIDSpectrum.h. Referenced by NCExtrapolationPID::ReadMCSpectra(), and NCExtrapolationPID::WriteSpectra(). { return fBeamInfo; }

TH2F* PIDSpectrum::GetDataHist (   )  const [inline]

Get the data histogram for this Spectrum. Definition at line 124 of file PIDSpectrum.h. Referenced by NCExtrapolationPID::WriteResources(). { return fData; }

double PIDSpectrum::GetEventWeight (  evtType  t, double  trueE, const NC::OscProb::OscPars *  coords )  const [protected]

Return the oscillation weight for event of type t with true energy trueE at oscillation parameters coords. Definition at line 297 of file PIDSpectrum.cxx. References kBeamNue, kNC, kNumu, kNutau, kOscNue, and NC::OscProb::OscPars::TransitionProbability(). Referenced by GetOscWeightHist(). { switch (t) { case PIDSpectrum::kNC: // NCExtrapolationModule only passes NC events from the nominal // files - NC events from the tau and electron files are // ignored. The rationale is that there's more stats in the // nominal files than the others. // This is the reason for the bizarre interaction type argument here. return coords->TransitionProbability(NCType::kNuMuToNuMu, NCType::kNC, NCType::kBaseLineFar, trueE); break; case PIDSpectrum::kNumu: return coords->TransitionProbability(NCType::kNuMuToNuMu, NCType::kCC, NCType::kBaseLineFar, trueE); break; case PIDSpectrum::kBeamNue: // Assume beam nue's don't oscillate. // TODO: Do this right // Apparently this is what's done in other places in NCUtils, so // maybe it can stay this way return 1.; break; case PIDSpectrum::kOscNue: return coords->TransitionProbability(NCType::kNuMuToNuE, NCType::kCC, NCType::kBaseLineFar, trueE); break; case PIDSpectrum::kNutau: return coords->TransitionProbability(NCType::kNuMuToNuTau, NCType::kCC, NCType::kBaseLineFar, trueE); break; default: // Something went wrong cerr << "GetOnePredicted got unknown evtType " << t << endl; assert(0); } }

PIDSpectrum::evtType PIDSpectrum::GetEvtType (  const ANtpTruthInfoBeam *  t  )  const [protected]

Get the event type for a given truth. TODO: This should go away and be replaced with however NCUtils does it Definition at line 464 of file PIDSpectrum.cxx. References abs(), evtType, ANtpTruthInfo::interactionType, ANtpTruthInfoBeam::nonOscNuFlavor, and ANtpTruthInfo::nuFlavor. Referenced by FillMC(). { if (t->interactionType==0) return kNC; else if (abs(t->nuFlavor)==14) return kNumu; else if (abs(t->nuFlavor)==12 && abs(t->nonOscNuFlavor)==12) return kBeamNue; else if (abs(t->nuFlavor)==12 && abs(t->nonOscNuFlavor)!=12) return kOscNue; else if (abs(t->nuFlavor)==16) return kNutau; else { cerr << "Unknown particle type " << t->nuFlavor << endl; assert(0); } }

const TH2F* PIDSpectrum::GetNearMC (   )  const [inline]

Get the near detector Monte Carlo spectrum. Definition at line 81 of file PIDSpectrum.h. {return fNearMC;}

int PIDSpectrum::GetNPIDBins (   )  const [inline]

Return the number of bins in PID. Definition at line 197 of file PIDSpectrum.h. Referenced by NCExtrapolationPID::ReadMCSpectra(). { return fNumPIDBins; }

TH2F * PIDSpectrum::GetOnePredicted (  evtType  t, const NC::OscProb::OscPars *  coords, const char *  histname = 0 )  const

Get the predicted histogram for one event type for the given coords. Parameters: t  The event type to evaluate for coords  The oscillation/systematic parameters to evaluate at histname  The name to give the histogram Definition at line 159 of file PIDSpectrum.cxx. References AddOnePredicted(), EvtTypeAsString(), Form(), and Nav::GetName(). Referenced by DrawSlices(), GetPseudoNCCCSpectrum(), and NCExtrapolationPID::WriteSpectra(). { // This spectrum string tmp = GetName(); tmp+="_pred"; tmp+=histname ? histname : ""; tmp+=EvtTypeAsString(t); TH2F* onePred=new TH2F(tmp.c_str(), Form("Predicted %s", EvtTypeAsString(t).c_str()), kNumEnergyBinsFar, kEnergyBinsFar, fNumPIDBins, fPIDBins); AddOnePredicted(onePred, t, coords); return onePred; }

TH1F * PIDSpectrum::GetOscWeightHist (  evtType  t, const NC::OscProb::OscPars *  coords, int  nBins, Double_t *  binEdges )  const [protected, virtual]

Get the histogram of oscillation weights in true energy for the given event type. Definition at line 348 of file PIDSpectrum.cxx. References GetEventWeight(). Referenced by AddOnePredicted(). { // Oscillation probabilities for this event type as a fn of energy static TH1F* oscProb=new TH1F("oscProb", "Oscillation probability", nBins, binEdges); // oscProb->Reset(); // oscProb->SetDirectory(0); // The "<=" means we include the overflow bin for (int i=1; i<=oscProb->GetNbinsX()+1; ++i) { // Poor man's integral over the bin double Emin=oscProb->GetXaxis()->GetBinLowEdge(i); double Emax=oscProb->GetXaxis()->GetBinUpEdge(i); double trueE=0.5*(Emin+Emax); oscProb->fArray[i] = GetEventWeight(t, trueE, coords); } return oscProb; }

TH2F * PIDSpectrum::GetPredicted (  const NC::OscProb::OscPars *  coords, const char *  histname, bool  doFN, TH1 *  ndSpectrum = 0 )  const

Get the predicted total histogram for the given coords. Parameters: coords  The oscillation/systematic parameters to evaluate at histname  The name to give the histogram doFN  Whether to do the Far-Near correction ndSpectrum  Optional spectrum to use for Far-Near correction Definition at line 129 of file PIDSpectrum.cxx. References AddOnePredicted(), FarNearCorrect(), fEvtTypes, fNearData, fNearMC, fNumPIDBins, fPIDBins, and kNumEnergyBinsFar. Referenced by DrawSlices(), DrawSpectrum(), and GetPseudoNCCCSpectrum(). { TH2F* pred=new TH2F(histname, "Predicted total", kNumEnergyBinsFar, kEnergyBinsFar, fNumPIDBins, fPIDBins); pred->SetDirectory(0); for (vector<evtType>::const_iterator it=fEvtTypes.begin(); it!=fEvtTypes.end(); ++it) { AddOnePredicted(pred, *it, coords); } // Because we set the bins manually using fArray, the number of // entries isn't calculated. This results in the histogram not // drawing properly, so set the number of entries manually to // something reasonable pred->SetEntries(int(pred->Integral())); // Multiply the predicted spectrum by the ND data/MC ratio if (doFN) FarNearCorrect(pred, fNearData, nearSpectrum ? nearSpectrum : fNearMC); return pred; }

TH1D * PIDSpectrum::GetPseudoNCCCComponent (  TH2F *  h, NCType::EEventType  nccc )  const

Return the 1D part of h that is type nccc in the manner of GetPseudoNCCCSpectrum. Definition at line 787 of file PIDSpectrum.cxx. References fCutValue. Referenced by GetPseudoNCCCSpectrum(), and NCExtrapolationPID::WriteResources(). { // Find the bin where we should cut to get "NC" one side and "CC" the other int cutBin=h->GetYaxis()->FindBin(fCutValue); TH1D* spec1d; if(nccc==NCType::kNC) spec1d=h->ProjectionX("pseudonc", 0, cutBin); else spec1d=h->ProjectionX("pseudocc", cutBin+1, -1); // Need to SetDirectory(0) because of a throughly idiotic behaviour // of TH2::ProjectionX, which will overwrite an existing histogram // if the name you request already exists in the current // directory. Idiots. spec1d->SetDirectory(0); return spec1d; }

TH1D * PIDSpectrum::GetPseudoNCCCSpectrum (  NCType::EEventType  nccc, evtType  t, const NC::OscProb::OscPars *  coords )  const

Return the "NC" or "CC" spectrum for true event type t at coords. This is a helper function for the benefit of filling the NCBeam plots, which assume an NC and CC spectrum. We can get something similar by taking events below and above the cut value as NC/CC respectively. It's a slight lie, but makes plots that people can understand. Definition at line 765 of file PIDSpectrum.cxx. References FarNearCorrect(), fNearData, fNearMC, GetOnePredicted(), GetPredicted(), and GetPseudoNCCCComponent(). Referenced by NCExtrapolationPID::WriteResources(). { TH2F* spec2d; if(t==kAll){ // The "true" makes it do far/near correction, so we don't need to do it manually //NC::OscProb::OscPars* coordsNoOsc=new NC::OscProb::NoOscillations; spec2d=GetPredicted(coords, "", true); } else{ spec2d=GetOnePredicted(t, coords); spec2d->SetDirectory(0); // Should make this configurable. Too lazy right now FarNearCorrect(spec2d, fNearData, fNearMC); } return GetPseudoNCCCComponent(spec2d, nccc); }

std::map<PIDSpectrum::evtType, TH2F*> PIDSpectrum::GetTrueEMap (   )  const [inline]

Get the true energy spectra for each event type. Definition at line 151 of file PIDSpectrum.h. {return fCmpTrue;};

std::map<PIDSpectrum::evtType, TH3F*> PIDSpectrum::GetTrueToRecoMap (   )  const [inline]

Get the true-to-reco maps for each event type. Definition at line 148 of file PIDSpectrum.h. {return fTrueToReco;} ;

Double_t * PIDSpectrum::MultiplyBins (  int  factor, int  nBinsOrig, const Double_t *  origBins )  const [private]

Make a binning based on origBins, but factor times finer. This is used to create the true energy binning for the true and true-to-reco hists, based on the NCUtils binning Definition at line 112 of file PIDSpectrum.cxx. Referenced by AddOnePredicted(). { Double_t* newBins=new Double_t[(nBinsOrig+1)*factor+2]; for (int iOrig=0; iOrig<nBinsOrig+1; ++iOrig) { for (int j=0; j<factor; ++j) { newBins[iOrig*factor+j]=origBins[iOrig]+ float(j)/factor*(origBins[iOrig+1]-origBins[iOrig]); } } return newBins; }

void PIDSpectrum::NormalizeTrueToRecos (   )  [virtual]

Normalize the true to reco matrices so that the cells represent the probability that a neutrino with given true energy produces an event with given reco energy and PID value Definition at line 433 of file PIDSpectrum.cxx. References fEvtTypes, and fTrueToReco. { // The definition of normalized here is that any sum at fixed trueE and PID // over all the recoE's should come to one. for(vector<evtType>::iterator it = fEvtTypes.begin(); it != fEvtTypes.end(); ++it){ TH3F* t2r = fTrueToReco[*it]; const int maxE = t2r->GetNbinsZ()+1; // The "<=" means we include the overflow bin for(int iTrueE = 1; iTrueE <= maxE; ++iTrueE){ for(int iPID = 1; iPID <= fNumPIDBins; ++iPID){ float integral = 0; for(int iRecoE = 1; iRecoE <= kNumEnergyBinsFar; ++iRecoE){ integral += t2r->GetBinContent(iRecoE, iPID, iTrueE); } if(integral == 0) integral = 1; // Don't divide by zero for(int iRecoE = 1; iRecoE <= kNumEnergyBinsFar; ++iRecoE){ const float content = t2r->GetBinContent(iRecoE, iPID, iTrueE); t2r->SetBinContent(iRecoE, iPID, iTrueE, content/integral); } } // end for PID } // end for trueE } // end for evtType }

double PIDSpectrum::PIDWithinRange (  double  origPID  )  const [protected]

Bring origPID within the range [PIDmin, PIDmax ] so that nothing goes into the overflow bins. This is something of a hack: the ANN gives special values (-1 and 1.5) to certain classes of events, with the rest getting values between 0 and 1 (ish). This function allows us to include all the events with special values without having a load of empty bins Definition at line 373 of file PIDSpectrum.cxx. References fPIDmax, and fPIDmin. Referenced by FillData(), and FillMC(). { if(origPID>fPIDmax) return fPIDmax-1e-3; if(origPID<fPIDmin) return fPIDmin+1e-3; return origPID; }

TH1D * PIDSpectrum::ProjectHist (  TH2 *  h, std::string  axis, std::string  histname, int  lobin = 0, int  hibin = -1 )  const [private]

Draw the energy or pid projection of histogram h. Definition at line 683 of file PIDSpectrum.cxx. References NC::Utility::CloneFast(), and Plot::Format(). Referenced by DrawNearSpectra(), and DrawSlices(). { TH1D* ret; histname=axis+histname; if(lobin) histname+=TString::Format("%d", lobin); else histname+="Integ"; if (axis=="energy") ret = h->ProjectionX(histname.c_str(), lobin, hibin); else if (axis=="pid") ret = h->ProjectionY(histname.c_str(), lobin, hibin); else assert(0 && "unknown projection requested"); ret->SetDirectory(0); TH1D* clone=(TH1D*)NC::Utility::CloneFast(ret); clone->SetDirectory(0); return clone; }

void PIDSpectrum::ResetData (   )  [inline]

Clears the data histogram. Definition at line 200 of file PIDSpectrum.h. References fData. {fData->Reset();}

void PIDSpectrum::ResetMC (   )

Clears all the Monte Carlo histograms. Definition at line 560 of file PIDSpectrum.cxx. References fCmpTrue, fTrueToReco, and ANtpTruthInfoBeam::Reset(). { typedef map<evtType, TH2F*>::iterator it2; for(it2 it = fCmpTrue.begin(); it != fCmpTrue.end(); ++it) it->second->Reset(); typedef map<evtType, TH3F*>::iterator it3; for(it3 it = fTrueToReco.begin(); it != fTrueToReco.end(); ++it) it->second->Reset(); }

void PIDSpectrum::ResetNearData (   )  [inline]

Clears the ND data histogram. Definition at line 206 of file PIDSpectrum.h. References fNearData. {fNearData->Reset();}

void PIDSpectrum::ResetNearMC (   )  [inline]

Clears the ND Monte Carlo histogram. Definition at line 209 of file PIDSpectrum.h. References fNearMC. {fNearMC->Reset();}

---

Member Data Documentation

NCBeam::Info PIDSpectrum::fBeamInfo [protected]

The beam info for this spectrum. Definition at line 286 of file PIDSpectrum.h. Referenced by DrawNearSpectra(), and DrawSlices().

std::map<evtType, TH2F*> PIDSpectrum::fCmpTrue [protected]

The true energy components. Definition at line 249 of file PIDSpectrum.h. Referenced by AddOnePredicted(), DrawSpectrum(), FillMC(), and ResetMC().

double PIDSpectrum::fCutValue [protected]

The NC/CC cut value. Only used for the NCBeam plot filling,. Definition at line 293 of file PIDSpectrum.h. Referenced by GetPseudoNCCCComponent().

TH2F* PIDSpectrum::fData [protected]

The data spectrum. Definition at line 253 of file PIDSpectrum.h. Referenced by DrawSlices(), DrawSpectrum(), FillData(), PIDSpectrum(), and ResetData().

std::vector<int> PIDSpectrum::fDrawColors [protected]

Colours for the components of the selected spectrum. Definition at line 289 of file PIDSpectrum.h. Referenced by DrawSlices().

std::vector<evtType> PIDSpectrum::fEvtTypes [protected]

The available event types so we can loop over them. Definition at line 261 of file PIDSpectrum.h. Referenced by DrawSlices(), GetPredicted(), and NormalizeTrueToRecos().

TH2F* PIDSpectrum::fNearData [protected]

The near data. Definition at line 258 of file PIDSpectrum.h. Referenced by DrawNearSpectra(), FillData(), GetPredicted(), GetPseudoNCCCSpectrum(), PIDSpectrum(), and ResetNearData().

TH2F* PIDSpectrum::fNearMC [protected]

The near MC. Definition at line 256 of file PIDSpectrum.h. Referenced by DrawNearSpectra(), FillMC(), GetPredicted(), GetPseudoNCCCSpectrum(), PIDSpectrum(), and ResetNearMC().

int PIDSpectrum::fNumBinEdgesDontUseThis [protected]

The number of bin edges. Need this to be able to stream the pidBins array: we need to tell ROOT how large the array is, but nPIDBins won't do it, since the number of bin *edges* is one larger. Therefore this should always be equal to nPIDBins+1 Definition at line 274 of file PIDSpectrum.h.

int PIDSpectrum::fNumPIDBins [protected]

Definition at line 264 of file PIDSpectrum.h. Referenced by AddOnePredicted(), DrawSlices(), and GetPredicted().

Double_t* PIDSpectrum::fPIDBins [protected]

The bin edges in PID - needed as a member function for the case where we're emulating the FarNear method, as it's used in GetPredicted Definition at line 280 of file PIDSpectrum.h. Referenced by GetPredicted(), and PIDSpectrum().

double PIDSpectrum::fPIDmax [protected]

Definition at line 265 of file PIDSpectrum.h. Referenced by PIDWithinRange().

double PIDSpectrum::fPIDmin [protected]

Definition at line 265 of file PIDSpectrum.h. Referenced by PIDWithinRange().

int PIDSpectrum::fTrueBinFactor [protected]

How much finer to make the true energy binning than the reco energy binning. Definition at line 283 of file PIDSpectrum.h. Referenced by AddOnePredicted().

std::map<evtType, TH3F*> PIDSpectrum::fTrueToReco [protected]

The true-to-reco matrices. Definition at line 251 of file PIDSpectrum.h. Referenced by AddOnePredicted(), FillMC(), NormalizeTrueToRecos(), and ResetMC().

---

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

• PIDSpectrum.h
• PIDSpectrum.cxx

---