loon: PIDSpectrum Class Reference

Enumerator:

kNC
kNumu
kBeamNue
kOscNue
kNutau
kAll

Definition at line 62 of file PIDSpectrum.h.

00062 { kNC=0, kNumu=1, kBeamNue=2, kOscNue=3, kNutau=4, kAll=5 };

PIDSpectrum::PIDSpectrum ( ) [inline]

Default constructor: Needed for persisting to disk.

Definition at line 44 of file PIDSpectrum.h.

00045     : fData(0), fNearMC(0), fNearMCCC(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`,
		bool	_drawFDData = `false`
	)

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

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

Definition at line 58 of file PIDSpectrum.h.

References fPIDBins.

00058 { if (fPIDBins) delete[] fPIDBins; }

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 187 of file PIDSpectrum.cxx.

References fCmpTrue, fTrueToReco, GetOscWeightHist(), kNumEnergyBinsFar, and MultiplyBins().

Referenced by GetOnePredicted(), and GetPredicted().

00188 {
00189   // Set up the binning for true energy axes
00190   Double_t* trueEBins=MultiplyBins(fTrueBinFactor, kNumEnergyBinsFar,
00191                                    kEnergyBinsFar);
00192   const int nTrueEBins=kNumEnergyBinsFar*fTrueBinFactor;
00193 
00194   TH1F* oscProb=GetOscWeightHist(t, coords, nTrueEBins, trueEBins);
00195   delete[] trueEBins;
00196 
00197   /*
00198 
00199   // This is some debugging code that might turn out to be useful again one day
00200 
00201   if (oscProb->Integral()==0 && (t==kNumu || t==kNC)) {
00202     cout << "oscprob empty for type " << EvtTypeAsString(t) << endl;
00203     cout << "Coords: ";
00204     for (int i=0; i<(int)coords.size(); ++i) cout << coords[i] << " ";
00205     cout << endl;
00206   }
00207 
00208   static bool doneAlready=false;
00209 
00210   if (t==kNC && !doneAlready) { oscProb->Write(); doneAlready=true; }
00211 
00212   */
00213   // Convolve it all together
00214 
00215   // Need to use find() because map::operator[] isn't const
00216   TH3F* t2r=fTrueToReco.find(t)->second;
00217   TH2F* trueCmp=fCmpTrue.find(t)->second;
00218   /*
00219   if ((t==kNC || t==kNumu) && hit->second->Integral()==0) {
00220     cout << "t2r empty for type " << EvtTypeAsString(t) << endl;
00221     cout << "Coords: ";
00222     for (int i=0; i<(int)coords.size(); ++i) cout << coords[i] << " ";
00223     cout << endl;
00224   }
00225 
00226   if ((t==kNC || t==kNumu) && hit2->second->Integral()==0) {
00227     cout << "fCmpTrue empty for type " << EvtTypeAsString(t) << endl;
00228     cout << "Coords: ";
00229     for (int i=0; i<(int)coords.size(); ++i) cout << coords[i] << " ";
00230     cout << endl;
00231   }
00232   */
00233 
00234   // TH2 and TH3 convert to 1D bin numbers when asked for a bin by 2D
00235   // number. It turns out that this is rather slow (probably since
00236   // it's all implemented by the TH1 base class). Since I know I won't
00237   // be accessing past the end of the range, and I know what dimension
00238   // the histogram has, I can calculate this manually for (hopefully)
00239   // a speedup
00240 
00241   // GetBinContent does bounds-checking, which we don't need,
00242   // since the loop limits ensure that we only access real
00243   // bins. Instead, use the internal array fArray, which for
00244   // some reason is public, to access the bin values without
00245   // incurring the time penalty from bounds-checking
00246   //
00247   // Similarly SetBinContent can be replaced by a direct write to fArray
00248   // The disadvantage here is that bin errors etc won't be set, but we
00249   // don't need them here
00250 
00251   // The arrangement of t2r in memory in fArray is my trueE, PID, recoE
00252   // most significant first. By putting the loops in this same order
00253   // we consider every point in the order they are laid out in memory.
00254   //
00255   // onePred has the same arrangement, except it doesn't have a trueE direction
00256   //
00257   // It is important we loop through every single bin (including underflow
00258   // and overflow) or the counting of the index will go wrong.
00259 
00260   // Begin at the start of the true-to-reco array
00261   float* pt2r = t2r->fArray;
00262 
00263   for(int iTrueE = 0; iTrueE <= nTrueEBins+1; ++iTrueE){
00264     const float oscWeight = oscProb->fArray[iTrueE];
00265 
00266     // For each new true energy go back to the start of the output array
00267     float* pOnePred = onePred->fArray;
00268 
00269     for(int iPID=0; iPID <= fNumPIDBins+1; ++iPID){
00270       // The trueCmp histogram has its axes the wrong way round, so have
00271       // to do this calculation instead of being able to use the same trick
00272       // as for t2r and onePred
00273       const int trueE_PIDBin = iTrueE+(nTrueEBins+2)*iPID;
00274 
00275       const float pidWeight = trueCmp->fArray[trueE_PIDBin];
00276       const float oscPidWeight = oscWeight * pidWeight;
00277 
00278       for(int iRecoE = 0; iRecoE <= kNumEnergyBinsFar+1; ++iRecoE){
00279 
00280         const float fillValue = (*pt2r)*oscPidWeight;
00281 
00282         *pOnePred += fillValue;
00283 
00284         // Advance to the next entry in each array
00285         ++pt2r;
00286         ++pOnePred;
00287 
00288       } // end for iRecoE
00289     } // end for iPID
00290   } // end for iTrueE
00291   /*
00292   if ((t==kNC || t==kNumu) && onePred->Integral()<1e-6) {
00293     cout << "predicted spectrum empty for type " << EvtTypeAsString(t) << endl;
00294     cout << "Coords: ";
00295     for (int i=0; i<(int)coords.size(); ++i) cout << coords[i] << " ";
00296     cout << endl;
00297   }
00298   */
00299   //  delete oscProb;
00300 }

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

std::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.

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

Draw the near spectra onto a canvas.

std::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

std::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, along.

std::vector<TCanvas*> PIDSpectrum::DrawSlicesNear	(	std::string	variable,
		TH1 *	systratios,
		std::string	suffix
	)			const `[private]`

Draw the ND energy or pid slices onto a canvas and return it,.

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

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

The list of event types available.

Definition at line 65 of file PIDSpectrum.h.

References fEvtTypes.

00065 { return fEvtTypes; }

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

Return event type as string.

Definition at line 487 of file PIDSpectrum.cxx.

References kBeamNue, kNC, kNumu, kNutau, and kOscNue.

Referenced by GetOnePredicted().

00488 {
00489 
00490   switch (t) {
00491   case kNC:
00492     return "NC";
00493   case kNumu:
00494     return "NuMu";
00495   case kBeamNue:
00496     return "BeamNuE";
00497   case kOscNue:
00498     return "OscNuE";
00499   case kNutau:
00500     return "NuTau";
00501   default:
00502     cerr << "EvtTypeAsString got unknown evtType " << t << endl;
00503     return "Unknown";
00504   }
00505 }

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 850 of file PIDSpectrum.cxx.

Referenced by GetPredicted(), and GetPseudoNCCCSpectrumFar().

00853 {
00854   // TODO: Replace the data and MC spectra with just the ratio, which
00855   // we'll only need to calculate once
00856 
00857   // Downcast ndMC - ugh. This should be OK because that histogram came
00858   // from us initially.
00859   TH2F* ndMC = (TH2F*)ndMC1D;
00860 
00861   // Use the now-standard trick of accessing bin contents directly to
00862   // avoid unnecessary bounds-checking
00863   for (int i=0; i<fdPred->fN; ++i) {
00864     // Don't divide by zero. This should really be left to crash, but
00865     // I'm too lazy at the moment
00866     if (ndMC->fArray[i]) fdPred->fArray[i]*=ndData->fArray[i]/ndMC->fArray[i];
00867   }
00868 }

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 418 of file PIDSpectrum.cxx.

References fData, fNearData, Detector::kFar, Detector::kNear, kNumEnergyBinsFar, min, PIDWithinRange(), and ANtpAnalysisInfo::separationParameter.

00423 {
00424 
00425   // Cap energy at 120 GeV
00426   // The -0.01 is to make sure the event goes in the 20-120 GeV
00427   // bin and not the overflow bin above
00428   E=std::min(E,kEnergyBinsFar[kNumEnergyBinsFar]-0.01);
00429 
00430   double pidVal=PIDWithinRange(ana->separationParameter);
00431   if (det==Detector::kFar) fData->Fill(E, pidVal, scale);
00432   else if (det==Detector::kNear) fNearData->Fill(E, pidVal, scale);
00433   else {
00434     cout << "Unknown detector " << det << ". Bailing" << endl;
00435     exit(1);
00436   }
00437 }

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 389 of file PIDSpectrum.cxx.

References fCmpTrue, fNearMC, fNearMCCC, fTrueToReco, GetEvtType(), ANtpTruthInfo::interactionType, NCType::kCC, Detector::kFar, Detector::kNear, kNumEnergyBinsFar, min, ANtpTruthInfo::nuEnergy, PIDWithinRange(), and ANtpAnalysisInfo::separationParameter.

00394 {
00395   evtType evt(GetEvtType(t));
00396 
00397   double pidVal=PIDWithinRange(ana->separationParameter);
00398 
00399   // Cap energy at 120 GeV
00400   // The -0.01 is to make sure the event goes in the 20-120 GeV
00401   // bin and not the overflow bin above
00402   E=std::min(E,kEnergyBinsFar[kNumEnergyBinsFar]-0.01);
00403 
00404   if (det==Detector::kFar) {
00405     fCmpTrue[evt]->Fill(t->nuEnergy, pidVal, scale);
00406     fTrueToReco[evt]->Fill(E, pidVal, t->nuEnergy, scale);
00407   } else if (det==Detector::kNear) {
00408     fNearMC->Fill(E, pidVal, scale);
00409     if(t->interactionType==NCType::kCC) fNearMCCC->Fill(E, pidVal, scale);
00410   } else {
00411     cout << "Unknown detector " << det << ". Bailing" << endl;
00412     exit(1);
00413   }
00414 }

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

Get the beam index for this spectrum.

Definition at line 170 of file PIDSpectrum.h.

References fBeamInfo.

Referenced by NCExtrapolationPID::ReadMCSpectra().

00170 { return fBeamInfo; }

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

Get the data histogram for this Spectrum.

Definition at line 137 of file PIDSpectrum.h.

References fData.

Referenced by NCExtrapolationPID::WriteResources().

00137 { 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 304 of file PIDSpectrum.cxx.

References NCType::kBaseLineFar, kBeamNue, NCType::kCC, NCType::kNC, kNC, kNumu, NCType::kNuMuToNuE, NCType::kNuMuToNuMu, NCType::kNuMuToNuTau, kNutau, kOscNue, and NC::OscProb::OscPars::TransitionProbability().

Referenced by GetOscWeightHist().

00306 {
00307   switch (t) {
00308   case PIDSpectrum::kNC:
00309     // NCExtrapolationModule only passes NC events from the nominal
00310     // files - NC events from the tau and electron files are
00311     // ignored. The rationale is that there's more stats in the
00312     // nominal files than the others.
00313 
00314     // This is the reason for the bizarre interaction type argument here.
00315     return coords->TransitionProbability(NCType::kNuMuToNuMu,
00316                                          NCType::kNC,
00317                                          NCType::kBaseLineFar,
00318                                          trueE);
00319     break;
00320   case PIDSpectrum::kNumu:
00321     return coords->TransitionProbability(NCType::kNuMuToNuMu,
00322                                          NCType::kCC,
00323                                          NCType::kBaseLineFar,
00324                                          trueE);
00325     break;
00326   case PIDSpectrum::kBeamNue:
00327     // Assume beam nue's don't oscillate.
00328 
00329     // TODO: Do this right
00330     // Apparently this is what's done in other places in NCUtils, so
00331     // maybe it can stay this way
00332     return 1.;
00333     break;
00334   case PIDSpectrum::kOscNue:
00335     return coords->TransitionProbability(NCType::kNuMuToNuE,
00336                                          NCType::kCC,
00337                                          NCType::kBaseLineFar,
00338                                          trueE);
00339     break;
00340   case PIDSpectrum::kNutau:
00341     return coords->TransitionProbability(NCType::kNuMuToNuTau,
00342                                          NCType::kCC,
00343                                          NCType::kBaseLineFar,
00344                                          trueE);
00345     break;
00346   default:
00347       // Something went wrong
00348     cerr << "GetOnePredicted got unknown evtType " << t << endl;
00349     assert(0);
00350   }
00351 }

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 472 of file PIDSpectrum.cxx.

References ANtpTruthInfo::interactionType, kBeamNue, kNC, kNumu, kNutau, kOscNue, ANtpTruthInfoBeam::nonOscNuFlavor, and ANtpTruthInfo::nuFlavor.

Referenced by FillMC().

00473 {
00474   if (t->interactionType==0)     return kNC;
00475   else if (abs(t->nuFlavor)==14) return kNumu;
00476   else if (abs(t->nuFlavor)==12 && abs(t->nonOscNuFlavor)==12) return kBeamNue;
00477   else if (abs(t->nuFlavor)==12 && abs(t->nonOscNuFlavor)!=12) return kOscNue;
00478   else if (abs(t->nuFlavor)==16) return kNutau;
00479   else {
00480     cerr << "Unknown particle type " << t->nuFlavor << endl;
00481     assert(0);
00482   }
00483 }

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

Get the near data histogram for this Spectrum.

Definition at line 140 of file PIDSpectrum.h.

References fNearData.

Referenced by NCExtrapolationPID::WriteResources().

00140 { return fNearData; }

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

Get the near detector Monte Carlo spectrum.

Definition at line 81 of file PIDSpectrum.h.

References fNearMC.

00081 {return fNearMC;}

int PIDSpectrum::GetNPIDBins ( ) const [inline]

Return the number of bins in PID.

Definition at line 213 of file PIDSpectrum.h.

References fNumPIDBins.

Referenced by NCExtrapolationPID::ReadMCSpectra().

00213 { 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 166 of file PIDSpectrum.cxx.

References AddOnePredicted(), EvtTypeAsString(), Form(), fPIDBins, Nav::GetName(), and kNumEnergyBinsFar.

Referenced by GetPseudoNCCCSpectrumFar().

00169 {
00170   // This spectrum
00171   string tmp = GetName();
00172   tmp+="_pred";
00173   tmp+=histname ? histname : "";
00174   tmp+=EvtTypeAsString(t);
00175   TH2F* onePred=new TH2F(tmp.c_str(),
00176                          Form("Predicted %s", EvtTypeAsString(t).c_str()),
00177                          kNumEnergyBinsFar, kEnergyBinsFar,
00178                          fNumPIDBins, fPIDBins);
00179 
00180   AddOnePredicted(onePred, t, coords);
00181 
00182   return onePred;
00183 }

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 355 of file PIDSpectrum.cxx.

References GetEventWeight().

Referenced by AddOnePredicted().

00358 {
00359   // Oscillation probabilities for this event type as a fn of energy
00360   static TH1F* oscProb=new TH1F("oscProb", "Oscillation probability", nBins, binEdges);
00361 
00362   //  oscProb->Reset();
00363   //  oscProb->SetDirectory(0);
00364 
00365 
00366   // The "<=" means we include the overflow bin
00367   for (int i=1; i<=oscProb->GetNbinsX()+1; ++i) {
00368     // Poor man's integral over the bin
00369     double Emin=oscProb->GetXaxis()->GetBinLowEdge(i);
00370     double Emax=oscProb->GetXaxis()->GetBinUpEdge(i);
00371     double trueE=0.5*(Emin+Emax);
00372     oscProb->fArray[i] = GetEventWeight(t, trueE, coords);
00373   }
00374 
00375   return oscProb;
00376 }

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 136 of file PIDSpectrum.cxx.

References AddOnePredicted(), FarNearCorrect(), fEvtTypes, fNearData, fNearMC, fPIDBins, and kNumEnergyBinsFar.

Referenced by GetPseudoNCCCSpectrumFar().

00140 {
00141   TH2F* pred=new TH2F(histname, "Predicted total;Reco Energy;PID",
00142                       kNumEnergyBinsFar, kEnergyBinsFar,
00143                       fNumPIDBins, fPIDBins);
00144   pred->SetDirectory(0);
00145 
00146   for (vector<evtType>::const_iterator it=fEvtTypes.begin(); it!=fEvtTypes.end(); ++it) {
00147     AddOnePredicted(pred, *it, coords);
00148   }
00149 
00150   // Because we set the bins manually using fArray, the number of
00151   // entries isn't calculated. This results in the histogram not
00152   // drawing properly, so set the number of entries manually to
00153   // something reasonable
00154   pred->SetEntries(int(pred->Integral()));
00155 
00156   // Multiply the predicted spectrum by the ND data/MC ratio
00157   if (doFN)
00158     FarNearCorrect(pred, fNearData, nearSpectrum ? nearSpectrum : fNearMC);
00159 
00160   return pred;
00161 
00162 }

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 940 of file PIDSpectrum.cxx.

References MuELoss::e, fCutValue, and NCType::kNC.

Referenced by GetPseudoNCCCSpectrumFar(), GetPseudoNCCCSpectrumNear(), and NCExtrapolationPID::WriteResources().

00941 {
00942   // Find the bin where we should cut to get "NC" one side and "CC" the other
00943   int cutBin=h->GetYaxis()->FindBin(fCutValue-1e-5);
00944   TH1D* spec1d;
00945   if(nccc==NCType::kNC)
00946     spec1d=h->ProjectionX("pseudonc", 0, cutBin);
00947   else
00948     spec1d=h->ProjectionX("pseudocc", cutBin+1, -1);
00949 
00950   // Need to SetDirectory(0) because of a throughly idiotic behaviour
00951   // of TH2::ProjectionX, which will overwrite an existing histogram
00952   // if the name you request already exists in the current
00953   // directory. Idiots.
00954   spec1d->SetDirectory(0);
00955 
00956   return spec1d;
00957 }

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

Return the "NC" or "CC" spectrum for true event type t at coords in detector det.

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 872 of file PIDSpectrum.cxx.

References GetPseudoNCCCSpectrumFar(), GetPseudoNCCCSpectrumNear(), Detector::kFar, and Detector::kNear.

Referenced by NCExtrapolationPID::WriteResources().

00876 {
00877   switch(det){
00878   case Detector::kFar:
00879     return GetPseudoNCCCSpectrumFar(nccc, t, coords);
00880   case Detector::kNear:
00881     return GetPseudoNCCCSpectrumNear(nccc, t);
00882   default:
00883     assert(0 && "Unknown detector");
00884   }
00885 }

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

Implementation of GetPseudoNCCCSpectrum for the FD

Definition at line 889 of file PIDSpectrum.cxx.

References FarNearCorrect(), GetOnePredicted(), GetPredicted(), GetPseudoNCCCComponent(), and kAll.

Referenced by GetPseudoNCCCSpectrum().

00892 {
00893   TH2F* spec2d;
00894   if(t==kAll){
00895     // The "true" makes it do far/near correction, so we don't need to do it manually
00896     //NC::OscProb::OscPars* coordsNoOsc=new NC::OscProb::NoOscillations;
00897     spec2d=GetPredicted(coords, "", true);
00898   }
00899   else{
00900     spec2d=GetOnePredicted(t, coords);
00901     spec2d->SetDirectory(0);
00902     // Should make this configurable. Too lazy right now
00903     FarNearCorrect(spec2d, fNearData, fNearMC);
00904   }
00905 
00906   return GetPseudoNCCCComponent(spec2d, nccc);
00907 }

TH1D * PIDSpectrum::GetPseudoNCCCSpectrumNear	(	NCType::EEventType	nccc,
		evtType	t
	)			const

Implementation of GetPseudoNCCCSpectrum for the ND

Definition at line 911 of file PIDSpectrum.cxx.

References fNearMCCC, GetPseudoNCCCComponent(), kAll, kBeamNue, kNumu, kNutau, and kOscNue.

Referenced by GetPseudoNCCCSpectrum().

00913 {
00914   TH2F* spec2d;
00915   switch(t){
00916   case kAll:
00917     spec2d=fNearMC;
00918     break;
00919   case kNumu:
00920     spec2d=fNearMCCC;
00921     break;
00922   case kBeamNue:
00923   case kOscNue:
00924   case kNutau:
00925     // For these three cases just return an empty histogram of appropriate size
00926     spec2d=(TH2F*)fNearMC->Clone("emptyhist");
00927     spec2d->Reset("ICE");
00928     spec2d->SetDirectory(0);
00929     break;
00930   default:
00931     cerr << "Don't have component " << t << " in the ND. Bailing" << endl;
00932     exit(1);
00933   }
00934 
00935   return GetPseudoNCCCComponent(spec2d, nccc);
00936 }

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

Get the true energy spectra for each event type.

Definition at line 167 of file PIDSpectrum.h.

References fCmpTrue.

00167 {return fCmpTrue;};

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

Get the true-to-reco maps for each event type.

Definition at line 164 of file PIDSpectrum.h.

References fTrueToReco.

00164 {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 119 of file PIDSpectrum.cxx.

Referenced by AddOnePredicted().

00121 {
00122   Double_t* newBins=new Double_t[(nBinsOrig+1)*factor+2];
00123 
00124   for (int iOrig=0; iOrig<nBinsOrig+1; ++iOrig) {
00125     for (int j=0; j<factor; ++j) {
00126       newBins[iOrig*factor+j]=origBins[iOrig]+
00127         float(j)/factor*(origBins[iOrig+1]-origBins[iOrig]);
00128     }
00129   }
00130 
00131   return newBins;
00132 }

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 441 of file PIDSpectrum.cxx.

References fEvtTypes, fTrueToReco, and kNumEnergyBinsFar.

00442 {
00443   // The definition of normalized here is that any sum at fixed trueE and PID
00444   // over all the recoE's should come to one.
00445 
00446   for(vector<evtType>::iterator it = fEvtTypes.begin();
00447       it != fEvtTypes.end(); ++it){
00448     TH3F* t2r = fTrueToReco[*it];
00449     const int maxE = t2r->GetNbinsZ()+1;
00450     // The "<=" means we include the overflow bin
00451     for(int iTrueE = 1; iTrueE <= maxE; ++iTrueE){
00452       for(int iPID = 1; iPID <= fNumPIDBins; ++iPID){
00453 
00454         float integral = 0;
00455         for(int iRecoE = 1; iRecoE <= kNumEnergyBinsFar; ++iRecoE){
00456           integral += t2r->GetBinContent(iRecoE, iPID, iTrueE);
00457         }
00458         if(integral == 0) integral = 1; // Don't divide by zero
00459 
00460         for(int iRecoE = 1; iRecoE <= kNumEnergyBinsFar; ++iRecoE){
00461           const float content = t2r->GetBinContent(iRecoE, iPID, iTrueE);
00462           t2r->SetBinContent(iRecoE, iPID, iTrueE, content/integral);
00463         }
00464 
00465       } // end for PID
00466     } // end for trueE
00467   } // end for evtType
00468 }

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 380 of file PIDSpectrum.cxx.

References MuELoss::e, fPIDmax, and fPIDmin.

Referenced by FillData(), and FillMC().

00381 {
00382   if(origPID>fPIDmax) return fPIDmax-1e-3;
00383   if(origPID<fPIDmin) return fPIDmin+1e-3;
00384   return origPID;
00385 }

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.

void PIDSpectrum::ResetData ( ) [inline]

Clears the data histogram.

Definition at line 216 of file PIDSpectrum.h.

References fData.

00216 {fData->Reset();}

void PIDSpectrum::ResetMC ( )

Clears all the Monte Carlo histograms.

Definition at line 568 of file PIDSpectrum.cxx.

References fCmpTrue, and fTrueToReco.

00569 {
00570   typedef map<evtType, TH2F*>::iterator it2;
00571   for(it2 it = fCmpTrue.begin(); it != fCmpTrue.end(); ++it)
00572     it->second->Reset();
00573 
00574   typedef map<evtType, TH3F*>::iterator it3;
00575   for(it3 it = fTrueToReco.begin(); it != fTrueToReco.end(); ++it)
00576     it->second->Reset();
00577 }

void PIDSpectrum::ResetNearData ( ) [inline]

Clears the ND data histogram.

Definition at line 222 of file PIDSpectrum.h.

References fNearData.

00222 {fNearData->Reset();}

void PIDSpectrum::ResetNearMC ( ) [inline]

Clears the ND Monte Carlo histogram.

Definition at line 225 of file PIDSpectrum.h.

References fNearMC, and fNearMCCC.

00225 {fNearMC->Reset(); fNearMCCC->Reset();}

NCBeam::Info PIDSpectrum::fBeamInfo [protected]

The beam info for this spectrum.

Definition at line 304 of file PIDSpectrum.h.

Referenced by GetBeamInfo().

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

The true energy components.

Definition at line 265 of file PIDSpectrum.h.

Referenced by AddOnePredicted(), FillMC(), GetTrueEMap(), and ResetMC().

double PIDSpectrum::fCutValue [protected]

The NC/CC cut value. Only used for the NCBeam plot filling,.

Definition at line 311 of file PIDSpectrum.h.

Referenced by GetPseudoNCCCComponent().

TH2F* PIDSpectrum::fData [protected]

The data spectrum.

Definition at line 269 of file PIDSpectrum.h.

Referenced by FillData(), GetDataHist(), and ResetData().

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

Colours for the components of the selected spectrum.

Definition at line 307 of file PIDSpectrum.h.

bool PIDSpectrum::fDrawFDData [protected]

Whether to draw the FD data on plots (for blinding purposes).

Definition at line 314 of file PIDSpectrum.h.

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

The available event types so we can loop over them.

Definition at line 279 of file PIDSpectrum.h.

Referenced by EventTypeList(), GetPredicted(), and NormalizeTrueToRecos().

TH2F* PIDSpectrum::fNearData [protected]

The near data.

Definition at line 276 of file PIDSpectrum.h.

Referenced by FillData(), GetNearDataHist(), GetPredicted(), and ResetNearData().

TH2F* PIDSpectrum::fNearMC [protected]

The near MC.

Definition at line 272 of file PIDSpectrum.h.

Referenced by FillMC(), GetNearMC(), GetPredicted(), and ResetNearMC().

TH2F* PIDSpectrum::fNearMCCC [protected]

The near MC CC component. Not used in analysis, but want it for drawing.

Definition at line 274 of file PIDSpectrum.h.

Referenced by FillMC(), GetPseudoNCCCSpectrumNear(), 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 292 of file PIDSpectrum.h.

int PIDSpectrum::fNumPIDBins [protected]

Definition at line 282 of file PIDSpectrum.h.

Referenced by GetNPIDBins().

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 298 of file PIDSpectrum.h.

Referenced by GetOnePredicted(), GetPredicted(), and ~PIDSpectrum().

double PIDSpectrum::fPIDmax [protected]

Definition at line 283 of file PIDSpectrum.h.

Referenced by PIDWithinRange().

double PIDSpectrum::fPIDmin [protected]

Definition at line 283 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 301 of file PIDSpectrum.h.

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

The true-to-reco matrices.

Definition at line 267 of file PIDSpectrum.h.

Referenced by AddOnePredicted(), FillMC(), GetTrueToRecoMap(), NormalizeTrueToRecos(), and ResetMC().

PIDSpectrum Class Reference

Public Types

Public Member Functions

Protected Member Functions

Protected Attributes

Private Member Functions

Detailed Description

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Types
	kNC = 0
	kNumu = 1
	kBeamNue = 2
	kOscNue = 3
	kNutau = 4
	kAll = 5
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, bool _drawFDData=false)
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 (Detector::Detector_t det, NCType::EEventType nccc, evtType t, const NC::OscProb::OscPars *coords) const
	Return the "NC" or "CC" spectrum for true event type t at coords in detector det.
TH1D *	GetPseudoNCCCSpectrumFar (NCType::EEventType nccc, evtType t, const NC::OscProb::OscPars *coords) const
TH1D *	GetPseudoNCCCSpectrumNear (NCType::EEventType nccc, evtType t) const
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.
TH2F *	GetNearDataHist () const
	Get the near 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 *	fNearMCCC
	The near MC CC component. Not used in analysis, but want it for drawing.
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,.
bool	fDrawFDData
	Whether to draw the FD data on plots (for blinding purposes).
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, along.
std::vector< TCanvas * >	DrawSlicesNear (std::string variable, TH1 *systratios, std::string suffix) const
	Draw the ND 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)