NuFCExperimentFactoryNSI Class Reference

Factory to generate pseudo-experiments for the NSI FC analysis. More...

#include <NuFCExperimentFactoryNSI.h>

List of all members.

Public Member Functions

 NuFCExperimentFactoryNSI (TRandom3 *rgen, TString nd_path, TString fd_path, TString tau_path)
 Construct the factory.
virtual ~NuFCExperimentFactoryNSI ()
const std::string & GetName () const
void GenerateNewExperiment (NuMMRunNSINu *mmRunSource, const NuMMParameters &mmPars, bool reroll=true)
 Generates a new experiment.
void GenerateNewExperiment (NuMMRunNSINubar *mmRunSource, const NuMMParameters &mmPars, bool reroll=true)
 Generates a new experiment.
const NuMatrixSpectrumGetSpectrum () const
 Get a reference to the FD Spectrum.
const NuMatrixSpectrumGetNewSpectrum () const
 Get a new copy (pointer) to the FD Spectrum.
const NuMatrixSpectrumGetSpectrumND () const
 Get a reference to the ND Spectrum.
const NuMatrixSpectrumGetNewSpectrumND () const
 Get a new copy (pointer) to the ND Spectrum.
void SetGausSyst (bool _gs)
 Set to true to have gaussian distributed syst values and false to have them be uniform.
void SetSystNorm (double size)
 Normalization systematic. Shift value is 1 +/- size. 0 = off.
void SetSystCurv (double size)
 Absolute Track Energy from curvature far systematic. Shift value is 1 +/- size. 0 = off.
void SetSystRelCurv (double size)
 Relative Track Energy from curvature far systematic. Shift value is 1 +/- size. 0 = off.
void SetSystRange (double size)
 Track Energy from range systematic. Shift value is 1 +/- size. 0 = off.
void SetSystShower (double size)
 Shower Energy overall systematic. Shift value is +/- size sigma. 0 = off.
void SetSystRelShower (double size)
 Relative shower energy (far only). Shift value is 1 +/- size. 0 = off.
bool IsDP (const NuFCEvent *nu) const
void SetDPCorr (bool _cor)
 Set to true to do the fully correlated DP systematic.
void SetOldDP (bool _old)
 Set to true to do the the old symmetric DP syst.
void SetSystDP (double size)
 Decay systematic. Shift value is 1 +/- size. 0 = off.
void SetSystNC (double size)
 NC and WS background systematic. Shift value is 1 +/- size. 0 = off.
void SetSystWS (double size)
 WS background systematic. Shift value is 1 +/- size. 0 = off.
void SetSystXSec (double size)
 All combined and NuBar cross-section systematics. Shift value is +/- size sigma. 0 = off.
void SetSystSKZP (double size)
 SKZP flux systematics. Shift value is +/- size sigma. 0 = off.
void NoSystematics ()
 Turn off all systematics.
void NoND (bool n=true)
 Don't generate a new ND spectrum each time. This means no ND systematics.
void DebugMode (bool deb=true)
 Fix systematics so the same is used in every experiment.
void PrintLoopProgress (Int_t entry, Float_t nEntries) const
 Print out a progress report in 5% increments.
int Counts (Sample_t s)
 Return the number of events in Sample s in last experiment.

Private Member Functions

void ScaleRejection (double &reject, int &nevents, double rand, double shift)
void FillMCEvents (TString mc_path)
void FillFromFilter (TTree *tin, int det)
void FillFromFCTree (TTree *tin, int det)
void FillFromDST (TString mc_path, TString xmlFileName, int det)
int Int (Sample_t s) const
void GenerateND (NuMMRunNSINu *mmRunSource)
void GenerateND (NuMMRunNSINubar *mmRunSource)
int EBin (double E) const
void Distro (Sample_t s, TString name, TString title)
TH1D * Distro (Sample_t s) const
int N (Sample_t s) const
void StoreEvent (Sample_t s, float *vals, int det)
void StoreEvent (Sample_t s, NuEvent ev, int det)
void StoreEvent (Sample_t s, NuFCEvent ev, int det)
double GetReject (Sample_t s, int i) const
void MakeSampler (Sample_t s, NuMatrixSpectrum pred)
double GetNormShift () const
double GetTrkEnergy (const NuFCEvent *ev) const
double GetShowerEnergy (const NuFCEvent *ev) const
double GetDPShift (const NuFCEvent *ev) const
double GetBGShift (const NuFCEvent *ev) const
double GetBGShift (const Sample_t s) const
double GetXSecShift (const NuFCEvent *ev) const
double GetSKZPShift (const NuFCEvent *ev) const
double GetShift (double size) const
void RollSystematics ()
void InvertShifts ()
 ClassDef (NuFCExperimentFactoryNSI, 0)

Private Attributes

bool isRHC
bool noND
const int IntND
double ndPoT
TH1D * spect [16]
TH1D * hSampler [16]
int counts [16]
std::vector< NuFCEvent * > events [16]
NuMatrixSpectrum fHist
NuMatrixSpectrum fHistND
Int_t numRecoBins
Int_t numTrueBins
std::vector< Double_t > vReco
std::vector< Double_t > vTrue
TRandom * fRandy
std::string fName
 Pointer to the random number generator.
bool fixedSyst
bool gausSyst
double norm_size
double norm_shift
double curv_size
double curv_shift
double relcurv_size
double relcurv_shift
double range_size
double range_shift
double shower_size
double shower_shift
double relshow_size
double relshow_shift
bool dp_corr
bool dp_old
double dp_size
double dp_shift
double nc_size
double nc_shift
double ws_size
double ws_shift
double xsec_size
double xsec_shift
double skzp_size
double skzp_shift

Static Private Attributes

static int counter = 0
 String name of this experiment.

Detailed Description

Factory to generate pseudo-experiments for the NSI FC analysis.

Author:
Alex Himmel, Last checkin
Author
zisvan
Version:
Revision
1.1
Date:
Date
2011/05/24 16:12:14

Created on: Sat Aug 16, 2008

Definition at line 38 of file NuFCExperimentFactoryNSI.h.


Constructor & Destructor Documentation

NuFCExperimentFactoryNSI::NuFCExperimentFactoryNSI ( TRandom3 *  rgen,
TString  nd_path,
TString  fd_path,
TString  tau_path 
)

Construct the factory.

The library files can be in one of three forms. They can be trees of NuFCEvent objects (these are available in the NuBar common area and can be generated with NuDSTAna::MakeFCTree), trees of the right variables made with NuFarNear, or simple NuDSTs.

Parameters:
rgen TRandom3* random engine.
nd_path Location of the near events library .root file.
fd_path Location of the far events library .root file.
tau_path Location of the tau events library .root file.

NuDSTAna::MakeFCTree

Definition at line 51 of file NuFCExperimentFactoryNSI.cxx.

References counter, curv_size, Distro(), dp_size, FillMCEvents(), Plot::Format(), fRandy, Msg::kInfo, Sample::kNCBackNSI, Sample::kSignalNSI, Sample::kTauNSI, Sample::kWSBackNSI, MAXMSG, nc_size, noND, norm_size, numRecoBins, numTrueBins, range_size, NuUtilities::RecoBins(), relcurv_size, relshow_size, shower_size, skzp_size, NuUtilities::TrueBins(), vReco, vTrue, ws_size, and xsec_size.

00052   : IntND(0), ndPoT(0), 
00053     fHist(), fHistND(), 
00054     fRandy(rgen),  
00055     fixedSyst(false), gausSyst(true), dp_corr(false), dp_old(false)
00056 {
00057   MAXMSG("NuFCExperimentFactoryNSI",Msg::kInfo,1) << "Creating a NuFCExperimentFactoryNSI" << endl;
00058     
00059   // Check if our random number generator is correctly setup
00060   if (!fRandy) {
00061     // If not, just use gRandom
00062     fRandy = gRandom;
00063   }
00064   noND = false;
00065   
00066   norm_size    = 0.016;
00067 
00068   curv_size    = 0.03;
00069   relcurv_size = 0.00;
00070   range_size   = 0.02;
00071   shower_size  = 1.0;
00072   relshow_size = 0.0212; // 1.85 q+ 1.05
00073   
00074   dp_size      = 0.00;
00075   //bg_size      = 0.50;
00076   nc_size      = 0.20;
00077   ws_size      = 0.30;
00078   
00079   xsec_size    = 1;
00080   skzp_size    = 0;
00081 
00082   // Get the binning scheme set up
00083   const NuUtilities cuts;
00084   NuBinningScheme::NuBinningScheme_t binningScheme =
00085   static_cast<NuBinningScheme::NuBinningScheme_t>(4);
00086   vReco = cuts.RecoBins(binningScheme);
00087   vTrue = cuts.TrueBins(binningScheme);
00088   
00089   numRecoBins = vReco.size() - 1;
00090   numTrueBins = vTrue.size() - 1;
00091   //recoBinsArray = &(vReco[0]);
00092   //trueBinsArray = &(vTrue[0]);
00093   
00094   TString suffix = TString::Format("_%i",counter);
00095 
00096   Distro(kSignalNSI, "signalSpect"+suffix, "Signal Distro");
00097   Distro(kWSBackNSI, "wsSpect"+suffix, "Wrong Sign Distro");
00098   Distro(kNCBackNSI, "ncSpect"+suffix, "NC Distro");
00099   Distro(kTauNSI, "tauSpect"+suffix, "Tau Distro");
00100 
00101   counter++;
00102   
00103   // Set up the FD MC chain
00104   FillMCEvents(nd_path);
00105   FillMCEvents(fd_path);
00106   FillMCEvents(tau_path);
00107     
00108   MAXMSG("NuFCExperimentFactoryNSI",Msg::kInfo,1) << "Done creating a NuFCExperimentFactoryNSI" << endl;
00109 }

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

Definition at line 44 of file NuFCExperimentFactoryNSI.h.

00044 {};


Member Function Documentation

NuFCExperimentFactoryNSI::ClassDef ( NuFCExperimentFactoryNSI  ,
 
) [private]
int NuFCExperimentFactoryNSI::Counts ( Sample_t  s  )  [inline]

Return the number of events in Sample s in last experiment.

Definition at line 111 of file NuFCExperimentFactoryNSI.h.

References counts, and Int().

00111 {return counts[Int(s)];};

void NuFCExperimentFactoryNSI::DebugMode ( bool  deb = true  ) 

Fix systematics so the same is used in every experiment.

Instead of randomly drawing the systematics for each new experiment, every systematic is always taken at its passed (maximum) value.

Parameters:
deb true = debug on, false = debug off. true by default.

Definition at line 992 of file NuFCExperimentFactoryNSI.cxx.

References fixedSyst, Msg::kInfo, and MSG.

00993 {
00994   fixedSyst = deb;
00995   
00996   if (fixedSyst) {
00997     MSG("NuFCExperimentFactoryNSI",Msg::kInfo) << "Now in debugging mode.  Systematics are not random -- they are fixed at the input values." << endl;
00998   }
00999   else {
01000     MSG("NuFCExperimentFactoryNSI",Msg::kInfo) << "Now out of debugging mode.  Systematics are randomly chosen based on input values." << endl;
01001   }
01002 }

TH1D * NuFCExperimentFactoryNSI::Distro ( Sample_t  s  )  const [private]

Return the unweighted spectrum of the library events for sample s

Definition at line 874 of file NuFCExperimentFactoryNSI.cxx.

References Int(), and spect.

00875 {
00876   return spect[Int(s)];
00877 }

void NuFCExperimentFactoryNSI::Distro ( Sample_t  s,
TString  name,
TString  title 
) [inline, private]

Initialize the histograms needed for sample s. This is the nweighted spectrum of events in the library and the rejection probability sampling historgram.

Definition at line 864 of file NuFCExperimentFactoryNSI.cxx.

References hSampler, Int(), numTrueBins, spect, and vTrue.

Referenced by MakeSampler(), NuFCExperimentFactoryNSI(), and StoreEvent().

00865 {
00866   spect[Int(s)] = new TH1D(name, title, numTrueBins, &(vTrue[0]));
00867   hSampler[Int(s)] = new TH1D();  // Prevents seg faults
00868 }

int NuFCExperimentFactoryNSI::EBin ( double  E  )  const [private]

A faster way of filling a histogram. Retuns the bin number that energy E refers to in binning scheme 4. Faster than Root's bin finder by more than a factor of 2.

Definition at line 786 of file NuFCExperimentFactoryNSI.cxx.

Referenced by GenerateND().

00787 {
00788   if (E < 0) return 0;
00789   if (E < 0.5) return 1;
00790   if (E < 20.) return static_cast<int>(E*4.);
00791   if (E < 30.) return static_cast<int>(E)+60;
00792   if (E < 50.) return static_cast<int>(E)/2+75;
00793   return 100;
00794 }

void NuFCExperimentFactoryNSI::FillFromDST ( TString  mc_path,
TString  xmlFileName,
int  det 
) [private]

Fill the event library for detector det from the NuDSTs in mc_path.

Definition at line 522 of file NuFCExperimentFactoryNSI.cxx.

References NuXMLConfig::AnaVersionString(), NuEvent::charge, NuInputEvents::GetEntriesNuEvent(), NuInputEvents::GetNextNuEvent(), NuEvent::iaction, NuInputEvents::InitialiseChains(), NuInputEvents::InitialiseNuEventBranch(), NuInputEvents::InputFileName(), NuEvent::inu, NuEvent::inunoosc, isRHC, Msg::kInfo, Sample::kNCBackNSI, Sample::kSignalNSI, Sample::kTauNSI, Sample::kWSBackNSI, MSG, NuInputEvents::ResetNuEventLoopPositionToStart(), and StoreEvent().

00523 {
00524   MSG("NuFCExperimentFactoryNSI",Msg::kInfo) << "Entering FillFromDST(" << mc_path << ")" << endl;
00525   
00526   // TChain holding all the FD MC
00527   NuInputEvents* mc_input = new NuInputEvents();
00528   mc_input->InputFileName(mc_path.Data());
00529   mc_input->InitialiseChains();
00530   mc_input->InitialiseNuEventBranch();  
00531   mc_input->ResetNuEventLoopPositionToStart();
00532   
00533   NuXMLConfig xmlConfig(xmlFileName);
00534   NuCutter cutter(xmlConfig.AnaVersionString(), 0);
00535   
00536   
00537   for (Int_t i=0; i < mc_input->GetEntriesNuEvent(); ++i) {
00538     NuEvent &nu=const_cast<NuEvent&>(mc_input->GetNextNuEvent(Msg::kInfo));
00539     
00540     //this->PrintLoopProgress(i,mc_input->GetEntriesNuEvent());
00541     
00542     //if (!IsGoodStdCuts(nu)) continue;
00543     cutter.MakeCuts(nu);
00544     if (cutter.Failed()) continue;
00545 
00546     if(isRHC){
00547       if (nu.charge == 1) {
00548         if (nu.inu == 16 || nu.inu == -16) {
00549           // Ignore taus from Nue's, NC's
00550           if (nu.inunoosc == 12 || nu.inunoosc==-12 || nu.iaction == 0) continue;
00551           StoreEvent(kTauNSI, nu, det);
00552         }
00553         else if (nu.iaction == 0) StoreEvent(kNCBackNSI, nu, det);
00554         else if (nu.inu == -14)   {
00555           StoreEvent(kSignalNSI, nu, det);
00556         }
00557         else StoreEvent(kWSBackNSI, nu, det);
00558       }
00559     }// end isRHC
00560 
00561     else {  
00562       if(nu.charge == -1){
00563         if (nu.inu == 16 || nu.inu == -16) {
00564           // Ignore taus from Nue's, NC's
00565           if (nu.inunoosc == 12 || nu.inunoosc==-12 || nu.iaction == 0) continue;
00566           StoreEvent(kTauNSI, nu, det);
00567         }
00568         else if (nu.iaction == 0) StoreEvent(kNCBackNSI, nu, det);
00569         else if (nu.inu == 14)   {
00570           StoreEvent(kSignalNSI, nu, det);
00571         }
00572         else StoreEvent(kWSBackNSI, nu, det);
00573       }
00574     }
00575   }
00576   MSG("NuFCExperimentFactoryNSI",Msg::kInfo) << "Leaving FillFromDST" << endl;  
00577 }

void NuFCExperimentFactoryNSI::FillFromFCTree ( TTree *  tin,
int  det 
) [private]

Fill the event library for detector det from the NuFCEvent tree tin.

Definition at line 474 of file NuFCExperimentFactoryNSI.cxx.

References NuFCEvent::charge, NuFCEvent::iaction, NuFCEvent::inu, NuFCEvent::inunoosc, isRHC, Msg::kInfo, Sample::kNCBackNSI, Sample::kSignalNSI, Sample::kTauNSI, Sample::kWSBackNSI, MSG, and StoreEvent().

Referenced by FillMCEvents().

00475 {
00476   MSG("NuFCExperimentFactoryNSI",Msg::kInfo) << "Entering FillFromFCTree" << endl;
00477   
00478   NuFCEvent *nu = new NuFCEvent();
00479   tin->SetBranchAddress("NuFCEvent", &nu);
00480 
00481   for (Int_t i=0; i < tin->GetEntries(); ++i) {
00482     tin->GetEntry(i);
00483     //this->PrintLoopProgress(i,tin->GetEntries());
00484 
00485     if (isRHC){
00486       if (nu->charge == 1) {
00487         if (nu->inu == 16 || nu->inu == -16) {
00488           // Ignore taus from Nue's, NC's
00489           if (nu->inunoosc == 12 || nu->inunoosc==-12 || nu->iaction == 0) continue;
00490           StoreEvent(kTauNSI, *nu, det);
00491         }
00492         else if (nu->iaction == 0) StoreEvent(kNCBackNSI, *nu, det);
00493         else if (nu->inu == -14)   {
00494           StoreEvent(kSignalNSI, *nu, det);
00495         }
00496         else StoreEvent(kWSBackNSI, *nu, det);
00497       }
00498     }// end isRHC
00499 
00500     else {
00501       if(nu->charge == -1){
00502         if (nu->inu == 16 || nu->inu == -16) {
00503           // Ignore taus from Nue's, NC's
00504           if (nu->inunoosc == 12 || nu->inunoosc==-12 || nu->iaction == 0) continue;
00505           StoreEvent(kTauNSI, *nu, det);
00506         }
00507         else if (nu->iaction == 0) StoreEvent(kNCBackNSI, *nu, det);
00508         else if (nu->inu == 14)   {
00509           StoreEvent(kSignalNSI, *nu, det);
00510         }
00511         else StoreEvent(kWSBackNSI, *nu, det);
00512       }
00513     }
00514   }
00515   MSG("NuFCExperimentFactoryNSI",Msg::kInfo) << "Leaving FillFromFCTree" << endl;  
00516 }

void NuFCExperimentFactoryNSI::FillFromFilter ( TTree *  tin,
int  det 
) [private]

Fill the event library for detector det from the NuFarNear Filter tree tin.

Definition at line 409 of file NuFCExperimentFactoryNSI.cxx.

References isRHC, Msg::kInfo, Sample::kNCBackNSI, Sample::kSignalNSI, Sample::kTauNSI, Sample::kWSBackNSI, MSG, and StoreEvent().

Referenced by FillMCEvents().

00410 {
00411   MSG("NuFCExperimentFactoryNSI",Msg::kInfo) << "Entering FillFromFilter" << endl;
00412   
00413   Float_t energy;       tin->SetBranchAddress("energy",       &energy);
00414   Float_t trkEn;        tin->SetBranchAddress("trkEn",        &trkEn);
00415   Float_t shwEn;        tin->SetBranchAddress("shwEn",        &shwEn);
00416   Float_t neuEnMC;      tin->SetBranchAddress("neuEnMC",      &neuEnMC);  
00417   Int_t   iaction;      tin->SetBranchAddress("iaction",      &iaction);    
00418   Int_t   inu;          tin->SetBranchAddress("inu",          &inu);
00419   Float_t charge;       tin->SetBranchAddress("charge",       &charge);
00420   Float_t ppvz;         tin->SetBranchAddress("ppvz",         &ppvz); 
00421   Int_t   containedTrk; tin->SetBranchAddress("containedTrk", &containedTrk); 
00422   Float_t fluxErr;      tin->SetBranchAddress("fluxErr",      &fluxErr); 
00423   
00424   float vals[10];
00425   for (Int_t i=0; i < tin->GetEntries(); ++i) {
00426     tin->GetEntry(i);
00427     //this->PrintLoopProgress(i,tin->GetEntries());
00428    
00429     vals[0] = energy;
00430     vals[1] = trkEn;
00431     vals[2] = shwEn;
00432     vals[3] = neuEnMC;
00433     vals[4] = (float)iaction;
00434     vals[5] = (float)inu;
00435     vals[6] = charge;
00436     vals[7] = ppvz;    
00437     vals[8] = (float)containedTrk;
00438     vals[9] = fluxErr;
00439 
00440     if(isRHC){
00441       if (charge == 1) {
00442         if (inu == 16 || inu == -16) {
00443           if (iaction == 0) continue;
00444           StoreEvent(kTauNSI, vals, det);
00445         }
00446         else if (iaction == 0) StoreEvent(kNCBackNSI, vals, det);
00447         else if (inu == -14) {
00448           StoreEvent(kSignalNSI, vals, det);
00449         }
00450         else StoreEvent(kWSBackNSI, vals, det);
00451       }
00452     }//end isRHC
00453 
00454     else {
00455       if (charge == -1)
00456         if (inu == 16 || inu == -16) {
00457           if (iaction == 0) continue;
00458           StoreEvent(kTauNSI, vals, det);
00459         }
00460         else if (iaction == 0) StoreEvent(kNCBackNSI, vals, det);
00461         else if (inu == 14)   {
00462           StoreEvent(kSignalNSI, vals, det);
00463         }
00464         else StoreEvent(kWSBackNSI, vals, det);
00465     }
00466   }
00467   MSG("NuFCExperimentFactoryNSI",Msg::kInfo) << "Leaving FillFromFilter" << endl;  
00468 }

void NuFCExperimentFactoryNSI::FillMCEvents ( TString  mc_path  )  [private]

Fill the event library for the files in mc_path. What type of files and what detector they come from is determined automatically.

Definition at line 361 of file NuFCExperimentFactoryNSI.cxx.

References det, FillFromFCTree(), FillFromFilter(), Msg::kError, Msg::kInfo, MSG, and ndPoT.

Referenced by NuFCExperimentFactoryNSI().

00362 {
00363   MSG("NuFCExperimentFactoryNSI",Msg::kInfo) << "Entering FillMCEvents(\"" << mc_path << "\")" << endl;
00364   TFile *fin = new TFile(mc_path,"READ");
00365 
00366   // Test the file opened, and fail otherwise
00367   if (fin->IsZombie()) {
00368     MSG("NuFCExperimentFactoryNSI",Msg::kError) << "Failed to read in events file " << mc_path << endl;
00369     throw runtime_error("Failed to read MC events file");
00370   }
00371   
00372   // Determine which Detector this file is for
00373   TH1D *hDetector = (TH1D*)gROOT->FindObject("hDetector");
00374   int det = (int)hDetector->GetMean();
00375  
00376   // Get the ND Pot
00377   if (det == 1) {
00378     TH1D *hTotalPot = (TH1D*)gROOT->FindObject("hTotalPot");
00379     ndPoT += hTotalPot->Integral();
00380     //    cout << "ND Pot: " << ndPoT << endl;
00381   }
00382   
00383   // Attempt to load the known tree types
00384   TTree *t1 = (TTree*)fin->Get("FDFilter");
00385   TTree *t2 = (TTree*)fin->Get("TAUFilter");
00386   TTree *t3 = (TTree*)fin->Get("NDFilter");
00387   TTree *t4 = (TTree*)fin->Get("FCTree");
00388   TTree *t5 = (TTree*)fin->Get("s");
00389 
00390   // If the tree exists in the file, load it
00391   if (t1)      FillFromFilter(t1, det);
00392   else if (t2) FillFromFilter(t2, det);
00393   else if (t3) FillFromFilter(t3, det);
00394   else if (t4) FillFromFCTree(t4, det);
00395   else if (t5) {
00396     MSG("NuFCExperimentFactoryNSI",Msg::kError) << "Use FillFromDST(path, xmlfile, detector) directly -- cannot use FillMCEvents for DST sources." << endl;
00397   }
00398   else {
00399     MSG("NuFCExperimentFactoryNSI",Msg::kError) << "No recognized tree in file: " << mc_path << endl;      
00400   }
00401   
00402   MSG("NuFCExperimentFactoryNSI",Msg::kInfo) << "Leaving FillMCEvents" << endl;  
00403 }  

void NuFCExperimentFactoryNSI::GenerateND ( NuMMRunNSINubar mmRunSource  )  [private]

Generate the ND RHC spectra.

Definition at line 720 of file NuFCExperimentFactoryNSI.cxx.

References counter, EBin(), events, fHistND, GetBGShift(), GetDPShift(), NuMMRunNSINubar::GetNDDataNSI(), GetShowerEnergy(), GetSKZPShift(), GetTrkEnergy(), GetXSecShift(), IntND, it, Msg::kInfo, MAXMSG, ndPoT, noND, numRecoBins, NuMatrix::ResetPOT(), NuMatrixSpectrum::ResetSpectrum(), and vReco.

00721 {
00722   if (noND) {
00723     MAXMSG("NuFCExperimentFactoryNSI",Msg::kInfo,1) << "Not Generating the ND" << endl;
00724     fHistND = *(mmRunSource->GetNDDataNSI());
00725     return;
00726   }
00727 
00728   else {
00729     MAXMSG("NuFCExperimentFactoryNSI",Msg::kInfo,1) << "Entering GenerateND" << endl;
00730     TString name = "ndSpect";
00731     name += counter;
00732     TH1D newND(name,"ND Spect", numRecoBins, &(vReco[0]));
00733     TH1D newDP(name+"_dp","ND Spect, DP", numRecoBins, &(vReco[0]));      
00734     /*
00735      unsigned int max = event[IntND].size();
00736      if (max == 0) {
00737      MSG("NuFCExperimentFactor",Msg::kError) << "No ND events loaded.  Cannot make an ND spectrum."
00738      }
00739      */
00740     
00741     double E;
00742     double rw;
00743     vector<NuFCEvent*>::const_iterator it;
00744     for (it = events[IntND].begin(); it != events[IntND].end(); ++it) {
00745       E = GetTrkEnergy(*it) + GetShowerEnergy(*it);
00746       rw = (*it)->rw;
00747       rw *= GetBGShift(*it);
00748       rw *= GetXSecShift(*it);
00749       rw *= GetSKZPShift(*it);
00750       
00751 //      if (dp_corr) {
00752 //        if (IsDP(*it)) newDP.AddBinContent(EBin(E), rw);
00753 //        else           newND.AddBinContent(EBin(E), rw);
00754 //      }
00755 //      else {
00756       rw *= GetDPShift(*it);
00757       newND.AddBinContent(EBin(E), rw);
00758 //      }
00759     }
00760     
00761 //    if (dp_corr) {
00762 //      // Set the decay pipe systematic size
00763 //      int lowb = 1;
00764 //      int highb = newND.GetXaxis()->FindFixBin(13.0);
00765 //      double U = newND.Integral(lowb, highb) / ndPoT * 2.92878e+20;
00766 //      double D = newDP.Integral(lowb, highb)/ ndPoT * 2.92878e+20;
00767 //      double N = 198064.;
00768 //      dp_shift = (N - U) / D;
00769 //      if (dp_shift < 0) dp_shift = 0;
00770 //      newDP.Scale(dp_shift);
00771 //      newND.Add(&newDP);
00772 //    }
00773     
00774     fHistND.ResetSpectrum(newND);
00775     fHistND.ResetPOT(ndPoT);
00776     MAXMSG("NuFCExperimentFactoryNSI",Msg::kInfo,1) << "Leaving GenerateND RHC" << endl;
00777   }
00778 }

void NuFCExperimentFactoryNSI::GenerateND ( NuMMRunNSINu mmRunSource  )  [private]

Just a copy of the standard cuts used in NuDSTana. It is needed in order to get the library events from NuDSTs. Generate the ND FHC spectra.

Definition at line 656 of file NuFCExperimentFactoryNSI.cxx.

References counter, EBin(), events, fHistND, GetBGShift(), GetDPShift(), NuMMRunNSINu::GetNDDataNSI(), GetShowerEnergy(), GetSKZPShift(), GetTrkEnergy(), GetXSecShift(), IntND, it, Msg::kInfo, MAXMSG, ndPoT, noND, numRecoBins, NuMatrix::ResetPOT(), NuMatrixSpectrum::ResetSpectrum(), and vReco.

Referenced by GenerateNewExperiment().

00657 {
00658   if (noND) {
00659     MAXMSG("NuFCExperimentFactoryNSI",Msg::kInfo,1) << "Not Generating the ND" << endl;
00660     fHistND = *(mmRunSource->GetNDDataNSI());
00661     return;
00662   }
00663  
00664   else {
00665     MAXMSG("NuFCExperimentFactoryNSI",Msg::kInfo,1) << "Entering GenerateND" << endl;
00666     
00667     TString name = "ndSpect";
00668     name += counter;
00669     TH1D newND(name,"ND Spect", numRecoBins, &(vReco[0]));
00670     TH1D newDP(name+"_dp","ND Spect, DP", numRecoBins, &(vReco[0]));      
00671     /*
00672      unsigned int max = event[IntND].size();
00673      if (max == 0) {
00674      MSG("NuFCExperimentFactor",Msg::kError) << "No ND events loaded.  Cannot make an ND spectrum."
00675      }
00676      */
00677     
00678     double E;
00679     double rw;
00680     vector<NuFCEvent*>::const_iterator it;
00681     for (it = events[IntND].begin(); it != events[IntND].end(); ++it) {
00682       E = GetTrkEnergy(*it) + GetShowerEnergy(*it);
00683       rw = (*it)->rw;
00684       rw *= GetBGShift(*it);
00685       rw *= GetXSecShift(*it);
00686       rw *= GetSKZPShift(*it);
00687       
00688 //      if (dp_corr) {
00689 //        if (IsDP(*it)) newDP.AddBinContent(EBin(E), rw);
00690 //        else           newND.AddBinContent(EBin(E), rw);
00691 //      }
00692 //      else {
00693       rw *= GetDPShift(*it);
00694       newND.AddBinContent(EBin(E), rw);
00695 //      }
00696     }
00697     
00698 //    if (dp_corr) {
00699 //      // Set the decay pipe systematic size
00700 //      int lowb = 1;
00701 //      int highb = newND.GetXaxis()->FindFixBin(13.0);
00702 //      double U = newND.Integral(lowb, highb) / ndPoT * 2.92878e+20;
00703 //      double D = newDP.Integral(lowb, highb)/ ndPoT * 2.92878e+20;
00704 //      double N = 198064.;
00705 //      dp_shift = (N - U) / D;
00706 //      if (dp_shift < 0) dp_shift = 0;
00707 //      newDP.Scale(dp_shift);
00708 //      newND.Add(&newDP);
00709 //    }
00710     
00711     fHistND.ResetSpectrum(newND);
00712     fHistND.ResetPOT(ndPoT);
00713     MAXMSG("NuFCExperimentFactoryNSI",Msg::kInfo,1) << "Leaving GenerateND FHC" << endl;
00714   }
00715 }

void NuFCExperimentFactoryNSI::GenerateNewExperiment ( NuMMRunNSINubar mmRunSource,
const NuMMParameters mmPars,
bool  reroll = true 
)

Generates a new experiment.

The internal Near, Far spectra are replaced with newly generated histograms with a new set of random systematics. The parameters are used to create the FD prediction that is the target distribution of the rejection sampling.

Parameters:
mmRunSource A NuMMRunNSINubar object
mmPars A set of parameters with the desired oscillations and/or transitions.

Definition at line 235 of file NuFCExperimentFactoryNSI.cxx.

References counter, counts, events, fHist, fRandy, GenerateND(), GetBGShift(), GetNormShift(), GetReject(), GetShowerEnergy(), GetSKZPShift(), GetTrkEnergy(), GetXSecShift(), Int(), isRHC, Msg::kInfo, Sample::kWSBackNSI, MakeSampler(), MAXMSG, MSG, NuMatrixSpectrum::Multiply(), N(), numRecoBins, NuMMRun::QuietModeOn(), NuMatrixSpectrum::ResetSpectrum(), RollSystematics(), ScaleRejection(), NuMatrixSpectrum::Spectrum(), NuMMRunNSINubar::TrueComponents(), and vReco.

00236 {
00237   MAXMSG("NuFCExperimentFactoryNSI",Msg::kInfo,1) << "Entering GenerateNewExperiment" << endl;
00238 
00239   // This is an RHC experiment
00240   isRHC = true;  
00241 
00242   TString name = "fdExp";
00243   name += (counter++);
00244 
00245   NuMatrixSpectrum fdTemp;
00246   mmRunSource->QuietModeOn();
00247 
00248   if (reroll) RollSystematics();
00249   GenerateND(mmRunSource);
00250   
00251   static int allDraws = 0;
00252   static int gthn1s = 0;
00253   static double maxProb = 0;
00254   
00255   TH1D *fdExp = new TH1D(name,"Single FC Experiment", numRecoBins, &(vReco[0]));
00256   for (int i = 11; i <= 14; i++) {
00257     Sample_t s = (Sample_t)i;
00258     
00259     fdTemp = mmRunSource->TrueComponents(mmPars, s);
00260       
00261     // Apply Normalization to the Prediction
00262     fdTemp.Multiply(GetNormShift());
00263     fdTemp.Multiply(GetBGShift(s));
00264     
00265     // Get the average events by integrating the PDF
00266     Double_t average_events = fdTemp.Spectrum()->Integral();
00267     counts[Int(s)] = 0;
00268     if (average_events == 0) continue;
00269     
00270     // Get the number of events to generate
00271     Int_t nevents = fRandy->Poisson(average_events);
00272     
00273     // Create sampling distribution
00274     MakeSampler(s, fdTemp);
00275     
00276     double numdraws = 0;
00277 
00278     // Now start sampling
00279     int ev;
00280     double rand, reject, E;
00281     bool seeking;
00282     
00283     while (nevents > 0) {
00284       // Get an event -- takes care of rejection sampling
00285       seeking = true;
00286       while (seeking) {
00287         ev = fRandy->Integer(N(s));
00288         ++numdraws;
00289         
00290         rand = fRandy->Uniform();
00291         reject = GetReject(s, ev);
00292           
00293         if (s != kWSBackNSI) {
00294             ScaleRejection(reject, nevents, rand, GetXSecShift(events[Int(s)][ev]));
00295         }
00296         ScaleRejection(reject, nevents, rand, GetSKZPShift(events[Int(s)][ev]));
00297 
00298         ++allDraws;
00299         if (reject > maxProb) maxProb = reject;
00300         if (reject > 1.) ++gthn1s;
00301         seeking = rand > reject;
00302       }
00303     
00304 
00305       // Get the new reco energy -- apply the systematic
00306       E = GetTrkEnergy(events[Int(s)][ev]) + GetShowerEnergy(events[Int(s)][ev]);
00307       
00308       // Fill our final reco spectrum
00309       fdExp->Fill(E);
00310       ++counts[Int(s)];
00311       
00312       --nevents;
00313     }
00314   }
00315   fHist.ResetSpectrum(*fdExp);
00316   delete fdExp;
00317   //fHist.ResetPOT(fdTemp.GetPOT());
00318 //   cout << "fHist before POT reset: POT, integral: " << fHist.GetPOT()
00319 //        << ", " << fHist.Integral() << endl 
00320 //        << "Setting fHist POT to: " << fdTemp.GetPOT() << endl
00321 //        << "fHist after POT reset: POT, integral: " << fHist.GetPOT()
00322 //        << ", " << fHist.Integral() << endl;
00323 
00324   float frac = gthn1s;
00325   frac /= (float)allDraws;
00326   if (frac > 0.001) {
00327     MSG("NuFCExperimentFactoryNSI",Msg::kInfo) << "So far, max wt=" << maxProb
00328                                                << ", and " << gthn1s << "/" << allDraws << " had prob > 1" << endl;
00329   }
00330   
00331   MAXMSG("NuFCExperimentFactoryNSI",Msg::kInfo,1) << "Leaving GenerateNewExperiment" << endl;
00332 }

void NuFCExperimentFactoryNSI::GenerateNewExperiment ( NuMMRunNSINu mmRunSource,
const NuMMParameters mmPars,
bool  reroll = true 
)

Generates a new experiment.

The internal Near, Far spectra are replaced with newly generated histograms with a new set of random systematics. The parameters are used to create the FD prediction that is the target distribution of the rejection sampling.

Parameters:
mmRunSource A NuMMRunNSINu object
mmPars A set of parameters with the desired oscillations and/or transitions.

Definition at line 122 of file NuFCExperimentFactoryNSI.cxx.

References counter, counts, events, fHist, fRandy, GenerateND(), GetBGShift(), GetNormShift(), GetReject(), GetShowerEnergy(), GetSKZPShift(), GetTrkEnergy(), GetXSecShift(), Int(), isRHC, Msg::kInfo, Sample::kWSBackNSI, MakeSampler(), MAXMSG, MSG, NuMatrixSpectrum::Multiply(), N(), numRecoBins, NuMMRun::QuietModeOn(), NuMatrixSpectrum::ResetSpectrum(), RollSystematics(), ScaleRejection(), NuMatrixSpectrum::Spectrum(), NuMMRunNSINu::TrueComponents(), and vReco.

00123 {
00124   MAXMSG("NuFCExperimentFactoryNSI",Msg::kInfo,1) << "Entering GenerateNewExperiment" << endl;
00125 
00126   // This is an FHC experiment
00127   isRHC = false;  
00128 
00129   TString name = "fdExp";
00130   name += (counter++);
00131 
00132   NuMatrixSpectrum fdTemp;
00133   mmRunSource->QuietModeOn();
00134 
00135   if (reroll) RollSystematics();
00136   GenerateND(mmRunSource);
00137   
00138   static int allDraws = 0;
00139   static int gthn1s = 0;
00140   static double maxProb = 0;
00141   
00142   TH1D *fdExp = new TH1D(name,"Single FC Experiment", numRecoBins, &(vReco[0]));
00143   for (int i = 11; i <= 14; i++) {
00144     Sample_t s = (Sample_t)i;
00145     
00146     //    cout << "Generating component, " << s << endl;
00147     fdTemp = mmRunSource->TrueComponents(mmPars, s);
00148     //    cout << "Got component  " << s << ", integral: " << fdTemp.Spectrum()->Integral() << endl;
00149 
00150     // Apply Normalization to the Prediction
00151     fdTemp.Multiply(GetNormShift());
00152     fdTemp.Multiply(GetBGShift(s));
00153     //    cout << "Applied shifts  " << s << ", integral: " << fdTemp.Spectrum()->Integral() << endl;
00154 
00155     // Get the average events by integrating the PDF
00156     Double_t average_events = fdTemp.Spectrum()->Integral();
00157     counts[Int(s)] = 0;
00158     if (average_events == 0) continue;
00159     
00160     // Get the number of events to generate
00161     //    cout << "Events before fluct: " << average_events << endl;
00162     Int_t nevents = fRandy->Poisson(average_events);
00163     //    cout << "Events after fluct: " << nevents << endl;
00164     // Create sampling distribution
00165     MakeSampler(s, fdTemp);
00166     
00167     double numdraws = 0;
00168 
00169     // Now start sampling
00170     int ev;
00171     double rand, reject, E;
00172     bool seeking;
00173     
00174     while (nevents > 0) {
00175       // Get an event -- takes care of rejection sampling
00176       seeking = true;
00177       while (seeking) {
00178         ev = fRandy->Integer(N(s));
00179         ++numdraws;
00180         
00181         rand = fRandy->Uniform();
00182         reject = GetReject(s, ev);
00183           
00184         if (s != kWSBackNSI) {
00185             ScaleRejection(reject, nevents, rand, GetXSecShift(events[Int(s)][ev]));
00186         }
00187         ScaleRejection(reject, nevents, rand, GetSKZPShift(events[Int(s)][ev]));
00188 
00189         ++allDraws;
00190         if (reject > maxProb) maxProb = reject;
00191         if (reject > 1.) ++gthn1s;
00192         seeking = rand > reject;
00193       }
00194     
00195 
00196       // Get the new reco energy -- apply the systematic
00197       E = GetTrkEnergy(events[Int(s)][ev]) + GetShowerEnergy(events[Int(s)][ev]);
00198       
00199       // Fill our final reco spectrum
00200 //       cout << "E= " << E 
00201 //         << ", counts= " << counts[Int(s)] 
00202 //         << ", nevents= " << nevents << endl;
00203       fdExp->Fill(E);
00204       ++counts[Int(s)];
00205       
00206       --nevents;
00207     }
00208   }
00209   fHist.ResetSpectrum(*fdExp);
00210   delete fdExp;
00211   //fHist.ResetPOT(fdTemp.GetPOT());
00212   //  cout << "Generated full spectrum fHist, integral: " << fHist.Integral() << endl; 
00213   
00214 
00215   float frac = gthn1s;
00216   frac /= (float)allDraws;
00217   if (frac > 0.001) {
00218     MSG("NuFCExperimentFactoryNSI",Msg::kInfo) << "So far, max wt=" << maxProb
00219                                                << ", and " << gthn1s << "/" << allDraws << " had prob > 1" << endl;
00220   }
00221   
00222   MAXMSG("NuFCExperimentFactoryNSI",Msg::kInfo,1) << "Leaving GenerateNewExperiment" << endl;
00223 }

double NuFCExperimentFactoryNSI::GetBGShift ( const Sample_t  s  )  const [inline, private]

Return the wrong sign/nc background systematic shift for the total number of background events if s is a background sample. Return 1 (no effect) otherwise.

Definition at line 1195 of file NuFCExperimentFactoryNSI.cxx.

References Sample::kNCBackNSI, Sample::kWSBackNSI, nc_shift, and ws_shift.

01196 {
01197   //if (s == kWrongSignPQ || s == kNCPQ) return bg_shift;
01198   if (s == kNCBackNSI) return nc_shift;
01199   if (s == kWSBackNSI) return ws_shift;
01200   return 1;
01201 }

double NuFCExperimentFactoryNSI::GetBGShift ( const NuFCEvent nu  )  const [inline, private]

Return the wrong sign/nc background systematic weight shift for this event. Return 1 (no effect) if nu is not a background event.

Definition at line 1181 of file NuFCExperimentFactoryNSI.cxx.

References NuFCEvent::charge, NuFCEvent::iaction, NuFCEvent::inu, nc_shift, and ws_shift.

Referenced by GenerateND(), and GenerateNewExperiment().

01182 {
01183   //if (nu->iaction == 0 || nu->inu == 14) return bg_shift;
01184   if (nu->iaction == 0) return nc_shift;
01185   if (nu->inu == 14 && nu->charge == 1) return ws_shift;
01186   if (nu->inu == -14 && nu->charge == -1) return ws_shift;
01187   return 1;
01188 }

double NuFCExperimentFactoryNSI::GetDPShift ( const NuFCEvent nu  )  const [inline, private]

Return the decay pipe systematic weight shift. Return 1 (no effect) if nu is not a decay pipe event.

Definition at line 1171 of file NuFCExperimentFactoryNSI.cxx.

References dp_shift, and IsDP().

Referenced by GenerateND().

01172 {
01173   if (IsDP(nu)) return dp_shift;
01174   return 1;  
01175 }

const std::string& NuFCExperimentFactoryNSI::GetName ( void   )  const [inline]

Definition at line 45 of file NuFCExperimentFactoryNSI.h.

References fName.

00045 { return fName; }

const NuMatrixSpectrum* NuFCExperimentFactoryNSI::GetNewSpectrum (  )  const [inline]

Get a new copy (pointer) to the FD Spectrum.

Definition at line 53 of file NuFCExperimentFactoryNSI.h.

References fHist.

00053 { return new NuMatrixSpectrum(fHist); }

const NuMatrixSpectrum* NuFCExperimentFactoryNSI::GetNewSpectrumND (  )  const [inline]

Get a new copy (pointer) to the ND Spectrum.

Definition at line 58 of file NuFCExperimentFactoryNSI.h.

References fHistND.

00058 { return new NuMatrixSpectrum(fHistND); }

double NuFCExperimentFactoryNSI::GetNormShift (  )  const [inline, private]

Return the normalization systematic shift. Use this to scale the total number of events being drawn.

Definition at line 1110 of file NuFCExperimentFactoryNSI.cxx.

References norm_shift.

Referenced by GenerateNewExperiment().

01111 {
01112   return norm_shift;
01113 }

double NuFCExperimentFactoryNSI::GetReject ( Sample_t  s,
int  i 
) const [private]

Get the rejection probability of event number i of sample s. Essentially, determine which energy bin this event is in and return the ratio in that bin.

MakeSampler for sample s needs to have been run for the results to be sensible.

MakeSampler

Definition at line 951 of file NuFCExperimentFactoryNSI.cxx.

References events, hSampler, Int(), and Sample::kNCBackNSI.

Referenced by GenerateNewExperiment().

00952 {
00953   double Esel, reject;
00954   if (s == kNCBackNSI) Esel = events[Int(s)][i]->energy;
00955   else                 Esel = events[Int(s)][i]->neuEnMC;
00956   
00957   int bin = hSampler[Int(s)]->GetXaxis()->FindFixBin(Esel);
00958   reject = hSampler[Int(s)]->GetBinContent(bin);
00959   return reject;
00960 }

double NuFCExperimentFactoryNSI::GetShift ( double  size  )  const [inline, private]

Get a systematic shift of the size passed. This can be random or or fixed depending on the debugging mode used.

Definition at line 1041 of file NuFCExperimentFactoryNSI.cxx.

References fixedSyst, fRandy, and gausSyst.

Referenced by RollSystematics().

01042 { 
01043   if (fixedSyst)     return size;
01044   else if (gausSyst) return fRandy->Gaus(0, size);
01045   else               return fRandy->Uniform(-size, size); 
01046 }

double NuFCExperimentFactoryNSI::GetShowerEnergy ( const NuFCEvent nu  )  const [inline, private]

Return the systematically shifted shower energy. Accounts for both overall and relative shower energy scales.

Definition at line 1134 of file NuFCExperimentFactoryNSI.cxx.

References NuFCEvent::detector, relshow_shift, shower_shift, and NuFCEvent::shwEn.

Referenced by GenerateND(), and GenerateNewExperiment().

01135 {
01136   // Energy-dependent systematic as described in doc-7173
01137   double offset = 6.6;
01138   double scale = 3.5;
01139   double eLife = 1.44; //GeV
01140   double enShift = 1.0 + shower_shift*0.01*
01141                    (offset + scale*TMath::Exp(-1.0*nu->shwEn/eLife));
01142 
01143   double E = nu->shwEn * enShift;
01144   if (nu->detector == 2) E *= relshow_shift;
01145 //  cout << "GetShowerEnergy:" 
01146 //  << " shwEn=" << nu->shwEn 
01147 //  << ", shower_shift=" << shower_shift
01148 //  << ", shower_size=" << shower_size
01149 //  << ", enShift=" << enShift 
01150 //  << ", relshow_shift=" << relshow_shift
01151 //  << ", returning E=" << E << endl;
01152   return E;
01153 }

double NuFCExperimentFactoryNSI::GetSKZPShift ( const NuFCEvent nu  )  const [inline, private]

Return the skzp systematic shift for this event. It is based on the fluxErr stored in nu.

Definition at line 1231 of file NuFCExperimentFactoryNSI.cxx.

References NuFCEvent::fluxErr, and skzp_shift.

Referenced by GenerateND(), and GenerateNewExperiment().

01232 {
01233   return 1. + skzp_shift*(nu->fluxErr - 1.);
01234 }

const NuMatrixSpectrum& NuFCExperimentFactoryNSI::GetSpectrum (  )  const [inline]

Get a reference to the FD Spectrum.

Definition at line 51 of file NuFCExperimentFactoryNSI.h.

References fHist.

00051 { return fHist; }

const NuMatrixSpectrum& NuFCExperimentFactoryNSI::GetSpectrumND (  )  const [inline]

Get a reference to the ND Spectrum.

Definition at line 56 of file NuFCExperimentFactoryNSI.h.

References fHistND.

00056 { return fHistND; }

double NuFCExperimentFactoryNSI::GetTrkEnergy ( const NuFCEvent nu  )  const [inline, private]

Return the systematically shifted track energy. Accounts for both range and curvature as required for this event.

Definition at line 1119 of file NuFCExperimentFactoryNSI.cxx.

References NuFCEvent::containedTrk, curv_shift, NuFCEvent::detector, range_shift, relcurv_shift, and NuFCEvent::trkEn.

Referenced by GenerateND(), and GenerateNewExperiment().

01120 {
01121   double E = nu->trkEn;
01122   if (nu->containedTrk == 1) E *= range_shift;
01123   else {
01124     E *= curv_shift;
01125     if (nu->detector == 2) E *= relcurv_shift;
01126   }
01127   return E;
01128 }

double NuFCExperimentFactoryNSI::GetXSecShift ( const NuFCEvent ev  )  const [inline, private]

Return the cross-section systematic weight shift for CC events. This includes both combined cross-section errors and the NuBar-Nu systematic ratio errors.

Definition at line 1207 of file NuFCExperimentFactoryNSI.cxx.

References MuELoss::a, NuFCEvent::charge, NuFCEvent::iaction, NuFCEvent::neuEnMC, and xsec_shift.

Referenced by GenerateND(), and GenerateNewExperiment().

01208 {
01209   if (ev->iaction != 1) return 1;
01210   if (ev->charge == -1) return 1; // X Sec is specific to NuBars
01211   
01212   static double xp = 25.;
01213   static double a = 0.895;
01214   static double xc = 4.0;
01215   static double b = 7.59e-3;
01216   static double c = -8.05e-4;
01217   double E = ev->neuEnMC;
01218   double shift = -0.0617;
01219   
01220   if (E < 25) shift += a - 1 + 2.*(1.-a)*E/xp + (a - 1.)*E*E/(xp*xp);
01221   if (E < xc) shift += (c*(xc-E)*(xc-E)*(xc-E)+b*(xc-E)*(xc-E));
01222   
01223   shift = 1 + xsec_shift * shift;
01224   return shift;
01225 }

int NuFCExperimentFactoryNSI::Int ( Sample_t  s  )  const [private]

Return the integer value corresponding to sample s. Used to access the appropriate array elements.

Definition at line 822 of file NuFCExperimentFactoryNSI.cxx.

References Sample::kAppearedNQ, Sample::kAppearedPQ, Msg::kError, Sample::kNCBackNSI, Sample::kNCNQ, Sample::kNCPQ, Sample::kSignalNQ, Sample::kSignalNSI, Sample::kSignalPQ, Sample::kTauNQ, Sample::kTauNSI, Sample::kTauPQ, Sample::kWrongSignNQ, Sample::kWrongSignPQ, Sample::kWSBackNSI, and MSG.

Referenced by Counts(), Distro(), GenerateNewExperiment(), GetReject(), MakeSampler(), N(), and StoreEvent().

00823 {
00824   switch (s) {
00825   case kSignalPQ:
00826     return 1;
00827   case kWrongSignPQ:
00828     return 2;
00829   case kNCPQ:
00830     return 3;
00831   case kTauPQ:
00832     return 4;
00833   case kAppearedPQ:
00834     return 5;
00835   case kSignalNQ:
00836     return 6;
00837   case kWrongSignNQ:
00838     return 7;
00839   case kNCNQ:
00840     return 8;
00841   case kTauNQ:
00842     return 9;
00843   case kAppearedNQ:
00844     return 10;
00845   case kSignalNSI:
00846     return 11;
00847   case kWSBackNSI:
00848     return 12;
00849   case kNCBackNSI:
00850     return 13;
00851   case kTauNSI:
00852     return 14;
00853   default:
00854     MSG("NuFCExperimentFactor",Msg::kError) << s << " is not a known Sample_t.  Returning -1" << endl;
00855     return -1;
00856   }
00857 }

void NuFCExperimentFactoryNSI::InvertShifts (  )  [private]

Set every shift (the randomly drawn ones) to its negative. No longer in use.

Definition at line 1086 of file NuFCExperimentFactoryNSI.cxx.

References curv_shift, dp_shift, nc_shift, norm_shift, range_shift, relcurv_shift, relshow_shift, shower_shift, skzp_shift, ws_shift, and xsec_shift.

01087 {
01088   // Get a single copy of these systematics
01089   norm_shift = 2 - norm_shift;
01090   
01091   curv_shift = 2 - curv_shift;
01092   relcurv_shift = 2 - curv_shift;
01093   range_shift = 2 - range_shift;
01094   shower_shift = 2 - shower_shift;
01095   relshow_shift = 2 - relshow_shift;
01096   
01097   dp_shift = 2 - dp_shift;
01098   //bg_shift = 2 - bg_shift;
01099   nc_shift = 2 - nc_shift;
01100   ws_shift = 2 - ws_shift;
01101   
01102   xsec_shift *= -1;
01103   skzp_shift *= -1;
01104 }

bool NuFCExperimentFactoryNSI::IsDP ( const NuFCEvent nu  )  const

Return true if nu is a decay pipe event.

Definition at line 1159 of file NuFCExperimentFactoryNSI.cxx.

References NuFCEvent::ppvz, and NuFCEvent::ptype.

Referenced by GetDPShift().

01160 {
01161   if (nu->ppvz < 4500) return false;
01162   if (abs(nu->ptype) == 13) return false;
01163   return true;
01164   
01165 }

void NuFCExperimentFactoryNSI::MakeSampler ( Sample_t  s,
NuMatrixSpectrum  pred 
) [private]

Create a sampling distribution (the ratio of target distribution to the distribution we currently have) for rejection sampling. This ratio will be normalized so its maximum value is 1.

Parameters:
s The sample were creating the sampler for.
pred The target distribution.

Definition at line 971 of file NuFCExperimentFactoryNSI.cxx.

References Distro(), hSampler, Int(), and NuMatrixSpectrum::Spectrum().

Referenced by GenerateNewExperiment().

00972 {
00973   if (hSampler[Int(s)]) delete hSampler[Int(s)];
00974   hSampler[Int(s)] = new TH1D(*(pred.Spectrum()));
00975   hSampler[Int(s)]->Divide(Distro(s));
00976   //  cout << "hSampler[Int(s)]->GetMaximum(): " << hSampler[Int(s)]->GetMaximum() << endl;
00977   if(hSampler[Int(s)]->GetMaximum() != 0)
00978     hSampler[Int(s)]->Scale(0.87/hSampler[Int(s)]->GetMaximum());
00979   // Scale to 0.87 (empirical #) instead of 1.0 to avoid 
00980   // clipping problems from having syst. shifted prob > 1
00981 }

int NuFCExperimentFactoryNSI::N ( Sample_t  s  )  const [inline, private]

Definition at line 147 of file NuFCExperimentFactoryNSI.h.

References events, and Int().

Referenced by GenerateNewExperiment().

00147 {return events[Int(s)].size();};

void NuFCExperimentFactoryNSI::NoND ( bool  n = true  )  [inline]

Don't generate a new ND spectrum each time. This means no ND systematics.

Definition at line 105 of file NuFCExperimentFactoryNSI.h.

References n, and noND.

00105 {noND = n;};

void NuFCExperimentFactoryNSI::NoSystematics (  ) 

Turn off all systematics.

Turns off all systematics and resets all the random shifts to 0. To test just one systematic, turn off all of them and then set the one you want to the desired value with SetSyst*.

Definition at line 1012 of file NuFCExperimentFactoryNSI.cxx.

References curv_size, dp_size, Msg::kInfo, MAXMSG, nc_size, norm_size, range_size, relcurv_size, relshow_size, RollSystematics(), shower_size, skzp_size, ws_size, and xsec_size.

01013 {
01014   MAXMSG("NuFCExperimentFactoryNSI",Msg::kInfo,1) << "Turning off systematics" << endl;
01015   norm_size = 0;
01016   
01017   curv_size = 0;
01018   relcurv_size = 0;
01019   range_size = 0;
01020   shower_size = 0;
01021   relshow_size = 0;
01022   
01023   dp_size = 0;
01024   //bg_size = 0;
01025   nc_size = 0;
01026   ws_size = 0;
01027   
01028   xsec_size = 0;
01029   skzp_size = 0;
01030   
01031   // Zero out the shifts
01032   RollSystematics();
01033   MAXMSG("NuFCExperimentFactoryNSI",Msg::kInfo,1) << "Done turning off systematics" << endl;
01034 }

void NuFCExperimentFactoryNSI::PrintLoopProgress ( Int_t  entry,
Float_t  nEntries 
) const

Print out a progress report in 5% increments.

Parameters:
entry The current count.
nEntries The total count.

Definition at line 804 of file NuFCExperimentFactoryNSI.cxx.

References Msg::kInfo, and MSG.

00805 {
00806   Float_t fract=ceil(nEntries/20.);  
00807   if (ceil(((Float_t)entry)/fract)==((Float_t)entry)/fract){
00808     MSG("NuFCExperimentFactoryNSI",Msg::kInfo) 
00809     <<"Fraction of loop complete: "<<entry 
00810     <<"/"<<nEntries<<"  ("
00811     <<(Int_t)(100.*entry/nEntries)<<"%)"<<endl;
00812   }
00813 }

void NuFCExperimentFactoryNSI::RollSystematics (  )  [private]

Generate a new set of random systematic shifts.

Definition at line 1052 of file NuFCExperimentFactoryNSI.cxx.

References curv_shift, curv_size, dp_corr, dp_old, dp_shift, dp_size, GetShift(), nc_shift, nc_size, norm_shift, norm_size, range_shift, range_size, relcurv_shift, relcurv_size, relshow_shift, relshow_size, shower_shift, shower_size, skzp_shift, skzp_size, ws_shift, ws_size, xsec_shift, and xsec_size.

Referenced by GenerateNewExperiment(), and NoSystematics().

01053 {
01054   // Get a single copy of these systematics
01055   norm_shift = 1 + GetShift(norm_size);
01056   
01057   relcurv_shift = 1 + GetShift(relcurv_size);
01058 
01059   double tmp = GetShift(curv_size);
01060   double ratio = 0;
01061   if (curv_size != 0) ratio = range_size/curv_size;
01062   curv_shift = 1 + tmp;
01063   range_shift = 1 + tmp*ratio;
01064   
01065   shower_shift = GetShift(shower_size);
01066   relshow_shift = 1 + GetShift(relshow_size);
01067 
01068   if (!dp_corr) {
01069     if (dp_size == 0) dp_shift = 1;
01070     else if (dp_old)  dp_shift = 1.0 + GetShift(dp_size);
01071     else              dp_shift = 0.8 + GetShift(dp_size);
01072     if (dp_shift < 0) dp_shift = 0;
01073   }
01074   //bg_shift = 1 + GetShift(bg_size);
01075   nc_shift = 1 + GetShift(nc_size);
01076   ws_shift = 1 + GetShift(ws_size);
01077   
01078   xsec_shift = GetShift(xsec_size);
01079   skzp_shift = GetShift(skzp_size);
01080 }

void NuFCExperimentFactoryNSI::ScaleRejection ( double &  reject,
int &  nevents,
double  rand,
double  shift 
) [inline, private]

This is the rejection-sampling equivalent to reweighting an events. A new, shifted rejection probability is calculated -- this will be used to determine whether to keep or reject the event. But, if this new rejection probability will keep an event that would have been rejected, increase the total number we draw to account for this change in normalization. If we now reject an event we would have kept, decrease the total number of events.

Parameters:
reject The original rejection probability. This may be modified to the new rejection probability.
nevents The number of events remaining to be drawn for this experiment. This may be modified.
rand The random number drawn for rejection purposes.
shift The shift in event "weight."

Definition at line 347 of file NuFCExperimentFactoryNSI.cxx.

Referenced by GenerateNewExperiment().

00348 { 
00349   double reject_new = reject*shift;
00350   if (rand > reject && rand < reject_new) nevents++;
00351   if (rand < reject && rand > reject_new) nevents--;
00352   
00353   reject = reject_new;
00354 }

void NuFCExperimentFactoryNSI::SetDPCorr ( bool  _cor  )  [inline]

Set to true to do the fully correlated DP systematic.

Definition at line 85 of file NuFCExperimentFactoryNSI.h.

References dp_corr.

00085 {dp_corr = _cor;};

void NuFCExperimentFactoryNSI::SetGausSyst ( bool  _gs  )  [inline]

Set to true to have gaussian distributed syst values and false to have them be uniform.

Definition at line 61 of file NuFCExperimentFactoryNSI.h.

References gausSyst.

00061 {gausSyst = _gs;}

void NuFCExperimentFactoryNSI::SetOldDP ( bool  _old  )  [inline]

Set to true to do the the old symmetric DP syst.

Definition at line 87 of file NuFCExperimentFactoryNSI.h.

References dp_old.

00087 {dp_old = _old;};

void NuFCExperimentFactoryNSI::SetSystCurv ( double  size  )  [inline]

Absolute Track Energy from curvature far systematic. Shift value is 1 +/- size. 0 = off.

Definition at line 72 of file NuFCExperimentFactoryNSI.h.

References curv_size.

00072 {curv_size = size;};

void NuFCExperimentFactoryNSI::SetSystDP ( double  size  )  [inline]

Decay systematic. Shift value is 1 +/- size. 0 = off.

Definition at line 89 of file NuFCExperimentFactoryNSI.h.

References dp_size.

00089 {dp_size = size;};

void NuFCExperimentFactoryNSI::SetSystNC ( double  size  )  [inline]

NC and WS background systematic. Shift value is 1 +/- size. 0 = off.

NC background systematic. Shift value is 1 +/- size. 0 = off.

Definition at line 93 of file NuFCExperimentFactoryNSI.h.

References nc_size.

00093 {nc_size = size;};

void NuFCExperimentFactoryNSI::SetSystNorm ( double  size  )  [inline]

Normalization systematic. Shift value is 1 +/- size. 0 = off.

SetSyst* Functions -- set the size of these systematic shifts The actual shifted value is varied from experiment to experiment Shift values are defined as 1 +/- size Set to 0 to turn off that systematic

Definition at line 69 of file NuFCExperimentFactoryNSI.h.

References norm_size.

00069 {norm_size = size;};

void NuFCExperimentFactoryNSI::SetSystRange ( double  size  )  [inline]

Track Energy from range systematic. Shift value is 1 +/- size. 0 = off.

Definition at line 76 of file NuFCExperimentFactoryNSI.h.

References range_size.

00076 {range_size = size;};

void NuFCExperimentFactoryNSI::SetSystRelCurv ( double  size  )  [inline]

Relative Track Energy from curvature far systematic. Shift value is 1 +/- size. 0 = off.

Definition at line 74 of file NuFCExperimentFactoryNSI.h.

References relcurv_size.

00074 {relcurv_size = size;};

void NuFCExperimentFactoryNSI::SetSystRelShower ( double  size  )  [inline]

Relative shower energy (far only). Shift value is 1 +/- size. 0 = off.

Definition at line 80 of file NuFCExperimentFactoryNSI.h.

References relshow_size.

00080 {relshow_size = size;};

void NuFCExperimentFactoryNSI::SetSystShower ( double  size  )  [inline]

Shower Energy overall systematic. Shift value is +/- size sigma. 0 = off.

Definition at line 78 of file NuFCExperimentFactoryNSI.h.

References shower_size.

00078 {shower_size = size;};

void NuFCExperimentFactoryNSI::SetSystSKZP ( double  size  )  [inline]

SKZP flux systematics. Shift value is +/- size sigma. 0 = off.

Definition at line 100 of file NuFCExperimentFactoryNSI.h.

References skzp_size.

00100 {skzp_size = size;};

void NuFCExperimentFactoryNSI::SetSystWS ( double  size  )  [inline]

WS background systematic. Shift value is 1 +/- size. 0 = off.

Definition at line 95 of file NuFCExperimentFactoryNSI.h.

References ws_size.

00095 {ws_size = size;};

void NuFCExperimentFactoryNSI::SetSystXSec ( double  size  )  [inline]

All combined and NuBar cross-section systematics. Shift value is +/- size sigma. 0 = off.

Definition at line 98 of file NuFCExperimentFactoryNSI.h.

References xsec_size.

00098 {xsec_size = size;};

void NuFCExperimentFactoryNSI::StoreEvent ( Sample_t  s,
NuFCEvent  ev,
int  det 
) [private]

Store the NuFCEvent ev in the sample s cache. Set its detector type to det.

Definition at line 918 of file NuFCExperimentFactoryNSI.cxx.

References NuFCEvent::detector, Distro(), NuFCEvent::energy, events, Int(), IntND, Sample::kNCBackNSI, Msg::kWarning, MSG, and NuFCEvent::neuEnMC.

00919 {
00920   NuFCEvent * loc = new NuFCEvent(ev);
00921   loc->detector = det;
00922   
00923   if (loc->detector == 2) { // Far Det and Taus
00924     events[Int(s)].push_back(loc);
00925     
00926     double E;
00927     if (s == kNCBackNSI) E = loc->energy;
00928     else                 E = loc->neuEnMC;
00929   
00930   Distro(s)->Fill(E);
00931   }
00932   else if (loc->detector == 1) {
00933     events[IntND].push_back(loc);
00934   }
00935   else {
00936     MSG("NuFCExperimentFactoryNSI",Msg::kWarning) << "Detector " << det << " not recognized.  Event not stored." << endl;
00937   }
00938 }

void NuFCExperimentFactoryNSI::StoreEvent ( Sample_t  s,
NuEvent  ev,
int  det 
) [private]

Store the NuEvent ev in the sample s cache. Set its detector type to det. It will be converted to a NuFCEvent first.

Definition at line 907 of file NuFCExperimentFactoryNSI.cxx.

References StoreEvent().

00908 {
00909   NuFCEvent loc(ev);
00910   StoreEvent(s, loc, det);
00911 }

void NuFCExperimentFactoryNSI::StoreEvent ( Sample_t  s,
float *  vals,
int  det 
) [private]

Store the event described by the values in vals in the sample s cache. Set its detector typ to det. The expected order of fields is: energy, trkEn, shwEn, neuEnMC, iaction, inu, charge, ppvz, containedTrk, fluxErr.

Definition at line 885 of file NuFCExperimentFactoryNSI.cxx.

References NuFCEvent::charge, NuFCEvent::containedTrk, NuFCEvent::energy, NuFCEvent::fluxErr, NuFCEvent::iaction, NuFCEvent::inu, NuFCEvent::neuEnMC, NuFCEvent::ppvz, NuFCEvent::shwEn, and NuFCEvent::trkEn.

Referenced by FillFromDST(), FillFromFCTree(), FillFromFilter(), and StoreEvent().

00886 {
00887   NuFCEvent loc;
00888   loc.energy       = vals[0];
00889   loc.trkEn        = vals[1];
00890   loc.shwEn        = vals[2];
00891   loc.neuEnMC      = vals[3];
00892   loc.iaction      = (Int_t)vals[4];
00893   loc.inu          = (Int_t)vals[5];
00894   loc.charge       = (Int_t)vals[6];
00895   loc.ppvz         = vals[7];
00896   loc.containedTrk = (Int_t)vals[8];
00897   loc.fluxErr      = vals[9];    
00898   
00899   StoreEvent(s, loc, det);
00900 }


Member Data Documentation

int NuFCExperimentFactoryNSI::counter = 0 [static, private]

String name of this experiment.

Definition at line 188 of file NuFCExperimentFactoryNSI.h.

Referenced by GenerateND(), GenerateNewExperiment(), and NuFCExperimentFactoryNSI().

Definition at line 169 of file NuFCExperimentFactoryNSI.h.

Referenced by Counts(), and GenerateNewExperiment().

Definition at line 201 of file NuFCExperimentFactoryNSI.h.

Referenced by GetTrkEnergy(), InvertShifts(), and RollSystematics().

Definition at line 207 of file NuFCExperimentFactoryNSI.h.

Referenced by RollSystematics(), and SetDPCorr().

Definition at line 208 of file NuFCExperimentFactoryNSI.h.

Referenced by RollSystematics(), and SetOldDP().

Definition at line 209 of file NuFCExperimentFactoryNSI.h.

Referenced by GetDPShift(), InvertShifts(), and RollSystematics().

std::vector<NuFCEvent*> NuFCExperimentFactoryNSI::events[16] [private]

Definition at line 175 of file NuFCExperimentFactoryNSI.h.

Referenced by GenerateNewExperiment(), GetNewSpectrum(), and GetSpectrum().

Definition at line 176 of file NuFCExperimentFactoryNSI.h.

Referenced by GenerateND(), GetNewSpectrumND(), and GetSpectrumND().

Definition at line 196 of file NuFCExperimentFactoryNSI.h.

Referenced by DebugMode(), and GetShift().

std::string NuFCExperimentFactoryNSI::fName [private]

Pointer to the random number generator.

Definition at line 186 of file NuFCExperimentFactoryNSI.h.

Referenced by GetName().

Definition at line 197 of file NuFCExperimentFactoryNSI.h.

Referenced by GetShift(), and SetGausSyst().

Definition at line 168 of file NuFCExperimentFactoryNSI.h.

Referenced by Distro(), GetReject(), and MakeSampler().

const int NuFCExperimentFactoryNSI::IntND [private]

Definition at line 136 of file NuFCExperimentFactoryNSI.h.

Referenced by GenerateND(), and StoreEvent().

Definition at line 211 of file NuFCExperimentFactoryNSI.h.

Referenced by GetBGShift(), InvertShifts(), and RollSystematics().

Definition at line 137 of file NuFCExperimentFactoryNSI.h.

Referenced by FillMCEvents(), and GenerateND().

Definition at line 135 of file NuFCExperimentFactoryNSI.h.

Referenced by GenerateND(), NoND(), and NuFCExperimentFactoryNSI().

Definition at line 199 of file NuFCExperimentFactoryNSI.h.

Referenced by GetNormShift(), InvertShifts(), and RollSystematics().

Definition at line 180 of file NuFCExperimentFactoryNSI.h.

Referenced by Distro(), and NuFCExperimentFactoryNSI().

Definition at line 203 of file NuFCExperimentFactoryNSI.h.

Referenced by GetTrkEnergy(), InvertShifts(), and RollSystematics().

Definition at line 202 of file NuFCExperimentFactoryNSI.h.

Referenced by GetTrkEnergy(), InvertShifts(), and RollSystematics().

Definition at line 205 of file NuFCExperimentFactoryNSI.h.

Referenced by GetShowerEnergy(), InvertShifts(), and RollSystematics().

Definition at line 204 of file NuFCExperimentFactoryNSI.h.

Referenced by GetShowerEnergy(), InvertShifts(), and RollSystematics().

Definition at line 215 of file NuFCExperimentFactoryNSI.h.

Referenced by GetSKZPShift(), InvertShifts(), and RollSystematics().

TH1D* NuFCExperimentFactoryNSI::spect[16] [private]

Definition at line 167 of file NuFCExperimentFactoryNSI.h.

Referenced by Distro().

std::vector<Double_t> NuFCExperimentFactoryNSI::vReco [private]
std::vector<Double_t> NuFCExperimentFactoryNSI::vTrue [private]

Definition at line 182 of file NuFCExperimentFactoryNSI.h.

Referenced by Distro(), and NuFCExperimentFactoryNSI().

Definition at line 212 of file NuFCExperimentFactoryNSI.h.

Referenced by GetBGShift(), InvertShifts(), and RollSystematics().

Definition at line 214 of file NuFCExperimentFactoryNSI.h.

Referenced by GetXSecShift(), InvertShifts(), and RollSystematics().


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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1