NCExtrapolation Class Reference

Base class for the various extrapolations More...

#include <NCExtrapolation.h>

Inheritance diagram for NCExtrapolation:
NCExtrapolationBeamMatrix NCExtrapolationFarNear NCExtrapolationNone NCExtrapolationPID

List of all members.

Classes

struct  InterpCacheCompare
 Helper for fInterpCache. More...

Public Types

typedef std::map< NCBeam::Info,
double > 
POTMap
typedef std::map< NCBeam::Info,
double > 
ScaleMap

Public Member Functions

 NCExtrapolation ()
virtual ~NCExtrapolation ()
virtual void AddEvent (NCEventInfo eventInfo, bool useMCAsData, NCType::EFileType fileType, NCBeam::Info beamInfo)
virtual void DoneFilling ()
 Called after the final AddEvent() but before any calls to FindSpectraForPars().
virtual void Prepare (const Registry &r)
virtual void WriteResources (const NC::OscProb::OscPars *trueOscPars)
void DrawBestFitSpectra (Detector::Detector_t det)
void DrawBestFitRatios (Detector::Detector_t det)
virtual TString GetShortName () const =0
 This is the name used to name things in the output file etc.
virtual TString GetLongName () const =0
 This is the name the extrapolation is known under on plots and such.
virtual void Reset (bool nearData, bool farData, bool nearMC, bool farMC)
 Causes the extrapolation to forget everything it knows about the requested spectra.
NC::Fitter::CoordNDim GetBestFitCoordNDim () const
const NC::CoordinateConverterGetCoordConv () const
NC::OscProb::OscParsGetBestFitOscPars () const
 Get best fit oscillation parameters found by the fitter or set manually.
double FindChiSqrForPars (const NC::OscProb::OscPars *oscPars, const NC::SystPars &systPars, double *analyticNorm=0)
double ChiSqrAtBestFit () const
 The best $ \chi^2 $ value found by the fitter.
void SetBestFitCoordNDim (const NC::Fitter::CoordNDim minCoord)
 Set fMinCoord and prepare GetBestFitSysts and GetBestFitOscPars.
virtual bool EnableNearToFarExtrapolation (bool enable)
 Enable or disable corrections to FD spectra based on ND spectra.
virtual std::vector< const TH1 * > GetNearMCSpectra (std::vector< NCBeam::Info > beamsToUse)
void InitializeInterpolator (const NC::OscProb::OscPars *osc)
 Create fInterpolator and initialize it with shifted spectra.

Static Public Member Functions

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

Protected Types

typedef std::map< std::pair
< NCBeam::Info,
Detector::Detector_t >, NCBeam * > 
fBeams_t

Protected Member Functions

virtual void FindSpectraForPars (const NC::OscProb::OscPars *oscPars, const NC::SystPars &systPars, std::vector< NCBeam::Info > beamsToUse, std::vector< TH1 * > &exps, std::vector< TH1 * > &obss)=0
 Override this in the derived class.
virtual void CleanupSpectra (std::vector< TH1 * > exps, std::vector< TH1 * > obss)
 Called after FindSpectraForPars() to delete necessary spectra.
NC::SystPars GetBestFitSysts () const
 Get best fit systematics found by the fitter or set manually.
void SetBestFitSysts (NC::SystPars systs)
 Warning - doesn't affect fMinCoord.
void SetBestFitOscPars (NC::OscProb::OscPars *pars)
 Warning - doesn't affect fMinCoord.
NCBeamGetBeam (Detector::Detector_t det, NCBeam::Info beamInfo, bool add=false)
void FillResultHistograms (NCBeam *beam)
double LogLikelihood (const std::vector< TH1 * > &exp, const std::vector< TH1 * > &obs) const
 The log-likelihood formula from the PDG.
double CalculateBestNormalization (double M, double D) const
 Analytically calculate the best normalization.
std::vector< double > SystsAsSigmas (const NCBeam::Info &info, bool *simpleScale=0) const
 Convert details of systematic shifts into a list of sigmas.

Protected Attributes

fBeams_t fBeams
bool fUseCC
bool fUseNC
bool fHideFDData
 Preserve blinding by not showing FD data in plots.
NCType::EOscModel fOscillationModel
bool fSoftPenalty
double fEnergyThreshold
 Don't accept any events below this energy.
NC::CoordinateConverter fCoordConv
bool fDoSystematicInterpolation
 TODO - remove, is in FitMaster.
NC::ISpectrumInterpolatorfInterpolator
bool fUseSystPenalty
bool fAddEventsToEnergyBins
 Enables us to disable filling expensive vectors in NCEnergyBin.

Private Types

typedef std::map< const
NC::OscProb::OscPars
*, NC::ISpectrumInterpolator
*, InterpCacheCompare
InterpCache_t

Private Attributes

double fMinChiSqr
bool fPrintFDEvents
 Whether to print all FD data events.
NC::Fitter::CoordNDim fMinCoord
 This is the best fit point as found by the fitter.
NC::SystPars fBestFitSysts
NC::OscProb::OscParsfBestFitOscPars
InterpCache_t fInterpCache

Detailed Description

Base class for the various extrapolations

Definition at line 40 of file NCExtrapolation.h.


Member Typedef Documentation

typedef std::map<std::pair<NCBeam::Info, Detector::Detector_t>, NCBeam*> NCExtrapolation::fBeams_t [protected]

vector of objects describing each beam, ie LE-10 185 ND, LE-10 185 FD, pME ND, pHE ND, etc

Definition at line 190 of file NCExtrapolation.h.

Definition at line 289 of file NCExtrapolation.h.

typedef std::map<NCBeam::Info, double> NCExtrapolation::POTMap

Definition at line 51 of file NCExtrapolation.h.

typedef std::map<NCBeam::Info, double> NCExtrapolation::ScaleMap

Definition at line 52 of file NCExtrapolation.h.


Constructor & Destructor Documentation

NCExtrapolation::NCExtrapolation (  )  [inline]

Definition at line 43 of file NCExtrapolation.h.

00044   :  fEnergyThreshold(0),
00045      fInterpolator(0),
00046      fUseSystPenalty(true),
00047      fBestFitOscPars(0)
00048   {
00049   }

NCExtrapolation::~NCExtrapolation (  )  [virtual]

Definition at line 46 of file NCExtrapolation.cxx.

References fBeams, fBestFitOscPars, fInterpCache, and it.

00047 {
00048   if(fBestFitOscPars) delete fBestFitOscPars;
00049 
00050   for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it)
00051     delete it->second;
00052 
00053   // Includes fInterpolator somewhere
00054   for(InterpCache_t::iterator it = fInterpCache.begin();
00055       it != fInterpCache.end(); ++it){
00056     delete it->first;
00057     delete it->second;
00058   } // end for it
00059 }


Member Function Documentation

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

Reimplemented in NCExtrapolationNone, and NCExtrapolationPID.

Definition at line 121 of file NCExtrapolation.cxx.

References NCBeam::AddEvent(), NCEventInfo::analysis, MsgStream::AttachOStream(), NCEventInfo::beam, ANtpHeaderInfo::dataType, ANtpHeaderInfo::detector, fEnergyThreshold, Plot::Format(), fPrintFDEvents, GetBeam(), MsgService::GetStream(), NCEventInfo::header, ANtpRecoInfo::inFiducialVolume, MsgService::Instance(), ANtpAnalysisInfo::isCC, ANtpAnalysisInfo::isNC, ANtpRecoInfo::isSimpleCutsClean, SimFlag::kData, Detector::kFar, Msg::kInfo, Detector::kNear, ANtpHeaderInfo::month, MSG, ANtpRecoInfo::nuEnergy, NCEventInfo::reco, run(), ANtpHeaderInfo::run, ANtpHeaderInfo::snarl, ANtpHeaderInfo::subRun, NCEventInfo::truth, and ANtpHeaderInfo::year.

00125 {
00126   //method to add events to near detector beams.  all methods should see
00127   //the same near detector spectra, so there is no loss of generalization.
00128 
00129   //dont bother with events outside the fiducial volume
00130   if(info.reco->inFiducialVolume < 1 ||
00131      (info.header->detector == int(Detector::kNear) &&
00132       info.reco->isSimpleCutsClean < 1)
00133      )
00134      return;
00135 
00136 
00137   if(info.reco->nuEnergy < fEnergyThreshold) return;
00138 
00139 
00140   //need to check if this type of beam already exists in the vector.
00141   //if so add the event to that beam. if not, make it the beam then
00142   //the add event to it
00143   Detector::Detector_t detector = Detector::Detector_t(info.header->detector);
00144 
00145   NCBeam* beam = GetBeam(detector, beamInfo, true);
00146 
00147   beam->AddEvent(info.header, info.beam, info.reco,
00148                  info.analysis, info.truth, useMCAsData, fileType);
00149 
00150   // TODO - I don't know what this was achieving, but whatever it was it isn't anymore
00151   /*
00152   //if this is MC and not useMCAsData, look to see if there are
00153   //other run types for the same beam.  default run type for MC is
00154   //NCType::kRunI
00155   if(info.header->dataType == int(SimFlag::kMC) && !useMCAsData){
00156     for(int i = NCType::kRunII; i < NCType::kMaxRun+1; ++i){
00157       beam = FindBeamIndex(detector, beamType, NCType::ERunType(i));
00158       if(beam < 100)
00159         fBeams[beam]->AddEvent(info.header, info.beam, info.reco,
00160                               info.analysis, info.truth, useMCAsData, fileType);
00161     }//end loop over possible runs
00162   }//end if mc and needs to be added to other runs for this beam
00163   */
00164 
00165   //print out the far detector events if requested
00166   if(fPrintFDEvents && detector == Detector::kFar &&
00167      info.header->dataType == int(SimFlag::kData)){
00168 
00169     static bool once=true;
00170     if(once){
00171       once=false;
00172       MsgService::Instance()->GetStream("DataEventListAll")->AttachOStream(Msg::kInfo,"ncutils_events_all.txt");
00173 
00174       MsgService::Instance()->GetStream("DataEventListCC")->AttachOStream(Msg::kInfo,"ncutils_events_cc.txt");
00175 
00176       MsgService::Instance()->GetStream("DataEventListNC")->AttachOStream(Msg::kInfo,"ncutils_events_nc.txt");
00177     }
00178     TString stream = "";
00179     TString date = TString::Format("%d-%02d ", info.header->year, info.header->month);
00180     TString run = TString::Format("F000%d", info.header->run);
00181 
00182     if(info.analysis->isNC < 1 && info.analysis->isCC > 0) stream = "DataEventListCC";
00183     if(info.analysis->isNC > 0 && info.analysis->isCC < 1) stream = "DataEventListNC";
00184 
00185     if(stream!="")
00186       MSG(stream, Msg::kInfo)
00187         << date
00188         << run << " " << info.header->subRun
00189         << " " << info.header->snarl
00190         << " " << info.reco->nuEnergy
00191         << " " << info.analysis->isCC
00192         << " " << info.analysis->isNC << endl;
00193 
00194       MSG("DataEventListAll", Msg::kInfo)
00195         << date
00196         << run << " " << info.header->subRun
00197         << " " << info.header->snarl
00198         << " " << info.reco->nuEnergy
00199         << " " << info.analysis->isCC
00200         << " " << info.analysis->isNC << endl;
00201 
00202   }
00203 }

double NCExtrapolation::CalculateBestNormalization ( double  M,
double  D 
) const [protected]

Analytically calculate the best normalization.

Parameters:
M Total number of events in oscillated, shifted Monte Carlo
D Total number of events in data
Returns:
The normalization s that minimizes $ \Delta\chi^2 $ for the log likelihood formula plus quadratic penalty for s

See docdb-4628/normcalcwriteup.pdf for derivation

See also:
LogLikelihood

Definition at line 634 of file NCExtrapolation.cxx.

References fUseSystPenalty, NCType::kNormalization, NCType::kParams, and SQR.

Referenced by FindChiSqrForPars().

00635 {
00636   // sigma -> infinity
00637   if(!fUseSystPenalty) return D/M;
00638 
00639   const double sigmaSqr = SQR(NCType::kParams[kNormalization].sigma);
00640 
00641   const double x = (1-sigmaSqr*M)/2;
00642 
00643   return x + sqrt(SQR(x)+D*sigmaSqr);
00644 }

double NCExtrapolation::ChiSqrAtBestFit (  )  const [inline]

The best $ \chi^2 $ value found by the fitter.

If this function is called before a fit has been performed the result will be garbage

Definition at line 128 of file NCExtrapolation.h.

References fMinChiSqr.

00128 {return fMinChiSqr;}

void NCExtrapolation::CleanupSpectra ( std::vector< TH1 * >  exps,
std::vector< TH1 * >  obss 
) [protected, virtual]

Called after FindSpectraForPars() to delete necessary spectra.

Parameters:
exps The expected spectra just returned from FindSpectraForPars()
obss The observed spectra just returned from FindSpectraForPars()

Reimplemented in NCExtrapolationBeamMatrix, NCExtrapolationFarNear, NCExtrapolationNone, and NCExtrapolationPID.

Definition at line 591 of file NCExtrapolation.cxx.

Referenced by FindChiSqrForPars(), and InitializeInterpolator().

00593 {
00594   // Do nothing. This is for derived classes to, optionally, override
00595 }

const Registry & NCExtrapolation::DefaultConfig ( void   )  [static]

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

Reimplemented in NCExtrapolationBeamMatrix, NCExtrapolationFarNear, NCExtrapolationNone, and NCExtrapolationPID.

Definition at line 62 of file NCExtrapolation.cxx.

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

Referenced by NCExtrapolationModule::DefaultConfig().

00063 {
00064   static Registry r;
00065 
00066   r.UnLockValues();
00067 
00068   r.Set("EnergyThreshold",           0);
00069 
00070   r.Set("UseCCEvents",               true);
00071   r.Set("UseNCEvents",               true);
00072   r.Set("Prediction",                false);
00073   r.Set("PredictionWriteData",       false);
00074 
00075   r.Set("PrintFDEvents",             false);
00076   r.Set("HideFDData",                false);
00077 
00078   r.Set("DoSystematicInterpolation", false);
00079   r.Set("DisableSystPenalty",        false);
00080   r.Set("AddEventsToEnergyBins",     true);
00081 
00082   r.LockValues();
00083   return r;
00084 }

virtual void NCExtrapolation::DoneFilling (  )  [inline, virtual]

Called after the final AddEvent() but before any calls to FindSpectraForPars().

Reimplemented in NCExtrapolationBeamMatrix, and NCExtrapolationPID.

Definition at line 63 of file NCExtrapolation.h.

00063 {}

void NCExtrapolation::DrawBestFitRatios ( Detector::Detector_t  det  ) 

!! mc_hist->Add(beam->GetMCFitNuEToNuEHistogram(i));

Definition at line 319 of file NCExtrapolation.cxx.

References bfld::AsString(), fBeams, fHideFDData, NCBeam::GetDataGraph(), NCBeam::GetDetector(), NCBeam::GetInfo(), NCBeam::GetMCFitHistogram(), NCBeam::GetMCHistogram(), GetShortName(), it, NCType::kCC, Detector::kFar, NCType::kNC, n, and one().

Referenced by WriteResources().

00320 {
00321   for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
00322     NCBeam* beam = it->second;
00323     if(beam->GetDetector() != det) continue;
00324 
00325     TString name = "bestFitRatiosCanv"+GetShortName();
00326     name += Detector::AsString(det);
00327     name += beam->GetInfo().GetDescription();
00328     TCanvas *canv = new TCanvas(name, name, 150, 10, 1200, 300);
00329     TPad *pad1 = new TPad(name+"pad1", "", 0.00, 0.0, 0.50, 1.0, 10);
00330     TPad *pad2 = new TPad(name+"pad2", "", 0.50, 0.0, 1.00, 1.0, 10);
00331     pad1->Draw();
00332     pad2->Draw();
00333     pad1->cd();
00334 
00335     for(NCType::EEventType i = NCType::kNC;
00336         i <= NCType::kCC;
00337         i = NCType::EEventType(int(i)+1)){
00338       if(i == NCType::kCC) pad2->cd();
00339 
00340       TString title = Detector::AsString(det);
00341       if(i == NCType::kNC)
00342         title += " Detector Neutral Current Selected";
00343       else
00344         title += " Detector Charged Current Selected";
00345 
00346       TH1* mc_ratio = (TH1*)beam->GetMCFitHistogram(i)->Clone();
00347       TH1* mc_hist = (TH1*)beam->GetMCHistogram(i)->Clone();
00349       mc_ratio->Divide(mc_hist);
00350 
00351       mc_ratio->SetTitle(title);
00352       mc_ratio->GetYaxis()->SetRangeUser(0, 1.5);
00353       mc_ratio->GetYaxis()->SetTitle("Ratio to no oscillations");
00354       mc_ratio->Draw("][");
00355 
00356       // This option plots only the ND data when making a prediction.
00357       // This effectively preserves blinding while being able to make
00358       // a FD prediction using the actual data.
00359 
00360       TGraphAsymmErrors* data_ratio = 0;
00361 
00362       if(!fHideFDData || beam->GetDetector() != Detector::kFar){
00363         data_ratio = (TGraphAsymmErrors*)beam->GetDataGraph(i)->Clone();
00364 
00365         for(int n = 0; n < data_ratio->GetN(); ++n){
00366           double x, y;
00367           data_ratio->GetPoint(n, x, y);
00368           double hi = y+data_ratio->GetErrorYhigh(n);
00369           double lo = y-data_ratio->GetErrorYlow(n);
00370           double div = mc_hist->GetBinContent(mc_hist->FindBin(x));
00371           if(div == 0){ // Don't want to see this point...
00372             div = 1;
00373             x = -1;
00374           }
00375           y /= div; hi /= div; lo /= div;
00376           data_ratio->SetPoint(n, x, y);
00377           data_ratio->SetPointEYhigh(n, hi-y);
00378           data_ratio->SetPointEYlow(n, y-lo);
00379         }
00380 
00381         data_ratio->Draw("pe1same");
00382       }
00383 
00384       TLine* one = new TLine(0, 1, mc_ratio->GetXaxis()->GetXmax(), 1);
00385       one->SetLineStyle(7); // medium dashed
00386       one->Draw("same");
00387 
00388 
00389       TLegend* leg = new TLegend(0.45, 0.35, 0.85, 0.15);
00390       leg->SetBorderSize(0);
00391       if(data_ratio) leg->AddEntry(data_ratio, "Data", "l");
00392       leg->AddEntry(mc_ratio, "Best Fit", "l");
00393 
00394       leg->Draw();
00395 
00396     }//end loop over NC/CC
00397 
00398     canv->Update();
00399 
00400     // Append the canvas object to the current directory
00401     // otherwise it can't be found and saved later
00402     gDirectory->Append(canv);
00403 
00404   }//end loop over beams
00405 }

void NCExtrapolation::DrawBestFitSpectra ( Detector::Detector_t  det  ) 

Loop over all the beams. Draw pairs of plots for each beam in the chosen detector - NC and CC

Definition at line 207 of file NCExtrapolation.cxx.

References bfld::AsString(), fBeams, fHideFDData, NCBeam::GetDataGraph(), NCBeam::GetDetector(), NCBeam::GetInfo(), NCBeam::GetMCBackgroundHistogram(), NCBeam::GetMCFitHistogram(), NCBeam::GetMCFitNuEToNuEHistogram(), NCBeam::GetMCFitNuMuToNuEHistogram(), NCBeam::GetMCFitNuMuToNuTauHistogram(), NCBeam::GetMCHistogram(), NCBeam::GetMCNoFitBackgroundHistogram(), NCBeam::GetMCNoFitNuEToNuEHistogram(), GetShortName(), it, NCType::kCC, Detector::kFar, and NCType::kNC.

Referenced by WriteResources().

00208 {
00209   for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
00210     NCBeam* beam = it->second;
00211     if(beam->GetDetector() != det) continue;
00212 
00213     TString name = "bestFitSpectraCanv"+GetShortName();
00214     name += Detector::AsString(det);
00215     name += beam->GetInfo().GetDescription();
00216     TCanvas *canv = new TCanvas(name, name, 150, 10, 1200, 300);
00217     TPad *pad1 = new TPad(name+"pad1", "", 0.00, 0.0, 0.50, 1.0, 10);
00218     TPad *pad2 = new TPad(name+"pad2", "", 0.50, 0.0, 1.00, 1.0, 10);
00219     pad1->Draw();
00220     pad2->Draw();
00221     pad1->cd();
00222 
00223     /*
00224     TString bestDeltaM;
00225     bestDeltaM.Format("#Deltam^{2} = %5.4f", BestDeltaMSqr());
00226     TString bestUMu3Sqr;
00227     bestUMu3Sqr.Format("sin^{2}(2#theta) = %3.2f", BestUMu3Sqr());
00228     TString bestUS3Sqr;
00229     bestUS3Sqr.Format("f_{s} = %3.2f", BestUS3Sqr());
00230     TString minChiSqr;
00231     minChiSqr.Format("#chi^{2} = %3.2f", fMinChiSqr);
00232     */
00233 
00234     for(NCType::EEventType i = NCType::kNC;
00235         i <= NCType::kCC;
00236         i = NCType::EEventType(int(i)+1)){
00237       if(i == NCType::kCC) pad2->cd();
00238 
00239       double maxY = 1.3*TMath::Max(beam->GetMCHistogram(i)->GetMaximum(),
00240                                    beam->GetDataGraph(i)->GetHistogram()->GetMaximum());
00241 
00242       beam->GetMCHistogram(i)->SetMaximum(maxY);
00243       TString title = Detector::AsString(det);
00244       if(i == NCType::kNC)
00245         title += " Detector Neutral Current Selected";
00246       else
00247         title += " Detector Charged Current Selected";
00248       beam->GetMCHistogram(i)->SetTitle(title);
00249       beam->GetMCHistogram(i)->Draw();
00250       //if the fit histogram was filled, draw the fits
00251       if(beam->GetMCFitHistogram(i)->Integral() > 0.){
00252         beam->GetMCFitHistogram(i)->Draw("same");
00253         beam->GetMCBackgroundHistogram(i)->Draw("same");
00254         beam->GetMCFitNuMuToNuTauHistogram(i)->Draw("same");
00255         beam->GetMCFitNuEToNuEHistogram(i)->Draw("same");
00256         beam->GetMCFitNuMuToNuEHistogram(i)->Draw("same");
00257       }
00258       else{
00259         //no fit, no point drawing NuMuToNuTau or NuMuToNuE
00260         beam->GetMCNoFitBackgroundHistogram(i)->Draw("same");
00261         beam->GetMCNoFitNuEToNuEHistogram(i)->Draw("same");
00262       }
00263 
00264       // Don't show the annoying 20-120 GeV bin by default. It's still
00265       // there in case anyone wants it
00266       beam->GetMCHistogram(i)->GetXaxis()->SetRangeUser(0, 19.9);
00267 
00268       // This option plots only the ND data when making a prediction.
00269       // This effectively preserves blinding while being able to make
00270       // a FD prediction using the actual data.
00271 
00272       if(!fHideFDData || beam->GetDetector() != Detector::kFar){
00273         beam->GetDataGraph(i)->Draw("pe1same");
00274       }
00275 
00276       TLegend *leg = new TLegend(0.45, 0.65, 0.85, 0.85);
00277       leg->SetBorderSize(0);
00278       leg->AddEntry(beam->GetDataGraph(i), "Data", "lep");
00279       leg->AddEntry(beam->GetMCHistogram(i),
00280                     "Monte Carlo (Nominal)", "l");
00281       if(beam->GetMCFitHistogram(i)->Integral() > 0.){
00282         leg->AddEntry(beam->GetMCFitHistogram(i),
00283                       "Monte Carlo (Best Fit)", "l");
00284         leg->AddEntry(beam->GetMCBackgroundHistogram(i),
00285                       "Background", "bf");
00286         leg->AddEntry(beam->GetMCFitNuMuToNuTauHistogram(i),
00287                       "Background - #nu_{#tau} CC", "bf");
00288         leg->AddEntry(beam->GetMCFitNuEToNuEHistogram(i),
00289                       "Background - Beam #nu_{e} CC", "bf");
00290         leg->AddEntry(beam->GetMCFitNuMuToNuEHistogram(i),
00291                       "Background - #nu_{e} CC", "bf");
00292       }
00293       else{
00294         leg->AddEntry(beam->GetMCNoFitBackgroundHistogram(i),
00295                       "Background", "bf");
00296         leg->AddEntry(beam->GetMCNoFitNuEToNuEHistogram(i),
00297                       "Background - Beam #nu_{e} CC", "bf");
00298       }
00299 //       leg->AddEntry(beam->GetDataGraph(i), bestDeltaM, "");
00300 //       leg->AddEntry(beam->GetDataGraph(i), bestUS3Sqr, "");
00301 //       leg->AddEntry(beam->GetDataGraph(i), bestUMu3Sqr, "");
00302 //       leg->AddEntry(beam->GetDataGraph(i), minChiSqr, "");
00303       leg->Draw();
00304 
00305     }//end loop over NC/CC
00306 
00307 
00308     canv->Update();
00309 
00310     // Append the canvas object to the current directory
00311     // otherwise it can't be found and saved later
00312     gDirectory->Append(canv);
00313 
00314   }//end loop over beams
00315 }

bool NCExtrapolation::EnableNearToFarExtrapolation ( bool  enable  )  [virtual]

Enable or disable corrections to FD spectra based on ND spectra.

Parameters:
enable true to perform extrapolation. false to give produce unextrapolated spectra
Returns:
The previous value of this setting

The default value of this setting should be to true

Reimplemented in NCExtrapolationBeamMatrix, NCExtrapolationFarNear, and NCExtrapolationPID.

Definition at line 692 of file NCExtrapolation.cxx.

Referenced by InitializeInterpolator().

00693 {
00694   assert(0 && "EnableNearToFarExtrapolation needs to be overridden");
00695 }

void NCExtrapolation::FillResultHistograms ( NCBeam beam  )  [protected]

Definition at line 502 of file NCExtrapolation.cxx.

References NCBeam::FillResultHistograms(), GetBestFitOscPars(), and GetBestFitSysts().

Referenced by WriteResources().

00503 {
00504   NC::OscProb::OscPars* oscPars = GetBestFitOscPars();
00505   if(!oscPars) return;
00506   const NC::SystPars systPars = GetBestFitSysts();
00507   beam->FillResultHistograms(oscPars, &systPars);
00508 }

double NCExtrapolation::FindChiSqrForPars ( const NC::OscProb::OscPars oscPars,
const NC::SystPars systPars,
double *  analyticNorm = 0 
)
Parameters:
analyticNorm If nonzero then flags that the best normalization should be calculated analytically and applied and the normalization factor returned in this variable
See also:
CalculateBestNormalization

Definition at line 511 of file NCExtrapolation.cxx.

References CalculateBestNormalization(), CleanupSpectra(), fBeams, FindSpectraForPars(), fUseSystPenalty, NCBeam::GetDetector(), NCBeam::GetInfo(), NCBeam::GetRunType(), it, Detector::kFar, NC::RunUtil::kRunAll, LogLikelihood(), and NC::SystPars::PenaltyForShifts().

Referenced by NCExtrapolationModule::DoMockExperiments(), and NC::GetChiSqrFromDerived::EvalAtEx().

00514 {
00515   // Because profiling reveals we spend all our time in ~TH1 removing ourselves
00516   // from the current directory. We restore this at the end of the function,
00517   // need to remember not to take any early returns
00518   const bool addDirectoryRestore = TH1::AddDirectoryStatus();
00519   TH1::AddDirectory(false);
00520 
00521 
00522   vector<TH1*> exps, obss;
00523 
00524   // For the moment just use all the beams
00525   std::vector<NCBeam::Info> beamsToUse;
00526   for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
00527     NCBeam* beam = it->second;
00528     // Only want to ask for each physical beam once (arbitrarily FD)
00529     // Don't take the overall combined "kRunAll" beam
00530     if(beam->GetDetector() == Detector::kFar &&
00531        beam->GetRunType() != NC::RunUtil::kRunAll &&
00532        beam->GetInfo().IsNominal())
00533       beamsToUse.push_back(beam->GetInfo());
00534   }
00535 
00536   FindSpectraForPars(oscPars, systPars, beamsToUse, exps, obss);
00537 
00538 
00539   if(analyticNorm){
00540     assert(exps.size() == obss.size());
00541     assert(exps.size() == 1); // TODO, this shouldn't be necessary
00542     const double M = exps[0]->Integral();
00543     const double D = obss[0]->Integral();
00544     //    cout << M << " " << D;
00545     *analyticNorm = CalculateBestNormalization(M, D);
00546     //    cout << " -> " << *analyticNorm << endl;
00547     exps[0]->Scale(*analyticNorm);
00548   }
00549 
00550   const double ret = LogLikelihood(exps, obss)+
00551     (fUseSystPenalty ? systPars.PenaltyForShifts() : 0);
00552 
00553   CleanupSpectra(exps, obss);
00554 
00555 
00556   TH1::AddDirectory(addDirectoryRestore);
00557 
00558   return ret;
00559 }

virtual void NCExtrapolation::FindSpectraForPars ( const NC::OscProb::OscPars oscPars,
const NC::SystPars systPars,
std::vector< NCBeam::Info beamsToUse,
std::vector< TH1 * > &  exps,
std::vector< TH1 * > &  obss 
) [protected, pure virtual]

Override this in the derived class.

Parameters:
oscPars Oscillation parameters
systPars Systematic shifts
beamsToUse Only fill exps and obss with spectra from beams with these indices
[out] exps To be filled with expected spectra
[out] obss To be filled with the corresponding observed spectra

Implemented in NCExtrapolationBeamMatrix, NCExtrapolationFarNear, NCExtrapolationNone, and NCExtrapolationPID.

Referenced by FindChiSqrForPars(), and InitializeInterpolator().

NCBeam * NCExtrapolation::GetBeam ( Detector::Detector_t  det,
NCBeam::Info  beamInfo,
bool  add = false 
) [protected]

Definition at line 475 of file NCExtrapolation.cxx.

References bfld::AsString(), fAddEventsToEnergyBins, fBeams, fUseCC, fUseNC, NCBeam::Info::GetDescription(), it, Msg::kInfo, and MSG.

Referenced by AddEvent(), and NCExtrapolationFarNear::WriteResources().

00479 {
00480   pair<NCBeam::Info, Detector::Detector_t> key(beamInfo, det);
00481   fBeams_t::iterator it = fBeams.find(key);
00482   if(it != fBeams.end()) return it->second;
00483 
00484   if(!add) return 0;
00485 
00486   NCBeam* beam = new NCBeam(det, beamInfo, fUseCC, fUseNC, fAddEventsToEnergyBins);
00487   fBeams[key] = beam;
00488 
00489   MSG("NCExtrapolation", Msg::kInfo) << "adding new beam "
00490                                      << beamInfo.GetDescription()
00491                                      << " to "
00492                                      << Detector::AsString(det)
00493                                      << endl;
00494 
00495   return beam;
00496 }

NC::Fitter::CoordNDim NCExtrapolation::GetBestFitCoordNDim (  )  const [inline]

Definition at line 97 of file NCExtrapolation.h.

References fMinCoord.

00097 {return fMinCoord;}

NC::OscProb::OscPars* NCExtrapolation::GetBestFitOscPars (  )  const [inline]

Get best fit oscillation parameters found by the fitter or set manually.

Once MakeDeltaChiSqrMaps has found the minimum, or after SetBestFitCoordNDim or SetBestFitOscPars have been called returns the oscillation parameters component of the best fit point.

Before then it will return a null pointer.

Definition at line 110 of file NCExtrapolation.h.

References fCoordConv, fMinCoord, and NC::CoordinateConverter::OscParsFromCoordNDim().

Referenced by FillResultHistograms(), NCExtrapolationBeamMatrix::WriteResources(), NCExtrapolationPID::WriteResources(), and NCExtrapolationFarNear::WriteResources().

00111   {
00112     return fCoordConv.OscParsFromCoordNDim(fMinCoord);
00113   }

NC::SystPars NCExtrapolation::GetBestFitSysts (  )  const [inline, protected]

Get best fit systematics found by the fitter or set manually.

Once MakeDeltaChiSqrMaps has found the minimum, or after SetBestFitCoordNDim or SetBestFitSysts have been called returns the systematic component of the best fit point.

Before then it will return unshifted systematics.

Definition at line 179 of file NCExtrapolation.h.

References fCoordConv, fMinCoord, and NC::CoordinateConverter::SystParsFromCoordNDim().

Referenced by FillResultHistograms(), NCExtrapolationBeamMatrix::WriteResources(), NCExtrapolationPID::WriteResources(), and NCExtrapolationFarNear::WriteResources().

const NC::CoordinateConverter* NCExtrapolation::GetCoordConv (  )  const [inline]

Definition at line 99 of file NCExtrapolation.h.

References fCoordConv.

Referenced by NC::GetChiSqrFromDerived::EvalAtEx().

00099 {return &fCoordConv;}

virtual TString NCExtrapolation::GetLongName (  )  const [pure virtual]

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

Implemented in NCExtrapolationBeamMatrix, NCExtrapolationFarNear, NCExtrapolationNone, and NCExtrapolationPID.

virtual std::vector<const TH1*> NCExtrapolation::GetNearMCSpectra ( std::vector< NCBeam::Info beamsToUse  )  [virtual]
virtual TString NCExtrapolation::GetShortName (  )  const [pure virtual]

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

Implemented in NCExtrapolationBeamMatrix, NCExtrapolationFarNear, NCExtrapolationNone, and NCExtrapolationPID.

Referenced by DrawBestFitRatios(), DrawBestFitSpectra(), NC::FitMaster::WriteResources(), and WriteResources().

void NCExtrapolation::InitializeInterpolator ( const NC::OscProb::OscPars osc  ) 

Create fInterpolator and initialize it with shifted spectra.

Definition at line 707 of file NCExtrapolation.cxx.

References NC::ISpectrumInterpolator::AddInputSpectra(), CleanupSpectra(), NC::Utility::CloneFast(), det, NC::Utility::DivideFast(), EnableNearToFarExtrapolation(), fBeams, fDoSystematicInterpolation, FindSpectraForPars(), fInterpCache, fInterpolator, NCBeam::Info::GetDescription(), NC::ISpectrumInterpolator::GetInputSpectra(), GetNearMCSpectra(), NCBeam::Info::GetRunType(), info, it, Detector::kFar, NC::RunUtil::kRunAll, n, and SystsAsSigmas().

00708 {
00709   assert(fDoSystematicInterpolation);
00710 
00711   using namespace NC::Utility;
00712 
00713   // Ignore calls that come as a result of running this function
00714   static bool already = false;
00715   if(already) return;
00716   already = true;
00717 
00718   // If we already have an interpolator cached for these oscillation
00719   // parameters then use just that.
00720   if(fInterpCache.find(osc) != fInterpCache.end()){
00721     fInterpolator = fInterpCache[osc];
00722     already = false;
00723     return;
00724   }
00725 
00726   // Stop the cache getting too huge (interpolators hold lots of histograms)
00727   if(fInterpCache.size() > 100){
00728     const NC::OscProb::OscPars* oldKey = 0;
00729     for(InterpCache_t::iterator it = fInterpCache.begin();
00730         it != fInterpCache.end(); ++it){
00731       if(fInterpolator && it->second == fInterpolator) oldKey = it->first;
00732       if(it->second != fInterpolator){
00733         delete it->first;
00734         delete it->second;
00735       }
00736     } // end for it
00737     fInterpCache.clear();
00738     if(oldKey) fInterpCache[oldKey] = fInterpolator;
00739   }
00740 
00741   // We don't want any of the FindSpectraForPars to use the interpolator
00742   NC::ISpectrumInterpolator* oldInterp = fInterpolator;
00743   fInterpolator = 0;
00744 
00745   // Disable extrapolation so that we get independently shifted spectrum in FD
00746   const bool farNearRestore = EnableNearToFarExtrapolation(false);
00747 
00748   // Build a list of beam infos describing each separate running condition
00749   vector<NCBeam::Info> allBeams;
00750   for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
00751     const pair<NCBeam::Info, Detector::Detector_t> key = it->first;
00752     const NCBeam::Info info = key.first;
00753     const Detector::Detector_t det = key.second;
00754 
00755     // RunAll isn't a "real" beam
00756     if(info.GetRunType() == NC::RunUtil::kRunAll) continue;
00757     // Arbitrary - just need to get each beam description once
00758     if(det != Detector::kFar) continue;
00759 
00760     allBeams.push_back(info);
00761   }
00762 
00763   // beamsToUseNom = [b for b in allBeams if b.IsNominal()]
00764   vector<NCBeam::Info> beamsToUseNom;
00765   for(unsigned int n = 0; n < allBeams.size(); ++n){
00766     if(allBeams[n].IsNominal()) beamsToUseNom.push_back(allBeams[n]);
00767   }
00768 
00769   vector<TH1*> junk;    // We don't need to get data spectra
00770   vector<TH1*> noms_fd; // Nominal FD spectra
00771   FindSpectraForPars(osc, NC::SystPars(), beamsToUseNom, noms_fd, junk);
00772 
00773   vector<const TH1*> noms_nd = GetNearMCSpectra(beamsToUseNom);
00774 
00775   // The index order of the spectra is:
00776   // (Run I FD) (Run I ND) (Run II FD) (Run II ND) etc.
00777   vector<const TH1*> noms;
00778   assert(noms_nd.size() == noms_fd.size());
00779   noms.reserve(2*noms_fd.size());
00780   for(unsigned int n = 0; n < noms_fd.size(); ++n){
00781     noms.push_back(CloneFast(noms_fd[n]));
00782     noms.push_back(CloneFast(noms_nd[n]));
00783   }
00784 
00785   CleanupSpectra(noms_fd, junk);
00786 
00787   // Don't assign to fInterpolator yet, as we still have some more calls
00788   // to FindSpectraForPars to make that don't require interpolation.
00789   NC::ISpectrumInterpolator* interp = new NC::SpectrumInterpolatorSimplest;
00790 
00791   // If an entry in this map is filled, then this shift has already been
00792   // dealt with for all runs by the simpleScale optimization below.
00793   map<vector<double>, bool> alreadyFilled;
00794 
00795   // We accumulate expected spectra binned by what shifts they have.
00796   typedef map<vector<double>, vector<const TH1*> > ExpsMap;
00797   ExpsMap expsMap;
00798 
00799   for(unsigned int beamIdx = 0; beamIdx < allBeams.size(); ++beamIdx){
00800     const NCBeam::Info info = allBeams[beamIdx];
00801 
00802     bool simpleScale;
00803     const vector<double> shift = SystsAsSigmas(info, &simpleScale);
00804 
00805     // The simpleScale optimization already fired and there's no need to add
00806     // any more ratios for this shift setting
00807     if(alreadyFilled[shift]) continue;
00808 
00809     if(simpleScale && oldInterp){
00810       // These are the ratios for this shift for _all_ beams
00811       vector<TH1*> ratios = oldInterp->GetInputSpectra(shift);
00812       // So set the flag so that we don't do this again and double-count
00813       alreadyFilled[shift] = true;
00814       interp->AddInputSpectra(shift, ratios);
00815       continue;
00816     } // end if
00817 
00818     vector<NCBeam::Info> beamsToUse(1, info);
00819     vector<TH1*> exps_fd;
00820     FindSpectraForPars(osc, NC::SystPars(), beamsToUse, exps_fd, junk);
00821 
00822     vector<const TH1*> exps_nd = GetNearMCSpectra(beamsToUse);
00823 
00824     vector<const TH1*> exps; // shifted expectations
00825     assert(exps_fd.size() == exps_nd.size());
00826     exps.reserve(2*exps_fd.size());
00827     // The spectra go in the vector matching noms, indexed like
00828     // (Run I FD) (Run I ND) (Run II FD) (Run II ND)
00829     for(unsigned int n = 0; n < exps_fd.size(); ++n){
00830       TH1* exp_fd = CloneFast(exps_fd[n]);
00831       exp_fd->SetName((info.GetDescription()+" FD").Data());
00832       exps.push_back(exp_fd);
00833 
00834       TH1* exp_nd = CloneFast(exps_nd[n]);
00835       exp_nd->SetName((info.GetDescription()+" ND").Data());
00836       exps.push_back(exp_nd);
00837     }
00838 
00839     CleanupSpectra(exps_fd, junk);
00840 
00841     // expsMap[shift].append(exps)
00842     vector<const TH1*>& already = expsMap[shift];
00843     already.insert(already.end(), exps.begin(), exps.end());
00844   } // end for beamIdx
00845 
00846   for(ExpsMap::iterator it = expsMap.begin(); it != expsMap.end(); ++it){
00847     vector<double> shift = it->first;
00848     vector<const TH1*> exps = it->second;
00849     // The shifted spectra are stored in the interpolator as ratios over the
00850     // nominal spectra.
00851     vector<TH1*> ratios;
00852     ratios.reserve(exps.size());
00853     assert(exps.size() == noms.size());
00854     for(unsigned int n = 0; n < exps.size(); ++n){
00855       ratios.push_back(CloneFast(exps[n]));
00856       delete exps[n];
00857       exps[n] = 0;
00858       DivideFast(ratios[n], noms[n]);
00859     }
00860     interp->AddInputSpectra(shift, ratios);
00861   } // end for it
00862 
00863   // All of these were Clone'd and need deleting
00864   for(unsigned int n = 0; n < noms.size(); ++n) delete noms[n];
00865 
00866   // Put extrapolation setting back to however it was.
00867   EnableNearToFarExtrapolation(farNearRestore);
00868 
00869   // Have to wait till the last moment, because if fSpectrumInterpolator exists
00870   // then FindSpectraForPars tries to use it -> chicken and egg problem.
00871   fInterpolator = interp;
00872 
00873   assert(fInterpCache.find(osc) == fInterpCache.end());
00874   fInterpCache[const_cast<NC::OscProb::OscPars*>(osc)->Clone()] = fInterpolator;
00875 
00876   already = false;
00877 }

double NCExtrapolation::LogLikelihood ( const std::vector< TH1 * > &  exp,
const std::vector< TH1 * > &  obs 
) const [protected]

The log-likelihood formula from the PDG.

Parameters:
exp A set of expected spectra
obs The corresponding observed spectra
Returns:
The log-likelihood formula from the PDG

\[ \chi^2=2\sum_i^{\rm bins}\left(e_i-o_i+o_i\ln\left({o_i\over e_i}\right)\right) \]

Includes underflow and overflow bins and handles zero observed or expected events correctly.

Referenced by FindChiSqrForPars().

void NCExtrapolation::Prepare ( const Registry r  )  [virtual]

Read whatever values you need out of the registry to initialize yourself. Please remember to chain up to the NCExtrapolation implementation too.

Reimplemented in NCExtrapolationBeamMatrix, NCExtrapolationFarNear, NCExtrapolationNone, and NCExtrapolationPID.

Definition at line 91 of file NCExtrapolation.cxx.

References fAddEventsToEnergyBins, fCoordConv, fDoSystematicInterpolation, fEnergyThreshold, fHideFDData, fOscillationModel, fPrintFDEvents, fUseCC, fUseNC, fUseSystPenalty, Registry::Get(), Msg::kInfo, MSG, and NC::CoordinateConverter::Prepare().

00092 {
00093   MSG("NCExtrapolation", Msg::kInfo) << "NCExtrapolation::Prepare(const Registry& r)"
00094                                      << endl;
00095   int         tmpb;  // a temp bool.
00096   int         tmpi;  // a temp int.
00097   double      tmpd;  // a temp double.
00098 
00099   if(r.Get("OscillationModel",    tmpi)) fOscillationModel = NCType::EOscModel(tmpi);
00100 
00101   if(r.Get("EnergyThreshold",       tmpd)) fEnergyThreshold = tmpd;
00102 
00103   if(r.Get("UseCCEvents",           tmpb)) fUseCC           = tmpb;
00104   if(r.Get("UseNCEvents",           tmpb)) fUseNC           = tmpb;
00105 
00106   fCoordConv.Prepare(r);
00107 
00108   if(r.Get("PrintFDEvents", tmpb)) fPrintFDEvents = tmpb;
00109   if(r.Get("HideFDData",    tmpb)) fHideFDData    = tmpb;
00110 
00111   if(r.Get("DoSystematicInterpolation", tmpb))
00112     fDoSystematicInterpolation = tmpb;
00113 
00114   if(r.Get("DisableSystPenalty",    tmpb)) fUseSystPenalty = !tmpb;
00115 
00116   if(r.Get("AddEventsToEnergyBins", tmpb)) fAddEventsToEnergyBins = tmpb;
00117 }

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

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

This implementation clears the relevant NCBeam objects

Derived classes please override if you handle spectra yourself

Reimplemented in NCExtrapolationPID.

Definition at line 575 of file NCExtrapolation.cxx.

References fBeams, NCBeam::GetDetector(), it, Detector::kFar, and NCBeam::Reset().

Referenced by NCExtrapolationModule::CreateExtrapolations(), and NCExtrapolationBeamMatrix::FillNDHistsForXSectionFit().

00577 {
00578   for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
00579     NCBeam* beam = it->second;
00580     if(beam->GetDetector() == Detector::kFar){
00581       if(farData || farMC) beam->Reset(farData, farMC);
00582     }
00583     else{
00584       if(nearData || nearMC) beam->Reset(nearData, nearMC);
00585     }
00586   }
00587 }

void NCExtrapolation::SetBestFitCoordNDim ( const NC::Fitter::CoordNDim  minCoord  ) 

Set fMinCoord and prepare GetBestFitSysts and GetBestFitOscPars.

Definition at line 563 of file NCExtrapolation.cxx.

References fBestFitOscPars, fBestFitSysts, fCoordConv, fMinCoord, NC::CoordinateConverter::OscParsFromCoordNDim(), and NC::CoordinateConverter::SystParsFromCoordNDim().

Referenced by NC::FitMaster::Run().

00564 {
00565   fMinCoord = minCoord;
00566 
00567   if(fBestFitOscPars) delete fBestFitOscPars;
00568 
00569   fBestFitSysts = fCoordConv.SystParsFromCoordNDim(minCoord);
00570   fBestFitOscPars = fCoordConv.OscParsFromCoordNDim(minCoord);
00571 }

void NCExtrapolation::SetBestFitOscPars ( NC::OscProb::OscPars pars  )  [inline, protected]

Warning - doesn't affect fMinCoord.

Definition at line 185 of file NCExtrapolation.h.

References fBestFitOscPars.

00185 {fBestFitOscPars = pars;}

void NCExtrapolation::SetBestFitSysts ( NC::SystPars  systs  )  [inline, protected]

Warning - doesn't affect fMinCoord.

Definition at line 182 of file NCExtrapolation.h.

References fBestFitSysts.

00182 {fBestFitSysts = systs;}

vector< double > NCExtrapolation::SystsAsSigmas ( const NCBeam::Info info,
bool *  simpleScale = 0 
) const [protected]

Convert details of systematic shifts into a list of sigmas.

Parameters:
info Takes the shift to represent from here
[out] simpleScale Does the shift consist only of systematics that are simple scales, ie should be independent of oscillation parameters?
Returns:
Shifts represented as number of sigmas, for use as a key in NC::ISpectrumInterpolator

Definition at line 648 of file NCExtrapolation.cxx.

References fBeams, NCBeam::GetInfo(), NCBeam::Info::IsShifted(), it, NCType::kNormalization, NCType::kNumSystematicParameters, NCType::kParams, n, and NCType::ParamDef::sigma.

Referenced by InitializeInterpolator().

00650 {
00651   // Work out which systematics ever have shifts. Go through all beams and take
00652   // the combination of all shifts we see.
00653   bool shiftedVars[NCType::kNumSystematicParameters] = {false,};
00654   for(fBeams_t::const_iterator it = fBeams.begin(); it != fBeams.end(); ++it){
00655     NCBeam* beam = it->second;
00656     for(int n = 0; n < NCType::kNumSystematicParameters; ++n){
00657       if(beam->GetInfo().IsShifted(NCType::EFitParam(n)))
00658         shiftedVars[n] = true;
00659     }
00660   }
00661 
00662   // Does `shift` represent a systematic that is just a scale factor?
00663   // This is true for combinations of normalization and cleaning, at least
00664   // in PID and FarNear. It is possible that for some future extrapolation
00665   // more or less systematics can be treated this way.
00666   if(simpleScale) *simpleScale = true;
00667 
00668   // Translation of 'info' into number of sigma shifts in each parameter
00669   vector<double> shift;
00670   shift.reserve(NCType::kNumSystematicParameters);
00671   for(int n = 0; n < NCType::kNumSystematicParameters; ++n){
00672     if(shiftedVars[n]){
00673       double val;
00674       if(!info.IsShifted(NCType::EFitParam(n), &val)) val = 0;
00675       const double sigmas = val/NCType::kParams[n].sigma;
00676       // Currently the interpolator can't understand shifts that aren't +/-1
00677       assert(sigmas == 0 || sigmas == 1 || sigmas == -1);
00678       shift.push_back(sigmas);
00679 
00680       if(simpleScale &&
00681          sigmas != 0 &&
00682          n != NCType::kNormalization) *simpleScale = false; //&&
00683          //n != NCType::kNCNearClean) *simpleScale = false;
00684     } // end if shiftedVars
00685   } // end for n
00686 
00687   return shift;
00688 }

void NCExtrapolation::WriteResources ( const NC::OscProb::OscPars trueOscPars  )  [virtual]

gDirectory will point to the output file. Please remember to chain up to the NCExtrapolation implementation.

trueOscPars may be zero (eg obviously in case of fit to real data)

Reimplemented in NCExtrapolationBeamMatrix, NCExtrapolationFarNear, NCExtrapolationNone, and NCExtrapolationPID.

Definition at line 413 of file NCExtrapolation.cxx.

References DrawBestFitRatios(), DrawBestFitSpectra(), fBeams, fHideFDData, FillResultHistograms(), fInterpolator, NCBeam::GetDetector(), NCBeam::GetInfo(), GetShortName(), handler(), IgnoreErrors(), it, Detector::kFar, Msg::kInfo, Detector::kNear, NCBeam::MakeResultPlots(), MSG, NCBeam::WritePredictionResources(), NC::ISpectrumInterpolator::WriteResources(), and NCBeam::WriteResources().

00414 {
00415   MSG("NCExtrapolation", Msg::kInfo)
00416     << "NCExtrapolation::WriteResources() - "
00417     << GetShortName() << endl;
00418 
00419   // Make sure all those pretty NCBeam plots get filled
00420   for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it)
00421     FillResultHistograms(it->second);
00422 
00423   // get a pointer to the current directory
00424   // this is one of the output files
00425   TDirectory* saveDir = gDirectory;
00426 
00427   TDirectory* beamsDir = gDirectory->mkdir("beams", "beams");
00428   beamsDir->cd();
00429 
00430   //write out the spectra etc from the beam objects, while
00431   //checking on the prediction only flag
00432   for(fBeams_t::iterator it = fBeams.begin(); it != fBeams.end(); ++it){
00433     NCBeam* beam = it->second;
00434     TString dirName = beam->GetInfo().GetDescription();
00435 
00436     // Make the directory dirName, or if it already exists then cd into it.
00437     // There seems to be no way to do this without trigerring an error, so
00438     // tell ROOT to ignore them, and then restore the old state afterwards.
00439     ErrorHandlerFunc_t handler = GetErrorHandler();
00440     SetErrorHandler(IgnoreErrors);
00441     if(!beamsDir->cd(dirName)) beamsDir->mkdir(dirName)->cd();
00442     SetErrorHandler(handler);
00443 
00444     beam->MakeResultPlots();
00445     if(fHideFDData && beam->GetDetector() == Detector::kFar)
00446       beam->WritePredictionResources();
00447     else
00448       beam->WriteResources();
00449 
00450   }//end loop over beams
00451 
00452   saveDir->cd();
00453 
00454   DrawBestFitSpectra(Detector::kNear);
00455   DrawBestFitSpectra(Detector::kFar);
00456 
00457   DrawBestFitRatios(Detector::kNear);
00458   DrawBestFitRatios(Detector::kFar);
00459 
00460 
00461   if(fInterpolator){
00462     TDirectory* saveDir = gDirectory;
00463     gDirectory->mkdir("interp")->cd();
00464     fInterpolator->WriteResources();
00465     saveDir->cd();
00466   }
00467 
00468 
00469   saveDir->cd();
00470 }


Member Data Documentation

Enables us to disable filling expensive vectors in NCEnergyBin.

Definition at line 261 of file NCExtrapolation.h.

Referenced by GetBeam(), and Prepare().

Definition at line 277 of file NCExtrapolation.h.

Referenced by SetBestFitCoordNDim(), SetBestFitOscPars(), and ~NCExtrapolation().

Definition at line 276 of file NCExtrapolation.h.

Referenced by SetBestFitCoordNDim(), and SetBestFitSysts().

TODO - remove, is in FitMaster.

Definition at line 255 of file NCExtrapolation.h.

Referenced by InitializeInterpolator(), Prepare(), and NCExtrapolationBeamMatrix::WriteResources().

Don't accept any events below this energy.

Definition at line 227 of file NCExtrapolation.h.

Referenced by AddEvent(), NCExtrapolationPID::AddEvent(), and Prepare().

bool NCExtrapolation::fHideFDData [protected]

Preserve blinding by not showing FD data in plots.

Definition at line 200 of file NCExtrapolation.h.

Referenced by DrawBestFitRatios(), DrawBestFitSpectra(), Prepare(), and WriteResources().

Definition at line 291 of file NCExtrapolation.h.

Referenced by InitializeInterpolator(), and ~NCExtrapolation().

double NCExtrapolation::fMinChiSqr [private]

Definition at line 264 of file NCExtrapolation.h.

Referenced by ChiSqrAtBestFit().

This is the best fit point as found by the fitter.

Warning - the SetBestFitOscPars and SetBestFitSysts functions DO NOT modify this variable. It is only to be used if you are genuinely interested in the *fitter* output.

Definition at line 274 of file NCExtrapolation.h.

Referenced by GetBestFitCoordNDim(), GetBestFitOscPars(), GetBestFitSysts(), and SetBestFitCoordNDim().

Definition at line 202 of file NCExtrapolation.h.

Referenced by Prepare().

Whether to print all FD data events.

Definition at line 267 of file NCExtrapolation.h.

Referenced by AddEvent(), and Prepare().

Definition at line 225 of file NCExtrapolation.h.

bool NCExtrapolation::fUseCC [protected]
bool NCExtrapolation::fUseNC [protected]

Definition at line 258 of file NCExtrapolation.h.

Referenced by CalculateBestNormalization(), FindChiSqrForPars(), and Prepare().


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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1