NuFCFitterNSI Class Reference

#include <NuFCFitterNSI.h>

List of all members.

Public Member Functions

 NuFCFitterNSI (NuFCExperiment *experiment, NuMMHelperCPT *helper, NuMatrixSpectrum *ndNQData, NuMatrixSpectrum *ndPQData, const NuMMParameters &trueparams)
virtual ~NuFCFitterNSI ()
Double_t Fit (const NuMMParameters &fitParameters)
void AlwaysArchive (std::string filename)
 Call if you want to write out every experiment.
void ArchiveTo (std::string filename)
void Archive (std::string filename)

Private Member Functions

NuMMParameters GridSearch (TH3D &likesurf, NuMMParameters mmGrid, NuMMParameters *Ref=0)
 Does a grid search on an arbitrary surface.
NuMMParameters CoarseGridSearch (NuMMParameters start)
 Do a coarse grid search, to seed minuit.
NuMMParameters FineGridSearch (NuMMParameters start)
 Does a finer grid search, if minuit fails.
NuMMParameters MinuitFit (NuMMParameters start)
 Does a minuit fit.
Double_t DeltaChi (const NuMMParameters &bestfit)
 Calculates the Delta chi2 value.
NuMMParameters RecoverNegativeDeltaChi (const NuMMParameters &bestfit)
 Recovers a negative delta chi2 value.
void ClassifyFitPoint (const NuMMParameters &mmFit)
 Examines the fit point and classifies it.
void DoOutput (const NuMMParameters &mmFit, Double_t deltachi)
 Dos the summary output form the fitting.
void ShiftHistogram (TH3D &surface, Double_t shift)
 Shifts a histogram. Doesn't really belong here, but not sure where else.
 ClassDef (NuFCFitterNSI, 0)

Private Attributes

NuMMRunNSIfRun
 Pointer to the run object.
NuMatrixSpectrumffdNQData
NuMatrixSpectrumffdPQData
NuSystFitter fsystFitter
 Fitter object.
NuMMParametersfBestFit
 Copy of the best fit point.
Bool_t fFineGridSearch
Bool_t fNegativeDeltaChi
Bool_t fHighFitPoint
Bool_t fZeroFitPoint
Bool_t fArchive
std::string fArchiveFilename
const NuMMParameters fTrueParams

Detailed Description

Definition at line 31 of file NuFCFitterNSI.h.


Constructor & Destructor Documentation

NuFCFitterNSI::NuFCFitterNSI ( NuFCExperiment experiment,
NuMMHelperCPT helper,
NuMatrixSpectrum ndNQData,
NuMatrixSpectrum ndPQData,
const NuMMParameters trueparams 
)

Creates the fitter object. We have to accept pointers here, as everything else we want to use only accepts them, and we can't easily change without possibly breaking large amounts of existing code.

Parameters:
experiment The experiment (i.e. far data) that we are fitting.
helper The NuMMHelperCPT object that contains the helper data
ndNQData The near detector neutrino data
ndPQData The enar detector antineutrino data

Definition at line 25 of file NuFCFitterNSI.cxx.

References NuSystFitter::BatchModeOn(), ffdNQData, ffdPQData, fRun, fsystFitter, NuFCExperiment::GetNQSpectrum(), NuFCExperiment::GetPQSpectrum(), and NuSystFitter::push_back().

00031   : fRun(0),
00032     ffdNQData(0),
00033     ffdPQData(0),
00034     fsystFitter(),
00035     fBestFit(0),
00036     fFineGridSearch(false),
00037     fNegativeDeltaChi(false),
00038     fHighFitPoint(false),
00039     fZeroFitPoint(false),
00040     fArchive(false),
00041     fArchiveFilename("archive.root"),
00042     fTrueParams(trueparams)
00043 {
00044     // Firstly, get copies of the far detector spectrum
00045   //NuMatrixSpectrum *fdNQ = &experiment->GetSpectrum();
00046   ffdPQData = new NuMatrixSpectrum(experiment->GetPQSpectrum());
00047   ffdNQData = new NuMatrixSpectrum(experiment->GetNQSpectrum());
00048   
00049   // Create the run!
00050   fRun = new NuMMRunNSI(helper, ndNQData, ndPQData, ffdNQData, ffdPQData);
00051   fRun->PredictNus(false);
00052   fRun->SetMaximumEnergy(49.999);
00053   // fRun->QuietModeOn();
00054   
00055   // Now add this run to the fitter
00056   fsystFitter.push_back(fRun);
00057   //fsystFitter.QuietModeOn();
00058   fsystFitter.BatchModeOn();
00059 }

NuFCFitterNSI::~NuFCFitterNSI (  )  [virtual]

Definition at line 63 of file NuFCFitterNSI.cxx.

References fBestFit, ffdNQData, ffdPQData, and fRun.

00064 {
00065   // Delete our run object
00066   if (fRun) delete fRun;
00067   fRun = 0;
00068   
00069   // Delete our copies of the far detector data
00070   if (ffdNQData) delete ffdNQData;
00071   if (ffdPQData) delete ffdPQData;
00072   ffdNQData = ffdPQData = 0;
00073   
00074   // Delete the copy of the best fit point
00075   if (fBestFit) delete fBestFit;
00076   fBestFit = 0;
00077 }


Member Function Documentation

void NuFCFitterNSI::AlwaysArchive ( std::string  filename  )  [inline]

Call if you want to write out every experiment.

Definition at line 56 of file NuFCFitterNSI.h.

References fArchive, and fArchiveFilename.

00056                                          {
00057     fArchiveFilename = filename;
00058     fArchive = true;
00059   }

void NuFCFitterNSI::Archive ( std::string  filename  ) 

Definition at line 545 of file NuFCFitterNSI.cxx.

References MuELoss::e, fBestFit, ffdNQData, ffdPQData, NuSystFitter::FillLikelihoodSurfaceNSI3D(), fsystFitter, Msg::kError, Msg::kInfo, MSG, and NuMatrixSpectrum::Spectrum().

Referenced by Fit().

00546 {
00547   //MSG("NuFCFitterNSI",Msg::kInfo) << "Disabled all archiving." << endl;
00548   //return;
00549 
00550   if (!ffdPQData || !ffdNQData) {
00551     MSG("NUFCFitter",Msg::kError) << "No far detector spectrum: Cannot archive" << endl;
00552     return;
00553   }
00554   if (!fBestFit) {
00555     MSG("NuFCFitterNSI", Msg::kError)
00556       << "Error: Trying to archive before fit" << endl;
00557     return;
00558   }
00559   
00560   MSG("NuFCFitterNSI", Msg::kInfo)
00561     << "    Archiving to file " << filename << "..." << endl;
00562 
00563   // The actual experiment spectrum
00564 //  TH1D fdSpec(*(ffdPQData->Spectrum()));
00565   NuMatrixSpectrum fdSpec(*ffdPQData);
00566 //   // fdSpec.RebinTo4GeV();
00567 //   // Grab the name out of this
00568   string basename = fdSpec.Spectrum()->GetName();
00569 // 
00570 // // No longer do the likelihood surface - it takes too long. Instead just write out the spectrum,
00571 // // and we can regenerate the surface later if we wish
00572 //   
00573   // Now we want the likelihood surfaces. First do the linear portion
00574   string surfname = basename + "_surf";
00575   TH3D likesurface(surfname.c_str(), surfname.c_str(), 50, 0, 1, 100, 0, 20e-3, 25, 0, 0.5);
00576   fsystFitter.FillLikelihoodSurfaceNSI3D(likesurface, *fBestFit);
00577 //   
00578 //   // Now do the logarithmic plot
00579 //   vector<Double_t> bins;
00580 //   for (Int_t i = -204; i <= 0; i++) {
00581 //     bins.push_back(pow(10, (Double_t)i/100.0));
00582 //   }
00583 //   Int_t bincount = bins.size()-1;
00584 //   string logname = surfname + "_log";
00585 //   TH2D logsurf(logname.c_str(), logname.c_str(), 50, 0, 1, bincount, &(bins.front()));
00586 //   Double_t minlog = fsystFitter.FillLikelihoodSurfaceBar(logsurf, *fBestFit);
00587 //   
00588 //   // Now, shift both histograms by the minimum
00589 //   Double_t surfmin = minlin < minlog ? minlin : minlog;
00590 //   ShiftHistogram(likesurface, -surfmin);
00591 //   ShiftHistogram(logsurf, -surfmin);
00592   
00593   // Now, write these out!
00594   TFile tf(filename.c_str(), "UPDATE");
00595     fdSpec.Spectrum()->Write();
00596     likesurface.Write();
00597     // logsurf.Write();
00598   tf.Close();
00599 }

void NuFCFitterNSI::ArchiveTo ( std::string  filename  )  [inline]

Definition at line 61 of file NuFCFitterNSI.h.

References fArchiveFilename.

00061                                      {
00062     fArchiveFilename = filename;
00063   }

NuFCFitterNSI::ClassDef ( NuFCFitterNSI  ,
 
) [private]
void NuFCFitterNSI::ClassifyFitPoint ( const NuMMParameters mmFit  )  [private]

Examines the fit point and classifies it.

Definition at line 519 of file NuFCFitterNSI.cxx.

References NuMMParameters::Dm2Bar(), MuELoss::e, NuMMParameters::Epsilon(), fHighFitPoint, fTrueParams, fZeroFitPoint, Msg::kInfo, NuMMParameters::LowerDm2BarLimit(), NuMMParameters::LowerEpsilonLimit(), NuMMParameters::LowerSn2BarLimit(), MSG, and NuMMParameters::Sn2Bar().

Referenced by Fit().

00520 {
00521   // cout << "    Checking against minimum dm: " << fTrueParams.LowerDm2BarLimit() << ", sn: " << fTrueParams.LowerSn2BarLimit() << endl;
00522  
00523   // Check to see if the fit is to 'zero'
00524   if (mmFit.Sn2Bar() <= fTrueParams.LowerSn2BarLimit() + 1e-4 || 
00525       mmFit.Dm2Bar() <= fTrueParams.LowerDm2BarLimit() + 1e-6 ||
00526       mmFit.Epsilon() <= fTrueParams.LowerEpsilonLimit() + 1e-4)
00527   {
00528       MSG("NuFCFitterNSI", Msg::kInfo)
00529         << "    Best fit is to no oscillation case\n";
00530       fZeroFitPoint = true;
00531   }
00532     
00533   // Now check if the best fit is maxing out
00534   // This message is a little out of date, but will do for now
00535   if (mmFit.Dm2Bar() >= 30e-3 - 0.25e-4) {
00536     MSG("NuFCFitterNSI", Msg::kInfo)
00537      << "    Best fit is maxing out on dm2bar ("
00538      << mmFit.Dm2Bar() << ")" << endl;
00539     fHighFitPoint = true;
00540   }
00541 }

NuMMParameters NuFCFitterNSI::CoarseGridSearch ( NuMMParameters  start  )  [private]

Do a coarse grid search, to seed minuit.

Definition at line 208 of file NuFCFitterNSI.cxx.

References NuMMParameters::Dm2Bar(), MuELoss::e, NuMMParameters::Epsilon(), NuFCConfig::FixedSin(), fTrueParams, NuFCConfig::GetConfig(), GridSearch(), NuMMParameters::IsDm2BarConstrained(), NuMMParameters::IsEpsilonConstrained(), NuMMParameters::IsSn2BarConstrained(), Msg::kInfo, NuMMParameters::LowerDm2BarLimit(), NuMMParameters::LowerEpsilonLimit(), NuMMParameters::LowerSn2BarLimit(), max, MAXMSG, min, MSG, NuMMParameters::Sn2Bar(), NuMMParameters::UpperDm2BarLimit(), NuMMParameters::UpperEpsilonLimit(), and NuMMParameters::UpperSn2BarLimit().

Referenced by Fit().

00209 {
00210   // The width and height of the grid
00211   Int_t gridWidth = 5;
00212   Double_t sinMin = 0.2;
00213   Double_t sinMax = 1.0;
00214   Double_t dmMin = 2e-3;
00215   Double_t dmMax = 30e-3;
00216   Double_t epsilonMin = 0.0;
00217   Double_t epsilonMax = 0.5;
00218 
00219   // Constrain to the tighter of mmGrid constraints
00220   // and the defaults above
00221   if (mmGrid.IsSn2BarConstrained()) {
00222     sinMin = max(sinMin, mmGrid.LowerSn2BarLimit());
00223     sinMax = min(sinMax, mmGrid.UpperSn2BarLimit());
00224     MAXMSG("NuFCFitterNSI",Msg::kInfo,1) << "Performing constrained grid search in sin: " << sinMin << " - " << sinMax << endl;
00225   }
00226   if (mmGrid.IsDm2BarConstrained()) {
00227     dmMin = max(dmMin, mmGrid.LowerDm2BarLimit());
00228     dmMax = min(dmMax, mmGrid.UpperDm2BarLimit());
00229     MAXMSG("NuFCFitterNSI",Msg::kInfo,1) << "Performing constrained grid search in dm2: " << dmMin << " - " << dmMax << endl;
00230   }
00231   if (mmGrid.IsEpsilonConstrained()) {
00232     epsilonMin = max(epsilonMin, mmGrid.LowerEpsilonLimit());
00233     epsilonMax = min(epsilonMax, mmGrid.UpperEpsilonLimit());
00234     MAXMSG("NuFCFitterNSI",Msg::kInfo,1) << "Performing constrained grid search in epsilon: " << epsilonMin << " - " << epsilonMax << endl;
00235   }
00236   
00237   const Int_t gridHeight = 15;
00238   const Int_t gridDepth = 10;
00239   // Check if we are doing one-parameter fit, and if so constrain the grid search
00240   // to one-dimension
00241   NuFCConfig &cfg = NuFCConfig::GetConfig();
00242   if (cfg.FixedSin()) {
00243     gridWidth = 1;
00244     sinMin = sinMax = fTrueParams.Sn2Bar();
00245   }
00246   
00247   cout << "About to generate TH3D surface with epsilon: " << epsilonMin << " - " << epsilonMax << " and gridDepth: " << gridDepth << endl;
00248   // Generate the surface
00249   TH3D likesurf("likesurf", "likesurf", gridWidth, sinMin, sinMax, gridHeight, dmMin, dmMax, gridDepth, epsilonMin, epsilonMax);
00250   cout << "Generated likesurf. " << endl;
00251   mmGrid = GridSearch(likesurf, mmGrid);
00252 
00253   cout<< "After mmGrid = GridSearch()" << endl;
00254   // Now do a log search
00255   vector<Double_t> bins;
00256   for (Int_t i = -15; i <= -1; i++) {
00257     bins.push_back(pow(10, (Double_t)i/10.0));
00258   }
00259   Int_t bincount = bins.size()-1;
00260 
00261   vector<Double_t> sinbins;
00262   for(Int_t i = 0; i<=gridWidth; i++) {
00263     sinbins.push_back(i*(sinMax-sinMin)/gridWidth);
00264   }
00265   Int_t sinbincount = sinbins.size()-1;
00266 
00267   vector<Double_t> epsbins;
00268   for(Int_t i = 0; i<=gridDepth; i++) {
00269     epsbins.push_back(i*(epsilonMax-epsilonMin)/gridDepth);
00270   }
00271   Int_t epsbincount = epsbins.size()-1;
00272 
00273 
00274   TH3D logsurf("logsurf", "logsurf", sinbincount, &(sinbins.front()), 
00275                bincount, &(bins.front()), epsbincount, &(epsbins.front()));
00276   mmGrid = GridSearch(logsurf, mmGrid, &mmGrid);
00277 
00278   // Output the minimum
00279   MSG("NuFCFitterNSI", Msg::kInfo) << "    Coarse grid search found sn2: "
00280                                    << mmGrid.Sn2Bar() << ", dm: " << mmGrid.Dm2Bar() 
00281                                    << ", epsilon: " << mmGrid.Epsilon() << '\n';
00282   
00283   return mmGrid;
00284 }

Double_t NuFCFitterNSI::DeltaChi ( const NuMMParameters bestfit  )  [private]

Calculates the Delta chi2 value.

Definition at line 431 of file NuFCFitterNSI.cxx.

References NuMMParameters::Dm2Bar(), NuMMParameters::Epsilon(), fsystFitter, fTrueParams, Msg::kDebug, NuMMParameters::LowerDm2BarLimit(), NuMMParameters::LowerEpsilonLimit(), NuMMParameters::LowerSn2BarLimit(), MSG, NuMMParameters::Sn2Bar(), and NuMMParameters::VectorParameters().

Referenced by Fit().

00432 {
00433     // Calculate the two likelihood points
00434     Double_t fitlike = fsystFitter(bestfit.VectorParameters());
00435 
00436     Double_t truelikelihood = 0;  
00437 
00438     // cout << "bfdm: " << bestfit.Dm2Bar() << ", limit = " << fTrueParams.LowerDm2BarLimit() << endl;
00439 
00440     // Don't allow the parameters to go below constrained - this can
00441     // give us an artificially negative delta chi2
00442     if ( (fTrueParams.Dm2Bar() < fTrueParams.LowerDm2BarLimit()) ||
00443          (fTrueParams.Sn2Bar() < fTrueParams.LowerSn2BarLimit()) ||
00444          (fTrueParams.Epsilon() < fTrueParams.LowerEpsilonLimit()) ) {  
00445         MSG("NuFCFitterNSI", Msg::kDebug) << "    True point is lower than fit is allowed." << endl;
00446         NuMMParameters pc(bestfit);
00447         pc.Dm2Bar(fTrueParams.LowerDm2BarLimit());
00448         pc.Sn2Bar(fTrueParams.LowerSn2BarLimit());
00449         pc.Epsilon(fTrueParams.LowerEpsilonLimit());
00450         truelikelihood = fsystFitter(pc.VectorParameters());
00451     } else {
00452         truelikelihood = fsystFitter(fTrueParams.VectorParameters());
00453     }
00454     
00455 
00456     // fmmPars.Dm2Bar(fmmPars.LowerDm2BarLimit());
00457     //   MAXMSG("NuFCGriPoint", 5, Msg::kInfo) << "Moving true grid point "
00458     //   fmmPars.Sn2Bar(fmmPars.LowerSn2BarLimit());
00459     //   
00460     // Double_t truelikelihood = fsystFitter(fTrueParams.VectorParameters());
00461 
00462     // Return the difference
00463     return truelikelihood - fitlike;
00464 }

void NuFCFitterNSI::DoOutput ( const NuMMParameters mmFit,
Double_t  deltachi 
) [private]

Dos the summary output form the fitting.

Definition at line 125 of file NuFCFitterNSI.cxx.

References NuMMParameters::Dm2Bar(), NuMMParameters::Epsilon(), fFineGridSearch, fNegativeDeltaChi, fsystFitter, Msg::kInfo, MSG, NuMMParameters::Sn2Bar(), and NuMMParameters::VectorParameters().

Referenced by Fit().

00126 {
00127   // Output the fit results
00128   Double_t fitlike = fsystFitter(mmFit.VectorParameters());
00129   MSG("NuFCFitterNSI", Msg::kInfo)
00130     << "    Minimum is dm2bar: " << mmFit.Dm2Bar()
00131     << " sn2bar: " << mmFit.Sn2Bar()
00132     << " epsilon: " << mmFit.Epsilon() << '\n'
00133     << "    Fit Likelyhood: " << fitlike << '\n'
00134     << "    Delta Chi2:     " << deltachi << '\n';
00135   
00136   // Did we want to ask the control script to log this?
00137   if (fFineGridSearch || fNegativeDeltaChi) {
00138     MSG("NuFCFitterNSI", Msg::kInfo) << "    Please Log this experiment." << endl;
00139   }
00140   
00141 }

NuMMParameters NuFCFitterNSI::FineGridSearch ( NuMMParameters  start  )  [private]

Does a finer grid search, if minuit fails.

Definition at line 288 of file NuFCFitterNSI.cxx.

References NuMMParameters::Dm2Bar(), MuELoss::e, NuMMParameters::Epsilon(), NuFCConfig::FixedSin(), fTrueParams, NuFCConfig::GetConfig(), GridSearch(), Msg::kInfo, NuMMParameters::LowerDm2BarLimit(), NuMMParameters::LowerEpsilonLimit(), NuMMParameters::LowerSn2BarLimit(), MSG, and NuMMParameters::Sn2Bar().

Referenced by MinuitFit().

00289 {
00290   // Calculate the location of the fine grid search
00291   Double_t gridL, gridR, gridT, gridB, gridU, gridD;
00292   gridL = gridR = gridT = gridB = gridU = gridD = -1.0;
00293   
00294   // The width and height of the grid
00295   Int_t gridWidth = 19;
00296   const Int_t gridHeight = 19;
00297   const Int_t gridDepth = 19;
00298 
00299   // Configuration
00300   NuFCConfig &cfg = NuFCConfig::GetConfig();
00301   
00302   // change behavior if we want a one-parameter fit
00303   if (!cfg.FixedSin()) {
00304       // Go +- 0.2 in sin space
00305     gridL = mmStart.Sn2Bar() - 0.2;
00306     gridR = gridL + 0.4;
00307     // Now limit it to the physical realm
00308     if (gridL <= fTrueParams.LowerSn2BarLimit()) gridL = fTrueParams.LowerSn2BarLimit();
00309     if (gridR > 1.0) gridR = 1.0;
00310   } else {
00311     // One-parameter fit
00312     gridL = gridR = fTrueParams.Sn2Bar();
00313     gridWidth = 1;
00314   }
00315 
00316   // And go +-2e-3 in mass space
00317   gridB = mmStart.Dm2Bar() - 2e-3;
00318   gridT = gridB + 4e-3;
00319   // We only really need to worry about lower limit here
00320   if (gridB <= fTrueParams.LowerDm2BarLimit()) gridB = fTrueParams.LowerDm2BarLimit();
00321 
00322   // Go +-0.1 in epsilon space
00323   gridD = mmStart.Epsilon() - 0.1;
00324   gridU = gridD + 0.2;
00325   if(gridD <= fTrueParams.LowerEpsilonLimit()) gridD = fTrueParams.LowerEpsilonLimit();
00326 
00327   MSG("NuFCFitterNSI", Msg::kInfo) 
00328     << "    Doing fine grid search on sin = [ "
00329     << gridL << ", " << gridR << " ]\n"
00330     << "                              dm  = [ "
00331     << gridB << ", " << gridT << " ]\n"
00332     << "                              epsilon  = [ "
00333     << gridD << ", " << gridU << " ]\n";
00334 
00335   // Ensure that we haven't done anything stupid like set low above high
00336   if (gridL > gridR || gridT <= gridB || gridD > gridU) {
00337 //      cout << "Error: Negative range in fine grid search" << endl;
00338       throw runtime_error("Negative range in fine grid search");
00339   }
00340   
00341   // Make the surface
00342   
00343   // Generate the surface
00344   TH3D likesurf("likesurf", "likesurf", gridWidth, gridL, gridR,
00345     gridHeight, gridB, gridT, gridDepth, gridD, gridU);
00346     
00347   // And now do a grid search on this
00348   NuMMParameters mmGrid = GridSearch(likesurf, mmStart);
00349   
00350   // Temporarily, archive any fine grid surface
00351   // TFile tf("archives/finegrid.root", "UPDATE");
00352   //   likesurf.Write();
00353   //   tf.Close();
00354   
00355   // Output the minimum
00356   MSG("NuFCFitterNSI", Msg::kInfo) << "    Fine grid search found sn2: "
00357                                    << mmGrid.Sn2Bar() << ", dm: " 
00358                                    << mmGrid.Dm2Bar() <<", epsilon: " 
00359                                    << mmGrid.Epsilon() << '\n';
00360   
00361   
00362   return mmGrid;
00363 }

Double_t NuFCFitterNSI::Fit ( const NuMMParameters fitParameters  ) 

Fits the previously passed in FC experiment. Returns an std::pair with the delta chi2, and the best fit

Parameters:
fitParameters a parameters object, set to the true value

Definition at line 81 of file NuFCFitterNSI.cxx.

References Archive(), ClassifyFitPoint(), CoarseGridSearch(), DeltaChi(), DoOutput(), MuELoss::e, fArchive, fArchiveFilename, fBestFit, MinuitFit(), and RecoverNegativeDeltaChi().

00082 {
00083   // Start with a coarse grid search, to seed minuit
00084   const NuMMParameters mmGrid = CoarseGridSearch(fitParameters);    
00085   cout << "Finished CoarseGridSearch, starting MinuitFit..." << endl;
00086 
00087   // Now do a minuit fit
00088   NuMMParameters mmFit = MinuitFit(mmGrid);
00089   cout << "Finished MinuitFit" << endl;
00090 
00091   Double_t dchi = DeltaChi(mmFit);
00092 
00093   // If this gave a negative value, rescue it (allow a tiny gap)
00094   if (dchi < 0) {
00095     if (dchi < -1e-5) {
00096       mmFit = RecoverNegativeDeltaChi(mmFit);
00097       dchi = DeltaChi(mmFit);
00098     } else {
00099       // It is below zero, but within error. Set it to zero!
00100       dchi = 0;
00101     }
00102   }
00103 
00104   // Make a copy of the best fit point
00105   if (fBestFit) delete fBestFit;
00106   fBestFit = new NuMMParameters(mmFit);
00107 
00108   // Classify the best fit point i.e. is it zero, high etc
00109   // this is mainly for informational purposes (i.e. the run script)
00110   ClassifyFitPoint(mmFit);
00111 
00112   // Do we want to archive?
00113   if (fArchive) { // || fNegativeDeltaChi
00114     Archive(fArchiveFilename);
00115   }
00116 
00117   // Do the output
00118   DoOutput(mmFit, dchi);
00119   
00120   return dchi;
00121 }

NuMMParameters NuFCFitterNSI::GridSearch ( TH3D &  likesurf,
NuMMParameters  mmGrid,
NuMMParameters Ref = 0 
) [private]

Does a grid search on an arbitrary surface.

Definition at line 367 of file NuFCFitterNSI.cxx.

References NuMMParameters::Dm2Bar(), NuMMParameters::Epsilon(), NuSystFitter::FillLikelihoodSurfaceNSI3D(), NuFCConfig::FixedSin(), fsystFitter, fTrueParams, NuFCConfig::GetConfig(), Msg::kDebug, MSG, NuMMParameters::Sn2Bar(), and NuMMParameters::VectorParameters().

Referenced by CoarseGridSearch(), and FineGridSearch().

00368 {
00369   
00370     Double_t likelihood = fsystFitter.FillLikelihoodSurfaceNSI3D(likesurf, mmGrid);
00371     //cout << "Filled likelihood surface" << endl;
00372     // Find the point on this surface that is minimum
00373     Int_t minbinX = 0, minbinY = 0, minbinZ = 0;
00374     likesurf.GetMinimumBin(minbinX, minbinY, minbinZ);
00375     Double_t minsin = likesurf.GetXaxis()->GetBinCenter(minbinX);
00376     Double_t minmass = likesurf.GetYaxis()->GetBinCenter(minbinY);
00377     Double_t mineps = likesurf.GetZaxis()->GetBinCenter(minbinZ);
00378     cout << "minsin: " << minsin << ", minmass: " << minmass
00379          << ", mineps: " << mineps << endl;
00380     // Access the configuration
00381     NuFCConfig &cfg = NuFCConfig::GetConfig();
00382 
00383     // Check the zero point, to see if it is lower than anything else on the grid
00384     mmGrid.Dm2Bar(0);
00385     // Different 'zero' point if one parameter
00386     if (cfg.FixedSin()) {
00387       mmGrid.Sn2Bar(fTrueParams.Sn2Bar());
00388     } else {
00389       mmGrid.Sn2Bar(0);
00390     }
00391     
00392     cout << ", mmGrid.VectorParameters().size(): " 
00393          << mmGrid.VectorParameters().size() << endl;
00394 
00395    Double_t zeropoint = fsystFitter(mmGrid.VectorParameters());
00396     cout << "zeropoint: " << zeropoint << endl;
00397     if (zeropoint < likelihood) {
00398         // Then the no-oscillation case is the minimum
00399         minsin = 0.0;
00400         minmass = 0.0;
00401         mineps = 0.0;
00402         
00403         if (cfg.FixedSin()) {
00404           minsin = fTrueParams.Sn2Bar();
00405         }
00406     }
00407     
00408     cout << "zeropoint: " << zeropoint << endl;
00409 
00410     // If we have been given a reference point, check to see if this is lower
00411     if (Ref) {
00412       Double_t refpoint = fsystFitter(Ref->VectorParameters());
00413       if (refpoint < likelihood) {
00414         MSG("NuFCFitterNSI", Msg::kDebug) << "    Reference point is lower than grid search" << endl;
00415         minsin = Ref->Sn2Bar();
00416         minmass = Ref->Dm2Bar();
00417         mineps = Ref->Epsilon();
00418       }
00419     }
00420     
00421     // Set mmPars to the newly found point, then return it
00422     mmGrid.Sn2Bar(minsin);
00423     mmGrid.Dm2Bar(minmass);
00424     mmGrid.Epsilon(mineps);
00425 
00426     return mmGrid;
00427 }

NuMMParameters NuFCFitterNSI::MinuitFit ( NuMMParameters  start  )  [private]

Does a minuit fit.

Definition at line 145 of file NuFCFitterNSI.cxx.

References NuMMParameters::Dm2(), NuMMParameters::Dm2Bar(), NuMMParameters::Epsilon(), fFineGridSearch, FineGridSearch(), NuFCConfig::FixedSin(), NuMMParameters::FixSn2Bar(), fsystFitter, fTrueParams, NuFCConfig::GetConfig(), Msg::kDebug, Msg::kInfo, NuMMParameters::LowerDm2BarLimit(), NuMMParameters::LowerEpsilonLimit(), NuMMParameters::LowerSn2BarLimit(), NuSystFitter::Minimise(), NuMMParameters::MinuitPass(), MSG, NuMMParameters::Sn2(), NuMMParameters::Sn2Bar(), and NuMMParameters::VectorParameters().

Referenced by Fit().

00146 {
00147   // If it is trying to start us at zero, reset to the search minimum
00148   // fTrueParams.LowerDm2BarLimit()) {
00149   if (mmStart.Sn2Bar() == 0) {
00150     MSG("NuFCFitterNSI", Msg::kDebug) << "    Starting minuit at no-oscillation; shifting to lower limits" << endl;
00151     mmStart.Sn2Bar(fTrueParams.LowerDm2BarLimit());
00152     mmStart.Dm2Bar(fTrueParams.LowerSn2BarLimit());
00153     mmStart.Epsilon(fTrueParams.LowerEpsilonLimit());
00154   }
00155   
00156   // Fix the sin, if specified
00157   NuFCConfig &cfg = NuFCConfig::GetConfig();
00158   if (cfg.FixedSin()) {
00159     mmStart.FixSn2Bar();
00160     mmStart.Sn2Bar(fTrueParams.Sn2Bar());
00161   }
00162   
00163   // Now run minuit itself
00164   NuMMParameters mmFit = fsystFitter.Minimise(mmStart);
00165   
00166   // Did minuit succeed?
00167   if (mmFit.MinuitPass())
00168   {
00169     // Tell them what minuit found
00170     MSG("NuFCFitterNSI", Msg::kInfo)
00171       << "    Minuit found: dm2: " << mmFit.Dm2()
00172       << " sn2: " << mmFit.Sn2()
00173       << " dm2bar: " << mmFit.Dm2Bar()
00174       << " sn2bar: " << mmFit.Sn2Bar()
00175       << " epsilon: " << mmFit.Epsilon()
00176       << '\n';
00177   } else {
00178       // Minuit failed to fit. Use a more intensive grid search,
00179       // starting on the minuit starting point (that we get from
00180       // a grid search)
00181       MSG("NuFCFitterNSI", Msg::kInfo) << "    Minuit fit failed." << endl;
00182       MSG("NuFCFitterNSI", Msg::kInfo) << "    Searching fine grid around sn2: " << mmStart.Sn2Bar()
00183                                        << ", dm2: " << mmStart.Dm2Bar()
00184                                        << ", epsilon: " << mmStart.Epsilon() << endl;
00185       NuMMParameters mmFineGrid = FineGridSearch(mmStart);
00186       fFineGridSearch = true;
00187       
00188       // Test the fine grid search against the result from minuit = it could be possible
00189       // that minuit found a 'better' point, even if it failed
00190       Double_t minuit_chi = fsystFitter(mmFit.VectorParameters());
00191       Double_t finegrid_chi = fsystFitter(mmFineGrid.VectorParameters());
00192       
00193       if (minuit_chi > finegrid_chi) {
00194         // The new, fine grid point is better
00195         mmFit = mmFineGrid;
00196       } else {
00197         // The minuit fail point is actually better. Use that.
00198         MSG("NuFCFitterNSI", Msg::kInfo) << "    Failed Minuit fit is better than fine grid, using." << endl;
00199       }
00200       
00201   }
00202   
00203   return mmFit;
00204 }

NuMMParameters NuFCFitterNSI::RecoverNegativeDeltaChi ( const NuMMParameters bestfit  )  [private]

Recovers a negative delta chi2 value.

Definition at line 469 of file NuFCFitterNSI.cxx.

References NuMMParameters::Dm2Bar(), NuMMParameters::Epsilon(), NuFCConfig::FixedSin(), NuMMParameters::FixSn2Bar(), fNegativeDeltaChi, fsystFitter, fTrueParams, NuFCConfig::GetConfig(), Msg::kError, Msg::kInfo, NuSystFitter::Minimise(), NuMMParameters::MinuitPass(), MSG, NuMMParameters::Sn2Bar(), and NuMMParameters::VectorParameters().

Referenced by Fit().

00470 {
00471   // Set the failure flag
00472   fNegativeDeltaChi = true;
00473 
00474   // Recalculate the likelihoods
00475   Double_t truelike = fsystFitter(fTrueParams.VectorParameters());
00476   Double_t oldlike = fsystFitter(bestfit.VectorParameters());;
00477 
00478   // Make a copy of the true parameters, to ensure that we do a 1D fit if specified
00479   // (this is assuming we cannot guarantee that the true parameters has the parameter fixed)
00480   NuMMParameters tpcopy(fTrueParams);
00481   NuFCConfig &cfg = NuFCConfig::GetConfig();
00482   if (cfg.FixedSin()) {
00483     tpcopy.FixSn2Bar();
00484   }
00485   // NuMMParameters newFit = fsystFitter.Minimise(fTrueParams);
00486   NuMMParameters newFit = fsystFitter.Minimise(tpcopy);
00487   
00488   MSG("NuFCFitterNSI", Msg::kInfo)
00489     << "    Negative delta chi2 likelihood (" << truelike - oldlike
00490     << ") found. Refitting from true (" << truelike << ")" << endl;
00491   
00492   cout << "newFit.MinuitPass(): " << newFit.MinuitPass() << endl;
00493   
00494   // If this passed, and is lower than our old likelihood, then use it
00495   if (newFit.MinuitPass())
00496   {
00497       MSG("NuFCFitterNSI", Msg::kInfo) << "    Refit passed." << endl;
00498       Double_t newlike = fsystFitter(newFit.VectorParameters());
00499       if (newlike > oldlike) {
00500           MSG("NuFCFitterNSI", Msg::kInfo)
00501             << "    New fit point ( dm: " << newFit.Dm2Bar() << ", sn: " << newFit.Sn2Bar() << ", epsilon: " << newFit.Epsilon() << " )\n"
00502             << "    at " << newlike << " is higher than before ( " << oldlike << " )! strange!\n"
00503             << "    Returning true fit point, as is lowest we have seen."
00504             << endl;
00505           newFit = fTrueParams;
00506       }
00507   } else {
00508       MSG("NuFCFitterNSI", Msg::kError)
00509         << "    Unrecoverable error: Refit on negative deltachi failed."
00510         << endl;
00511       throw logic_error("NuFCFitterNSI Refit on negative deltachi failed.");
00512   }
00513 
00514   return newFit;
00515 }

void NuFCFitterNSI::ShiftHistogram ( TH3D &  surface,
Double_t  shift 
) [private]

Shifts a histogram. Doesn't really belong here, but not sure where else.

Definition at line 603 of file NuFCFitterNSI.cxx.

00604 {
00605     for (Int_t x = 1; x <= surface.GetNbinsX(); x++)
00606     {
00607         for (Int_t y = 1; y <= surface.GetNbinsY(); y++)
00608         {
00609           for (Int_t z = 1; z <= surface.GetNbinsZ(); z++)
00610             {
00611             Double_t contents = surface.GetBinContent(x, y, z);
00612             surface.SetBinContent(x, y, z, contents + shift);
00613             }
00614         }
00615     }
00616 }


Member Data Documentation

Bool_t NuFCFitterNSI::fArchive [private]

Definition at line 115 of file NuFCFitterNSI.h.

Referenced by AlwaysArchive(), and Fit().

std::string NuFCFitterNSI::fArchiveFilename [private]

Definition at line 116 of file NuFCFitterNSI.h.

Referenced by AlwaysArchive(), ArchiveTo(), and Fit().

Copy of the best fit point.

Definition at line 106 of file NuFCFitterNSI.h.

Referenced by Archive(), Fit(), and ~NuFCFitterNSI().

Pointers to copies of the far detector spectrum - the NuMMRunFC objects only accept pointers, but don't alter or copy them.

Definition at line 99 of file NuFCFitterNSI.h.

Referenced by Archive(), NuFCFitterNSI(), and ~NuFCFitterNSI().

Definition at line 100 of file NuFCFitterNSI.h.

Referenced by Archive(), NuFCFitterNSI(), and ~NuFCFitterNSI().

Definition at line 109 of file NuFCFitterNSI.h.

Referenced by DoOutput(), and MinuitFit().

Bool_t NuFCFitterNSI::fHighFitPoint [private]

Definition at line 112 of file NuFCFitterNSI.h.

Referenced by ClassifyFitPoint().

Definition at line 110 of file NuFCFitterNSI.h.

Referenced by DoOutput(), and RecoverNegativeDeltaChi().

Pointer to the run object.

Definition at line 95 of file NuFCFitterNSI.h.

Referenced by NuFCFitterNSI(), and ~NuFCFitterNSI().

Fitter object.

Definition at line 103 of file NuFCFitterNSI.h.

Referenced by Archive(), DeltaChi(), DoOutput(), GridSearch(), MinuitFit(), NuFCFitterNSI(), and RecoverNegativeDeltaChi().

Bool_t NuFCFitterNSI::fZeroFitPoint [private]

Definition at line 113 of file NuFCFitterNSI.h.

Referenced by ClassifyFitPoint().


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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1