RockFiller Class Reference

#include <RockFiller.h>

List of all members.

Public Member Functions

 RockFiller (const NuMatrix1D &baseTrueEnergy, const RockMatrix &baseRockMatrix)
 RockFiller (const NuMatrixCPT &baseTrueEnergy, const RockMatrix &baseRockMatrix)
 ~RockFiller ()
Int_t AddFile (TString filename)
Int_t AddFile (TString filename, Double_t POT)
Int_t GetEntries (void)
 Returns the number of rock events stored.
void Reduce (void)
 Attempt to reduce the array of events by looking for identical ones.
RockMatrixGetRockMatrix (const NuMMParameters &pars, NuMatrixCPT *True, int disappearanceModel, const bool display=false)
 Returns a true energy histogram of the stored events.
void SetPot (Double_t POT)
Double_t GetPot (void)
 Returns the internal POT being used.

Private Member Functions

void PrintProgress (Int_t entry, Int_t nEntries, Float_t percentage=20.)
 Prints the progress every percentage%.

Private Attributes

std::vector< RockEventfRockEvents
 A vector of stored Rock Events.
NuMatrix1D fBaseMatrix
RockMatrixfRockMatrix

Detailed Description

Holds a collection of rock events, and fills histograms at will

Definition at line 15 of file RockFiller.h.


Constructor & Destructor Documentation

RockFiller::RockFiller ( const NuMatrix1D baseTrueEnergy,
const RockMatrix baseRockMatrix 
)

Initialise a RockFiller class.

Parameters:
baseTrueEnergy The 'true energy' histogram to bin the events in
baseRockMatrix A rock matrix pointer for the 'base' matrix

Definition at line 31 of file RockFiller.cxx.

References NuMatrix::Complete(), NuMatrix1D::Complete(), NuMatrix::Copy(), fRockMatrix, Msg::kFatal, Msg::kWarning, and MSG.

00033           : fBaseMatrix(baseTrueEnergy),
00034             fRockMatrix(0)
00035 {
00036   // Ensure that the base matrix has a spectrum
00037   if (!baseTrueEnergy.Complete()) {
00038     MSG("RockFiller", Msg::kFatal) << "Have been passed a base matrix with no spectrum" << endl;
00039     throw runtime_error("Cannot proceed with empty base matrix");
00040   }
00041 
00042   // Throw a warning if the base rock matrix is non fully initialised
00043   if (!baseRockMatrix.Complete())
00044   {
00045     MSG("RockFiller",Msg::kWarning) << "Passed improperly initialised RockMatrix" << endl;
00046   }
00047 
00048   // Copy the base rock matrix
00049   fRockMatrix =  dynamic_cast<RockMatrix*>(baseRockMatrix.Copy());
00050 }

RockFiller::RockFiller ( const NuMatrixCPT baseTrueEnergy,
const RockMatrix baseRockMatrix 
)

Definition at line 52 of file RockFiller.cxx.

References base, NuMatrix::Complete(), NuMatrixCPT::Complete(), NuMatrix::Copy(), fBaseMatrix, fRockMatrix, Msg::kFatal, Msg::kWarning, MSG, and NuMatrixCPT::PQ().

00054           : fBaseMatrix(),
00055             fRockMatrix(0)
00056 {
00057   // Ensure that the base matrix has a spectrum
00058   if (!baseTrueEnergy.Complete()) {
00059     MSG("RockFiller", Msg::kFatal) << "Have been passed a base matrix with no spectrum" << endl;
00060     throw runtime_error("Cannot proceed with empty base matrix");
00061   }
00062 
00063   // Copy one of the histograms, for now
00064   const NuMatrix1D *base = dynamic_cast<const NuMatrix1D*>(baseTrueEnergy.PQ());
00065   if (!base) {
00066     MSG("RockFiller",Msg::kFatal) << "Cannot currently handle non-1d CPT matrices" << endl;
00067     throw runtime_error("Cannot proceed with non-1D base matrix");
00068   }
00069   fBaseMatrix = *base;
00070 
00071   // Throw a warning if the base rock matrix is non fully initialised
00072   if (!baseRockMatrix.Complete())
00073   {
00074     MSG("RockFiller",Msg::kWarning) << "Passed improperly initialised RockMatrix" << endl;
00075   }
00076 
00077   // Copy the base rock matrix
00078   fRockMatrix =  dynamic_cast<RockMatrix*>(baseRockMatrix.Copy());
00079 }

RockFiller::~RockFiller (  ) 

Definition at line 83 of file RockFiller.cxx.

References fRockMatrix.

00084 {
00085   if (fRockMatrix) delete fRockMatrix;
00086   fRockMatrix = 0;
00087 }


Member Function Documentation

Int_t RockFiller::AddFile ( TString  filename,
Double_t  POT 
)

Add a set of events to the event collection. This version allows explicit definition of POT

Parameters:
filename The RockEvent file to add
POT The POT in the file passed in
Returns:
The number of events added.

Definition at line 124 of file RockFiller.cxx.

References fBaseMatrix, fRockEvents, fRockMatrix, Msg::kError, Msg::kInfo, Msg::kWarning, MSG, RockEvent::neuEnMC, RockEvent::neuEnMCBin, RockMatrix::PrepareBins(), RockEvent::rw, and NuMatrix1D::Spectrum().

00125 {
00126   // We want to internally scale all events to the same POT.
00127   // Let's choose 1E20, because all analyses should be above this now!
00128   Double_t potscale = 1.0;
00129   // Check for POT count errors
00130   if (POT == 0.0) {
00131     MSG("RockFiller",Msg::kError) << "POT specified for rock muon file "
00132       << filename << " is zero. Not scaling (Assuming 1E20)" << endl;
00133   } else {
00134     potscale = 1.E20 / POT;
00135   }
00136 
00137 
00138   // Read out the trees via a TChain
00139   TChain r("r", "r");
00140   Int_t fcount = r.AddFile(filename, -1);
00141   // Check we actually read the file
00142   if (fcount == 0) {
00143     // We failed to read the file
00144     MSG("RockFiller", Msg::kWarning) << "Failed to read rock muon events from file " << filename << endl;
00145     return 0;
00146   }
00147 
00148   // Make a storage object and grab a branch
00149   RockEvent *nu = new RockEvent();
00150   TBranch *branch = r.GetBranch("r");
00151   branch->SetAddress(&nu);
00152 
00153   // Now, loop over all the rock events and copy them into a vector
00154   MSG("RockFiller",Msg::kInfo) << "Reading Rock Event file " << filename
00155       << " of " << POT << " POT" << endl;
00156   Int_t numEntries = r.GetEntries();
00157   for (Int_t entry=0; entry < numEntries; ++entry)
00158   {
00159     // Print progress
00160     // PrintProgress(entry, numEntries);
00161 
00162     // Get the event
00163     branch->GetEntry(entry);
00164 
00165     // Scale it to 1E20 POT
00166     nu->rw *= potscale;
00167 
00168     // Precalculate the bins for this NuEvent.
00169     // True energy is used by us, so do that here
00170     nu->neuEnMCBin = fBaseMatrix.Spectrum()->FindBin(nu->neuEnMC);
00171     // Delegate the other binning calculations to the RockMatrix object
00172     fRockMatrix->PrepareBins(nu);
00173 
00174     // // cout << "Energy bin is " << nu->neuEnMCBin << endl;
00175     // nu->trkEnBin    = fBaseMatrix.Spectrum()->FindBin(nu->trkEn);
00176     // nu->shwEnBin    = fBaseMatrix.Spectrum()->FindBin(nu->shwEn);
00177 
00178     // Store it
00179     fRockEvents.push_back(*nu);
00180   }
00181   // Return the number of entries just added
00182   return numEntries;
00183 }

Int_t RockFiller::AddFile ( TString  filename  ) 

Add a set of events to the event collection.

Parameters:
filename The RockEvent file to add
Returns:
The number of events added

Definition at line 92 of file RockFiller.cxx.

References exit(), Msg::kError, and MSG.

00093 {
00094   // Open the file and verify
00095   TFile tf(filename);
00096   if (!tf.IsOpen()) {
00097     MSG("RockFiller", Msg::kError) << "File " << filename << " does not exist or is bad." << endl;
00098     exit(1);
00099   }
00100 
00101   // Try to grab the POT histogram out of the file
00102   TH1F *pot = dynamic_cast<TH1F*>(tf.Get("hTotalPot"));
00103 
00104   // Grab the POT value
00105   Double_t filePot = 0.0;
00106   if (pot) {
00107     filePot = pot->Integral();
00108   } else {
00109     MSG("RockFiller", Msg::kError) << "File " << filename
00110       << " contains no POT info. Make an hTotalPot TH1D!" << endl;
00111 
00112   }
00113 
00114   // Close the file before moving on, as we open it again for AddFile
00115   tf.Close();
00116 
00117   // Pass on to AddFile, which will actually read the events
00118   return AddFile(filename, filePot);
00119 }

Int_t RockFiller::GetEntries ( void   )  [inline]

Returns the number of rock events stored.

Definition at line 40 of file RockFiller.h.

References fRockEvents.

00040 { return fRockEvents.size(); }

Double_t RockFiller::GetPot ( void   )  [inline]

Returns the internal POT being used.

Definition at line 70 of file RockFiller.h.

References fRockMatrix, and NuMatrix::GetPOT().

00070 { return fRockMatrix->GetPOT(); }

RockMatrix * RockFiller::GetRockMatrix ( const NuMMParameters pars,
NuMatrixCPT True,
int  disappearanceModel,
const bool  display = false 
)

Returns a true energy histogram of the stored events.

Get a filled RockMatrix object. This loops over all the stored events, correcting their weight according to the (optionally) passed in correction histogram.

Parameters:
TrueCorrection A true energy histogram to reweight events. CPT because it needs to reweight wrong sign correctly
Returns:
A filled copy of the rock matrix passed in the constructor Get a filled RockMatrix object. Using a parameters object, each event will be oscillated individually.

Definition at line 201 of file RockFiller.cxx.

References NuMatrixCPT::Complete(), NuMatrix::Copy(), NuMatrixCPT::DecayCC(), NuMatrixCPT::DecayMuToTau(), NuMatrixCPT::DecayNC(), NuMatrixCPT::Decohere(), fBaseMatrix, RockMatrix::Fill(), RockMatrix::FillJess(), fRockEvents, fRockMatrix, NuMatrix::GetPOT(), RockEvent::iaction, RockEvent::inu, Msg::kError, Msg::kInfo, Msg::kWarning, MAXMSG, MSG, RockEvent::neuEnMCBin, NuMatrixCPT::NQ(), NuMatrixCPT::Nu(), NuMatrixCPT::NuBar(), NuMatrixCPT::Oscillate(), NuMatrixCPT::PQ(), RockEvent::rw, NuMatrix::Scale(), NuMatrixCPT::SetValue(), and NuMatrix::Spectrum().

00202 {
00203   // MSG("RockFiller",Msg::kInfo) << "Oscillating per-event" << endl;
00204 
00205   if (!correction) {
00206     MSG("RockFiller",Msg::kWarning) << "No correction histogram being used." << endl;
00207   }
00208 
00209   // Validate the correction histogram
00210   if (correction && !correction->Complete())
00211   {
00212     MSG("RockFiller",Msg::kWarning) << "Rock muon spectrum correction passed in has no histogram. Ignoring." << endl;
00213 
00214     // Ignore the correction spectrum in this case
00215     correction = 0;
00216   }
00217 
00218   // Build base oscillation correction histograms
00219   NuMatrixCPT oschist(&fBaseMatrix, &fBaseMatrix);
00220   oschist.SetValue(1.0);
00221   NuMatrixCPT* oscNChist = 0; // No NChist means just use one everywher
00222   NuMatrixCPT* oscMuToTauhist = 0; // No MuToTau means use 1-osc
00223 
00224   switch(disMod){
00225   case 0:
00226     oschist.Oscillate(pars);
00227     break;
00228   case 1:
00229     oschist.DecayCC(pars);
00230 
00231     oscNChist = new NuMatrixCPT(&fBaseMatrix, &fBaseMatrix);
00232     oscNChist->SetValue(1.0);
00233     oscNChist->DecayNC(pars);
00234 
00235     oscMuToTauhist = new NuMatrixCPT(&fBaseMatrix, &fBaseMatrix);
00236     oscMuToTauhist->SetValue(1.0);
00237     oscMuToTauhist->DecayMuToTau(pars);
00238 
00239     break;
00240   case 2:
00241     oschist.Decohere(pars);
00242     break;
00243   default:
00244     assert(0 && "Badly configured disappearance model");
00245   }
00246 
00247   // If the base matrix has *no* POT, then we probably want to print a warning
00248   Double_t potScale = fRockMatrix->GetPOT() / 1.E20;
00249   if (potScale == 0.0) {
00250     MAXMSG("RockFiller",Msg::kWarning,5) << "Base matrix has no POT. You probably want to set some with SetPot()" << endl;
00251   }
00252 
00253   // Make a copy of the base RockMatrix
00254   // NuMatrix1D trkEn = fBaseMatrix;
00255   RockMatrix *rm =  dynamic_cast<RockMatrix*>(fRockMatrix->Copy());
00256   // Wipe this prediction out. Prevents adding on top of an existing one
00257   rm->Scale(0);
00258 
00259   // Now, loop over all events
00260   for (UInt_t i = 0; i < fRockEvents.size(); ++i)
00261   {
00262     // Get the event
00263     RockEvent &r = fRockEvents[i];
00264     // Calculate the reweighting for POT
00265     Double_t evtRw = r.rw * potScale;
00266 
00267     // Calculate any corrections if we are doing so
00268     // But, do so correctly for true charge
00269     if (correction) {
00270       if (r.inu > 0) {
00271         evtRw *= correction->NQ()->Spectrum()->GetBinContent(r.neuEnMCBin);
00272       } else if ( r.inu < 0) {
00273         evtRw *= correction->PQ()->Spectrum()->GetBinContent(r.neuEnMCBin);
00274       } else {
00275         MAXMSG("RockFiller",Msg::kError,10) << "Event " << i << "has inu=0! Cannot determine true charge-sign." << endl;
00276       }
00277     }
00278 
00279     // Need to fill jestograms before oscillations are applied. Otherwise taus
00280     // get zeroed out for no-oscillation case.
00281     rm->FillJess(r, evtRw);
00282 
00283     // Calculate the oscillation weighting, but only for true CC muons
00284     // Don't oscillate NC interactions
00285     if (r.iaction == 1) {
00286       // Oscillate differently depending on the particle type
00287       // Nue is deliberately ignored here - they will all be unoscillated
00288       if (r.inu == 14) {
00289         evtRw *= oschist.Nu()->Spectrum()->GetBinContent(r.neuEnMCBin);
00290       } else if (r.inu == -14) {
00291         evtRw *= oschist.NuBar()->Spectrum()->GetBinContent(r.neuEnMCBin);
00292       } else if (r.inu == 16) {
00293         // Inverse-oscillate tau neutrinos
00294         if(oscMuToTauhist){
00295           evtRw *= oscMuToTauhist->Nu()->Spectrum()->GetBinContent(r.neuEnMCBin);
00296         }
00297         else{
00298           evtRw *= 1. - oschist.Nu()->Spectrum()->GetBinContent(r.neuEnMCBin);
00299         }
00300       } else if (r.inu == -16) {
00301         // Inverse-oscillate tau neutrinos
00302         if(oscMuToTauhist){
00303           evtRw *= oscMuToTauhist->Nu()->Spectrum()->GetBinContent(r.neuEnMCBin);
00304         }
00305         else{
00306           evtRw *= 1. - oschist.NuBar()->Spectrum()->GetBinContent(r.neuEnMCBin);
00307         }
00308       } // end if antitau
00309     } else {
00310       // The particle is NC but also a tau - since we don't oscillate the
00311       // NC, ignore these completely to avoid duplication
00312       if (r.inu == 16 || r.inu == -16) {
00313         MAXMSG("RockFiller",Msg::kInfo,1)
00314           <<"Removing NC events from tau sample (nu.rw set to zero)"<<endl;
00315         evtRw = 0.0;
00316       }
00317 
00318       if(oscNChist) evtRw *= oscNChist->Nu()->Spectrum()->GetBinContent(r.neuEnMCBin);
00319     }
00320 
00321     // Finally, fill the rock matrix with this event
00322     rm->Fill(r, evtRw, display);
00323   }
00324 
00325   if(oscNChist) delete oscNChist;
00326   if(oscMuToTauhist) delete oscMuToTauhist;
00327 
00328   return rm;
00329 }

void RockFiller::PrintProgress ( Int_t  entry,
Int_t  nEntries,
Float_t  percentage = 20. 
) [private]

Prints the progress every percentage%.

Definition at line 188 of file RockFiller.cxx.

References Msg::kInfo, and MSG.

00189 {
00190   Float_t fract=ceil(nEntries/percentage);
00191   if (TMath::Ceil(((Float_t)entry)/fract)==((Float_t)entry)/fract){
00192     MSG("RockFiller", Msg::kInfo)
00193      << "Fraction of loop complete: "<<entry
00194      << "/"<<nEntries<<"  ("
00195      << (Int_t)(100.*entry/nEntries)<<"%)"<<endl;
00196   }
00197 }

void RockFiller::Reduce ( void   ) 

Attempt to reduce the array of events by looking for identical ones.

Definition at line 333 of file RockFiller.cxx.

References RockEvent::binNumber, fRockEvents, RockEvent::iaction, RockEvent::inu, Msg::kInfo, MSG, RockEvent::neuEnMCBin, and RockEvent::rw.

00334 {
00335   // Count of reduced events
00336   Int_t num = 0;
00337   // Count of skipped events
00338   Int_t numskipped = 0;
00339   // Map of [mcbin][fixbin][vector]
00340   map<Int_t, map<Int_t, vector<RockEvent> > > reduced;
00341 
00342   // Loop over every event we own
00343   for (UInt_t i = 0; i < fRockEvents.size(); ++i)
00344   {
00345     RockEvent &nu = fRockEvents[i];
00346 
00347     // Don't try to merge events with a bin number of -1
00348     if (nu.binNumber == -1) {
00349       numskipped++;
00350       num++;
00351       reduced[-1][-1].push_back(nu);
00352       continue;
00353     }
00354 
00355     // Look at every event in our so far 'reduced' sample for this category
00356     vector<RockEvent> &evtclass = reduced[nu.neuEnMCBin][nu.binNumber];
00357 
00358     // Find an event in here that matches!
00359     UInt_t j = 0;
00360     for (j = 0; j < evtclass.size(); ++j)
00361     {
00362       // Make sure the other variables match
00363       RockEvent &red = evtclass[j];
00364 
00365       if ( (nu.inu == red.inu) && (nu.iaction == red.iaction) )
00366       {
00367         break;
00368       }
00369     }
00370 
00371     // Did we find an event, or run out?
00372     // If we found one, merge with it. Otherwise, add as a new event!
00373     if (j < evtclass.size())
00374     {
00375       evtclass[j].rw += nu.rw;
00376     } else {
00377       // We didn't find a matching event. Add it!
00378       evtclass.push_back(nu);
00379       num++;
00380     }
00381   }
00382   // Print messages about how many events we reduced
00383   if (numskipped == num) {
00384     MSG("RockFiller",Msg::kInfo) << "Failed to reduce event set; No global bin numbers." << endl;
00385   } else {
00386     if (numskipped) {
00387       MSG("RockFiller",Msg::kInfo) << "Failed to reduce " << numskipped << " events due to no global bin numbers." << endl;
00388     }
00389     MSG("RockFiller",Msg::kInfo) << "Reduced " << fRockEvents.size() << " to " << num << " events." << endl;
00390   }
00391 
00392   // Now, copy these over to our storage object
00393   fRockEvents.clear();
00394   // Don't reserve, otherwise we tend to be grabbing very large chunks of memory (>200mb)
00395   // fRockEvents.reserve(num);
00396   // Loop over all MC bins
00397   map<Int_t, map<Int_t, vector<RockEvent> > >::iterator itMC = reduced.begin();
00398   for (; itMC != reduced.end(); ++itMC)
00399   {
00400     // Now, loop over all fitting bins
00401     map<Int_t, vector<RockEvent> >::iterator itBin = (*itMC).second.begin();
00402     for (; itBin != (*itMC).second.end(); ++itBin)
00403     {
00404       // Copy over this vector into the rockevents one
00405       vector<RockEvent> &evts = (*itBin).second;
00406       for (UInt_t i = 0; i < evts.size(); ++i)
00407       {
00408         fRockEvents.push_back(evts[i]);
00409       }
00410       // Now, clear this vector so that we clear space as we go
00411       evts.clear();
00412     } // reco bins
00413   } // mc bins
00414 
00415 }

void RockFiller::SetPot ( Double_t  POT  )  [inline]

Sets how many POT to use for the predictions This works by altering the POT count in the base histogram

Definition at line 68 of file RockFiller.h.

References fRockMatrix, and NuMatrix::ResetPOT().

00068 { fRockMatrix->ResetPOT(POT); }


Member Data Documentation

The base 'True Energy' matrix. This is used for corrections to the spectrum, so should be passed in for oscillations etc.

Definition at line 78 of file RockFiller.h.

Referenced by AddFile(), GetRockMatrix(), and RockFiller().

std::vector<RockEvent> RockFiller::fRockEvents [private]

A vector of stored Rock Events.

Definition at line 74 of file RockFiller.h.

Referenced by AddFile(), GetEntries(), GetRockMatrix(), and Reduce().

The base RockMatrix - the format that the events are returned as

Definition at line 81 of file RockFiller.h.

Referenced by AddFile(), GetPot(), GetRockMatrix(), RockFiller(), SetPot(), and ~RockFiller().


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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1