GhostFramework Class Reference

#include <GhostFramework.h>

List of all members.

Public Member Functions

 ~GhostFramework ()
void UseResolutionLikelihood (int reslik=1)
void UseStandardLikelihood (int stdlik=1)
void AddSample (GhostSample *sample)
void RunFit ()
Double_t LogLikelihood (Double_t sin2_par, Double_t dm2_par)
Double_t LogLikelihoodNuisance (Double_t *par)
void Pause ()
void SetChi2Plot (Bool_t plot)

Static Public Member Functions

static GhostFrameworkInstance (GhostInputs *inputs, GhostFakeData *fake_data)

Private Member Functions

void RebinRealData ()
void FindInitialNuisanceValues ()
void OscFitGridSearch ()
 GhostFramework (GhostInputs *inputs, GhostFakeData *fake_data)

Private Attributes

GhostInputsgInputs
GhostFakeDatagFakeData
Bool_t fFitDecoherence
Bool_t fFitDecay
Bool_t fResolutionLikelihood
Bool_t fStandardLikelihood
Bool_t fRunParallel
Bool_t fFitScales
Bool_t fMinuitInitialValues
Int_t fNumMockDataExpts
Int_t fNumSamples
Double_t fDm2BinWidth
Double_t fSin2BinWidth
Double_t fSin2MinOverall
Double_t fSin2MaxOverall
Double_t fDm2MinOverall
Double_t fDm2MaxOverall
Int_t fNumDm2BinsOverall
Int_t fNumSin2BinsOverall
Int_t fNum2DSysts
Int_t fNum1DSysts
Int_t fNumSin2Bins
Int_t fNumDm2Bins
Double_t fSin2Min
Double_t fSin2Max
Double_t fDm2Min
Double_t fDm2Max
Int_t fNumScalingFactors
Bool_t fNuisanceParameters
Double_t fCCDetSigma
Double_t fCCRockSigma
Double_t fNCSigma
Double_t fNuMuBarSigma
Double_t fCCDetScale
Double_t fCCRockScale
Double_t fNCScale
Double_t fNuMuBarScale
Bool_t fFitCCDet
Bool_t fFitCCRock
Bool_t fFitNC
Bool_t fFitNuMuBar
Double_t fSin2_par
Double_t fDm2_par
Double_t * fInitialScalingFactors
Double_t * fMinScale
Double_t * fMinScaleError
Double_t * f2DInitialSystematicShift
Double_t * f2DMinSystematicShift
Double_t * f2DMinErrorSystematicShift
Double_t * f1DInitialSystematicShift
Double_t * f1DMinSystematicShift
Double_t * f1DMinErrorSystematicShift
Double_t fMu2_Initial
Double_t fMu2_High
Double_t fMu2_par
Double_t fMu2_Error
Int_t fSurfacePair
Int_t fDecoPower
TMinuit * fInitialNMinuit
TMinuit * fFreeMinuit
Double_t * vstart
Double_t * step
string * parName
Double_t * lowphysbound
Double_t * highphysbound
Double_t fMinLik
Int_t fMinIndexDm2
Int_t fMinIndexSin2
Double_t fMinLikSin2
Double_t fMinLikSin2Error
Double_t fMinLikDm2
Double_t fMinLikDm2Error
TGraph * gBestFit
TFile * fOutputRootFile
vector< GhostSample * > vSamples
TH2D * hOverallSurface
TH2D * hDeltaOverallSurface
TH2D * hCCDetNormalisation
TH2D * hCCDetNormalisationError
TH2D * hCCRockNormalisation
TH2D * hCCRockNormalisationError
TH2D * hNuMuBarNormalisation
TH2D * hNuMuBarNormalisationError
TH2D * hNCNormalisation
TH2D * hNCNormalisationError
TH2D * hDecoherenceSurface
TH2D * hDecoherenceSurfaceError
vector< TH2D * > v2DSystematicSurfaces
vector< TH2D * > v2DSystematicErrorSurfaces
vector< TH2D * > v1DSystematicSurfaces
vector< TH2D * > v1DSystematicErrorSurfaces
Bool_t fPlot
Double_t fTrackGlobal
TH1D * hShwChi2
TH2D * hTrkChi2

Detailed Description

Definition at line 26 of file GhostFramework.h.


Constructor & Destructor Documentation

GhostFramework::~GhostFramework (  ) 

Definition at line 177 of file GhostFramework.cxx.

References f1DInitialSystematicShift, f1DMinErrorSystematicShift, f1DMinSystematicShift, f2DInitialSystematicShift, f2DMinErrorSystematicShift, f2DMinSystematicShift, fInitialScalingFactors, fMinScale, fMinScaleError, Msg::kInfo, and MSG.

00178 {
00179   MSG("GhostFramework",Msg::kInfo) << " *-* GhostFramework::~GhostFramework() *-* " << endl;
00180   delete[] fInitialScalingFactors;
00181   delete[] fMinScale;
00182   delete[] fMinScaleError;
00183   delete[] f2DInitialSystematicShift;
00184   delete[] f2DMinSystematicShift;
00185   delete[] f2DMinErrorSystematicShift;
00186   delete[] f1DInitialSystematicShift;
00187   delete[] f1DMinSystematicShift;
00188   delete[] f1DMinErrorSystematicShift;
00189 }

GhostFramework::GhostFramework ( GhostInputs inputs,
GhostFakeData fake_data 
) [private]

Definition at line 53 of file GhostFramework.cxx.

References f1DInitialSystematicShift, f1DMinErrorSystematicShift, f1DMinSystematicShift, f2DInitialSystematicShift, f2DMinErrorSystematicShift, f2DMinSystematicShift, fCCDetScale, GhostInputs::fCCDetScale, GhostInputs::fCCDetSigma, fCCDetSigma, fCCRockScale, GhostInputs::fCCRockScale, fCCRockSigma, GhostInputs::fCCRockSigma, fDecoPower, GhostInputs::fDecoPower, fDm2Max, GhostInputs::fDm2Max, GhostInputs::fDm2MaxOverall, fDm2MaxOverall, fDm2Min, GhostInputs::fDm2Min, GhostInputs::fDm2MinOverall, fDm2MinOverall, fFitCCDet, GhostInputs::fFitCCDet, GhostInputs::fFitCCRock, fFitCCRock, fFitDecay, GhostInputs::fFitDecay, GhostInputs::fFitDecoherence, fFitDecoherence, fFitNC, GhostInputs::fFitNC, GhostInputs::fFitNuMuBar, fFitNuMuBar, fFitScales, GhostInputs::fFitScales, fInitialScalingFactors, fMinScale, fMinScaleError, GhostInputs::fMinuitInitialValues, fMinuitInitialValues, fMu2_Error, fMu2_High, fMu2_Initial, fMu2_par, GhostInputs::fNCScale, fNCScale, fNCSigma, GhostInputs::fNCSigma, fNuisanceParameters, GhostInputs::fNuisanceParameters, fNum1DSysts, GhostInputs::fNum1DSysts, GhostInputs::fNum2DSysts, fNum2DSysts, fNumDm2Bins, GhostInputs::fNumDm2Bins, GhostInputs::fNumDm2BinsOverall, fNumDm2BinsOverall, GhostInputs::fNumMockDataExpts, fNumMockDataExpts, fNumSamples, GhostInputs::fNumSamples, fNumScalingFactors, GhostInputs::fNumScalingFactors, GhostInputs::fNumSin2Bins, fNumSin2Bins, fNumSin2BinsOverall, GhostInputs::fNumSin2BinsOverall, GhostInputs::fNuMuBarScale, fNuMuBarScale, fNuMuBarSigma, GhostInputs::fRunParallel, fRunParallel, GhostInputs::fSin2Max, fSin2Max, fSin2MaxOverall, GhostInputs::fSin2MaxOverall, GhostInputs::fSin2Min, fSin2Min, fSin2MinOverall, GhostInputs::fSin2MinOverall, GhostInputs::fSurfacePair, fSurfacePair, gFakeData, gInputs, and vSamples.

Referenced by Instance().

00053                                                                            :
00054   fFitDecoherence(false),
00055   fFitDecay(false),
00056   fRunParallel(0),
00057   fFitScales(0),
00058   fMinuitInitialValues(0),
00059   fNumMockDataExpts(1),
00060   fNumSamples(1),
00061   fDm2BinWidth(-999.9),
00062   fSin2BinWidth(-999.9),
00063   fSin2MinOverall(-999.9),
00064   fSin2MaxOverall(-999.9),
00065   fDm2MinOverall(-999.9),
00066   fDm2MaxOverall(-999.9),
00067   fNumDm2BinsOverall(-999),
00068   fNumSin2BinsOverall(-999),
00069   fNum2DSysts(-999),
00070   fNum1DSysts(-999),
00071   fNumSin2Bins(-999),
00072   fNumDm2Bins(-999),
00073   fSin2Min(-999.9),
00074   fSin2Max(-999.9),
00075   fDm2Min(-999.9),
00076   fDm2Max(-999.9),
00077   fNumScalingFactors(0),
00078   fNuisanceParameters(1),
00079   fCCDetSigma(0),
00080   fCCRockSigma(0),
00081   fNCSigma(0),
00082   fNuMuBarSigma(0),  
00083   fCCDetScale(0),
00084   fCCRockScale(0),
00085   fNCScale(0),
00086   fNuMuBarScale(0),
00087   fFitCCDet(0), 
00088   fFitCCRock(0),
00089   fFitNC(0),
00090   fFitNuMuBar(0),
00091   fSin2_par(-999.9),
00092   fDm2_par(-999.9),
00093   fMinLik(-999.9),
00094   fMinIndexDm2(-999),
00095   fMinIndexSin2(-999),
00096   fMinLikSin2(-999.9),
00097   fMinLikSin2Error(-999.9),
00098   fMinLikDm2(-999.9),
00099   fMinLikDm2Error(-999.9),
00100   fPlot(false)
00101 {
00102   gInputs = inputs;
00103   gFakeData = fake_data;
00104 
00105   fFitDecoherence = gInputs->fFitDecoherence;
00106   fFitDecay = gInputs->fFitDecay;
00107   fRunParallel = gInputs->fRunParallel;
00108   fFitScales = gInputs->fFitScales;
00109   fNumMockDataExpts = gInputs->fNumMockDataExpts;
00110   fNumSamples = gInputs->fNumSamples;
00111   fSin2MinOverall = gInputs->fSin2MinOverall;
00112   fSin2MaxOverall = gInputs->fSin2MaxOverall;
00113   fDm2MinOverall = gInputs->fDm2MinOverall;
00114   fDm2MaxOverall = gInputs->fDm2MaxOverall;
00115   fNumDm2BinsOverall = gInputs->fNumDm2BinsOverall;
00116   fNumSin2BinsOverall = gInputs->fNumSin2BinsOverall;
00117   fNum2DSysts = gInputs->fNum2DSysts;
00118   fNum1DSysts = gInputs->fNum1DSysts;
00119   fNumSin2Bins  = gInputs->fNumSin2Bins;
00120   fNumDm2Bins = gInputs->fNumDm2Bins;
00121   fSin2Min = gInputs->fSin2Min;
00122   fSin2Max = gInputs->fSin2Max;
00123   fDm2Min = gInputs->fDm2Min;
00124   fDm2Max = gInputs->fDm2Max;
00125   fNumScalingFactors = gInputs->fNumScalingFactors;
00126   fNuisanceParameters = gInputs->fNuisanceParameters;
00127   fMinuitInitialValues = gInputs->fMinuitInitialValues;
00128 
00129   fCCDetSigma = gInputs->fCCDetSigma;
00130   fCCRockSigma = gInputs->fCCRockSigma;
00131   fNCSigma = gInputs->fNCSigma;
00132   fNuMuBarSigma = gInputs-> fNuMuBarSigma; 
00133   fCCDetScale = gInputs->fCCDetScale;
00134   fCCRockScale = gInputs->fCCRockScale;
00135   fNCScale = gInputs->fNCScale;
00136   fNuMuBarScale = gInputs->fNuMuBarScale;
00137   fFitCCDet = gInputs->fFitCCDet;
00138   fFitCCRock = gInputs->fFitCCRock;
00139   fFitNC = gInputs->fFitNC;
00140   fFitNuMuBar = gInputs->fFitNuMuBar;
00141 
00142   vSamples.clear();
00143 
00144   //Set up arrays for holding fitted systematic shifts
00145   fInitialScalingFactors = new Double_t[fNumScalingFactors];
00146   fMinScale = new Double_t[fNumScalingFactors];
00147   fMinScaleError = new Double_t[fNumScalingFactors];
00148   f2DInitialSystematicShift = new Double_t[fNum2DSysts]; 
00149   f2DMinSystematicShift = new Double_t[fNum2DSysts];
00150   f2DMinErrorSystematicShift = new Double_t[fNum2DSysts];
00151   f1DInitialSystematicShift = new Double_t[fNum1DSysts];
00152   f1DMinSystematicShift = new Double_t[fNum1DSysts];
00153   f1DMinErrorSystematicShift = new Double_t[fNum1DSysts];
00154   for(Int_t i=0;i<fNum2DSysts;i++){
00155     f2DInitialSystematicShift[i] = 0.0;
00156     f2DMinSystematicShift[i] = 0.0;
00157     f2DMinErrorSystematicShift[i] = 0.0;
00158   }
00159   for(Int_t i=0;i<fNum1DSysts;i++){
00160     f1DInitialSystematicShift[i] = 0.0;
00161     f1DMinSystematicShift[i] = 0.0;
00162     f1DMinErrorSystematicShift[i] = 0.0;
00163   }
00164   for(Int_t i=0;i<fNumScalingFactors;i++){
00165     fInitialScalingFactors[i] = 0.0;
00166     fMinScale[i] = 0.0;
00167     fMinScaleError[i] = 0.0;
00168   }
00169   fMu2_Initial = 0.0;
00170   fMu2_High = 1.0;
00171   fMu2_par = 0.0;
00172   fMu2_Error = 0.0;
00173   fSurfacePair = gInputs->fSurfacePair;
00174   fDecoPower = gInputs->fDecoPower;
00175 }


Member Function Documentation

void GhostFramework::AddSample ( GhostSample sample  )  [inline]

Definition at line 37 of file GhostFramework.h.

References vSamples.

00037 { vSamples.push_back(sample); }

void GhostFramework::FindInitialNuisanceValues (  )  [private]

Definition at line 1107 of file GhostFramework.cxx.

References MuELoss::e, f1DInitialSystematicShift, f2DInitialSystematicShift, fCCDetSigma, fCCRockSigma, fcnNuisance(), fFitCCDet, fFitCCRock, fFitDecoherence, fFitNC, fFitNuMuBar, GhostInputs::fGlobalSystematicName, fInitialNMinuit, fInitialScalingFactors, fMu2_High, fMu2_Initial, fNCSigma, fNum1DSysts, fNum2DSysts, fNumScalingFactors, fNuMuBarSigma, gInputs, highphysbound, Msg::kDebug, lowphysbound, MSG, n, parName, step, and vstart.

Referenced by LogLikelihood().

01108 {
01109   TString name;
01110 
01111   MSG("GhostFramework",Msg::kDebug) << " *-* GhostFramework::FindInitialNuisanceValues() *-*" << endl;
01112  
01113   // ============================
01114   // Loop for normalisation systs
01115   // ============================
01116   Int_t ierflg = 0;
01117   Int_t num_par = fNumScalingFactors + fNum2DSysts + fNum1DSysts + fFitDecoherence;
01118   Double_t error = 0.0;
01119   Double_t value = 0.0;
01120   Double_t strategy = 0.;
01121   
01122   TMinuit* temp;
01123   temp = gMinuit;
01124 
01125   Int_t counter = 0;
01126 
01127   for(int i=0;i<num_par;i++){
01128     
01129     value = -999.9; error = -999.9;
01130 
01131     fInitialNMinuit = new TMinuit(num_par);
01132     fInitialNMinuit->SetPrintLevel(-1);
01133     fInitialNMinuit->mnexcm("SET STR",&strategy,1,ierflg);
01134     fInitialNMinuit->SetFCN(fcnNuisance);
01135     
01136     vstart = new Double_t[num_par];  step = new Double_t[num_par];  parName = new string[num_par];
01137     lowphysbound = new Double_t[num_par];  highphysbound = new Double_t[num_par];
01138 
01139     counter = 0;
01140 
01141     //Scaling systematics
01142     if(fFitCCDet){
01143       name.Resize(0); name.Append("cc_det");
01144       parName[counter] = name.Data(); vstart[counter]=0.0; step[counter] = 0.005*fCCDetSigma; lowphysbound[counter] = -1.0; highphysbound[counter] = +2.0; counter++;
01145     }
01146     if(fFitCCRock){
01147       name.Resize(0); name.Append("cc_rock");
01148       parName[counter] = name.Data(); vstart[counter]=0.0; step[counter] = 0.005*fCCRockSigma; lowphysbound[counter] = -1.0; highphysbound[counter] = +2.0; counter++;
01149     }
01150     if(fFitNC){
01151       name.Resize(0); name.Append("nc");
01152       parName[counter] = name.Data(); vstart[counter]=0.0; step[counter] = 0.005*fNCSigma; lowphysbound[counter] = -1.0; highphysbound[counter] = +2.0; counter++;
01153     }
01154     if(fFitNuMuBar){
01155       name.Resize(0); name.Append("numubar");
01156       parName[counter] = name.Data(); vstart[counter]=0.0; step[counter] = 0.005*fNuMuBarSigma; lowphysbound[counter] = -1.0; highphysbound[counter] = +2.0; counter++;
01157     }
01158 
01159     //Energy systematics
01160     for(int j=fNumScalingFactors;j<num_par-fFitDecoherence;j++){
01161       parName[j] = ((gInputs->fGlobalSystematicName).at(j-fNumScalingFactors)).Data(); vstart[j] = 0.0; step[j] = 0.005; lowphysbound[j] = -2.; highphysbound[j] = 2.;
01162     }
01163     if(fFitDecoherence){
01164       parName[num_par-1] = "Mu2"; vstart[num_par-1]= fMu2_Initial; step[num_par-1] = 1e-5; lowphysbound[num_par-1] = 0.0; highphysbound[num_par-1] = fMu2_High; 
01165     }
01166     
01167     for (Int_t n=0 ; n<num_par ; n++) {
01168       fInitialNMinuit->mnparm(n, parName[n].c_str(), vstart[n], step[n], lowphysbound[n], highphysbound[n], ierflg);
01169       fInitialNMinuit->FixParameter(n); //Fix every parameter
01170     }
01171     
01172     fInitialNMinuit->Release(i);  //Release parameter i
01173     
01174     fInitialNMinuit->Migrad();
01175     
01176     if(i<fNumScalingFactors){fInitialNMinuit->GetParameter(i,fInitialScalingFactors[i],error);}
01177     if(i<(fNumScalingFactors + fNum2DSysts) && i>=fNumScalingFactors){fInitialNMinuit->GetParameter(i,f2DInitialSystematicShift[i-fNumScalingFactors],error);}
01178     if(i>=(fNumScalingFactors + fNum2DSysts)&&i<num_par-fFitDecoherence) fInitialNMinuit->GetParameter(i,f1DInitialSystematicShift[i-(fNumScalingFactors + fNum2DSysts)],error);
01179     if(i==num_par-1&&fFitDecoherence) fInitialNMinuit->GetParameter(i,fMu2_Initial,error);
01180 
01181     fInitialNMinuit->Clear();
01182     delete fInitialNMinuit;
01183     delete[] vstart;
01184     delete[] step;
01185     delete[] parName;
01186     delete lowphysbound;
01187     delete highphysbound;
01188   }
01189   gMinuit = temp;
01190 }

GhostFramework * GhostFramework::Instance ( GhostInputs inputs,
GhostFakeData fake_data 
) [static]

Definition at line 45 of file GhostFramework.cxx.

References fGhostFramework, and GhostFramework().

00046 { 
00047   if(!fGhostFramework){
00048     fGhostFramework = new GhostFramework(inputs,fake_data);
00049   }
00050   return fGhostFramework;
00051 }

Double_t GhostFramework::LogLikelihood ( Double_t  sin2_par,
Double_t  dm2_par 
)

Definition at line 627 of file GhostFramework.cxx.

References MuELoss::e, f1DInitialSystematicShift, f1DMinErrorSystematicShift, f1DMinSystematicShift, GhostInputs::f1DSystematicName, f2DInitialSystematicShift, f2DMinErrorSystematicShift, f2DMinSystematicShift, GhostInputs::f2DSystematicName, fCCDetSigma, fCCRockSigma, fcnNuisance(), fDm2_par, fFitCCDet, fFitCCRock, fFitDecoherence, fFitNC, fFitNuMuBar, fFreeMinuit, GhostInputs::fGlobalSystematicName, FindInitialNuisanceValues(), fInitialScalingFactors, fMinScale, fMinScaleError, fMinuitInitialValues, fMu2_Error, fMu2_High, fMu2_Initial, fMu2_par, fNCSigma, fNum1DSysts, fNum2DSysts, fNumSamples, fNumScalingFactors, fNuMuBarSigma, fPlot, fResolutionLikelihood, fSin2_par, fSurfacePair, fTrackGlobal, gInputs, hCCDetNormalisation, hCCDetNormalisationError, hCCRockNormalisation, hCCRockNormalisationError, hDecoherenceSurface, hDecoherenceSurfaceError, highphysbound, hNCNormalisation, hNCNormalisationError, hNuMuBarNormalisation, hNuMuBarNormalisationError, hOverallSurface, Msg::kDebug, Msg::kInfo, lowphysbound, MSG, n, parName, size, step, v1DSystematicErrorSurfaces, v1DSystematicSurfaces, v2DSystematicErrorSurfaces, v2DSystematicSurfaces, GhostInputs::vFinalData3D, vSamples, and vstart.

Referenced by OscFitGridSearch().

00628 {
00629   MSG("GhostFramework",Msg::kDebug) << " *-* GhostFramework::LogLikelihood() *-* " << endl;
00630   
00631   Double_t lnL = 0.;
00632 
00633   fSin2_par = sin2_par;
00634   fDm2_par = dm2_par;
00635 
00636   // ===============
00637   // Systematics Fit
00638   // ===============
00639   if((fNumScalingFactors + fNum2DSysts + fNum1DSysts)>0||fFitDecoherence){
00640 
00641     Int_t counter = 0;
00642 
00643     for(Int_t i=0;i<fNumScalingFactors;i++){
00644       fInitialScalingFactors[i] = 0.0;
00645       fMinScale[i] = 0.0;
00646       fMinScaleError[i] = 0.0;
00647     }
00648     for(Int_t i=0;i<fNum2DSysts;i++){
00649       f2DInitialSystematicShift[i] = 0.0;
00650       f2DMinSystematicShift[i] = 0.0;
00651       f2DMinErrorSystematicShift[i] = 0.0;
00652     }
00653     for(Int_t i=0;i<fNum1DSysts;i++){
00654       f1DInitialSystematicShift[i] = 0.0;
00655       f1DMinSystematicShift[i] = 0.0;
00656       f1DMinErrorSystematicShift[i] = 0.0;
00657     }
00658 
00659     switch(fSurfacePair){
00660       case 1: fMu2_Initial = 2.4e-3; fMu2_High = 1.0; break;
00661       case 2: fMu2_Initial = 0.95; fMu2_High = 1.0; break;
00662       default: fMu2_Initial = 0.0; fMu2_High = 1.0;
00663     }
00664     fMu2_par = 0.0;
00665     fMu2_Error = 0.0;
00666 
00667     TMinuit* old = gMinuit;
00668     if(fMinuitInitialValues) this->FindInitialNuisanceValues();
00669     gMinuit = old;
00670 
00671 
00672     //Run main "free" Minuit
00673     Int_t ierflg = 0;
00674     Int_t num_par = fNumScalingFactors + fNum2DSysts + fNum1DSysts + fFitDecoherence;
00675     Double_t strategy = 0.;
00676 
00677     fFreeMinuit = new TMinuit(num_par);
00678     fFreeMinuit->SetPrintLevel(-1);
00679     fFreeMinuit->SetFCN(fcnNuisance);
00680     fFreeMinuit->mnexcm("SET STR",&strategy,1,ierflg);
00681     
00682     //Set limits on parameters
00683     vstart = new Double_t[num_par];  step = new Double_t[num_par];  parName = new string[num_par];
00684     lowphysbound = new Double_t[num_par];  highphysbound = new Double_t[num_par];
00685 
00686     TString name;
00687 
00688     if(fFitCCDet){
00689       name.Resize(0); name.Append("cc_det");
00690       parName[counter] = name.Data(); vstart[counter] = fInitialScalingFactors[counter]; step[counter] = 0.005*fCCDetSigma; lowphysbound[counter] = -1.0; highphysbound[counter] = +2.0; counter++;
00691     }
00692     if(fFitCCRock){
00693       name.Resize(0); name.Append("cc_rock");
00694       parName[counter] = name.Data(); vstart[counter] = fInitialScalingFactors[counter]; step[counter] = 0.005*fCCRockSigma; lowphysbound[counter] = -1.0; highphysbound[counter] = +2.0; counter++;
00695     }
00696     if(fFitNC){
00697       name.Resize(0); name.Append("nc");
00698       parName[counter] = name.Data(); vstart[counter] = fInitialScalingFactors[counter]; step[counter] = 0.005*fNCSigma; lowphysbound[counter] = -1.0; highphysbound[counter] = +2.0; counter++;
00699     }
00700     if(fFitNuMuBar){
00701       name.Resize(0); name.Append("numubar");
00702       parName[counter] = name.Data(); vstart[counter] = fInitialScalingFactors[counter]; step[counter] = 0.005*fNuMuBarSigma; lowphysbound[counter] = -1.0; highphysbound[counter] = +2.0; counter++;
00703     }
00704 
00705     counter = 0;
00706     
00707     //2D Energy systematics
00708     for(int j=(fNumScalingFactors);j<(fNumScalingFactors + fNum2DSysts);j++){
00709       parName[j] = ((gInputs->fGlobalSystematicName).at(j-fNumScalingFactors)).Data(); vstart[j] = f2DInitialSystematicShift[j-fNumScalingFactors]; step[j] = 0.005; lowphysbound[j] = -2.; highphysbound[j] = 2.;
00710     }
00711     //1D Energy systematics    
00712     for(int j=(fNumScalingFactors + fNum2DSysts);j<(fNumScalingFactors + fNum2DSysts + fNum1DSysts);j++){
00713       parName[j] = ((gInputs->fGlobalSystematicName).at(j-fNumScalingFactors)).Data(); vstart[j] = f1DInitialSystematicShift[j-(fNumScalingFactors + fNum2DSysts)]; step[j] = 0.005; lowphysbound[j] = -2.; highphysbound[j] = 2.;
00714     }
00715 
00716     if(fFitDecoherence){
00717       parName[num_par-1] = "Mu2"; vstart[num_par-1] = fMu2_Initial; step[num_par-1] = 1e-5; lowphysbound[num_par-1] = 0.0; highphysbound[num_par-1] = fMu2_High;
00718     }
00719     
00720     if(fPlot){
00721       vstart[3] = fTrackGlobal;
00722     }
00723 
00724     for (Int_t n=0 ; n<num_par ; n++) {
00725       fFreeMinuit->mnparm(n, parName[n].c_str(), vstart[n], step[n], lowphysbound[n], highphysbound[n], ierflg);
00726     }
00727     
00728     if(fPlot){
00729       //Fix track parameter to global
00730       fFreeMinuit->FixParameter(3);
00731     }
00732     
00733     fFreeMinuit->Migrad();
00734     MSG("GhostInputs",Msg::kDebug) << " *-* GhostInputs::GetLikelihood() Got Likelihood Value *-*" << endl;
00735     
00736     //Get parameters
00737     for(Int_t i=0;i<num_par-fFitDecoherence;i++){
00738       if(i<fNumScalingFactors) fFreeMinuit->GetParameter(i,fMinScale[i],fMinScaleError[i]);
00739       if(i>=fNumScalingFactors && i<(fNumScalingFactors + fNum2DSysts)) fFreeMinuit->GetParameter(i,f2DMinSystematicShift[i-fNumScalingFactors],f2DMinErrorSystematicShift[i-fNumScalingFactors]);
00740       if(i>=(fNumScalingFactors + fNum2DSysts)) fFreeMinuit->GetParameter(i,f1DMinSystematicShift[i-(fNumScalingFactors+fNum2DSysts)],f1DMinErrorSystematicShift[i-(fNumScalingFactors+fNum2DSysts)]);
00741     }
00742     if(fFitDecoherence) fFreeMinuit->GetParameter(num_par-1,fMu2_par,fMu2_Error); 
00743     
00744     lnL+=fFreeMinuit->fAmin;
00745     
00746     
00747     //-----------------------------------------------------------------------------------------------------------------------------------------------------
00748     
00749     //Fill appropriate histograms for this point in parameter space
00750     hOverallSurface->Fill(sin2_par,dm2_par,lnL);
00751     
00752     //(1) Scaling; loop over samples and fill accordingly based on whether scaling factor exists
00753     Int_t arraycounter = 0;
00754     counter = 0;
00755 
00756     if(fFitCCDet){
00757       hCCDetNormalisation->Fill(sin2_par,dm2_par,fMinScale[arraycounter]/fCCDetSigma);
00758       hCCDetNormalisationError->Fill(sin2_par,dm2_par,fMinScaleError[arraycounter]/fCCDetSigma);
00759       counter++; arraycounter++;}
00760     if(fFitCCRock){
00761       hCCRockNormalisation->Fill(sin2_par,dm2_par,fMinScale[arraycounter]/fCCRockSigma);
00762       hCCRockNormalisationError->Fill(sin2_par,dm2_par,fMinScaleError[arraycounter]/fCCRockSigma);
00763       counter++; arraycounter++;}
00764     if(fFitNC){
00765       hNCNormalisation->Fill(sin2_par,dm2_par,fMinScale[arraycounter]/fNCSigma);
00766       hNCNormalisationError->Fill(sin2_par,dm2_par,fMinScaleError[arraycounter]/fNCSigma);
00767       counter++; arraycounter++;}
00768     if(fFitNuMuBar){
00769       hNuMuBarNormalisation->Fill(sin2_par,dm2_par,fMinScale[arraycounter]/fNuMuBarSigma);
00770       hNuMuBarNormalisationError->Fill(sin2_par,dm2_par,fMinScaleError[arraycounter]/fNuMuBarSigma);
00771       counter++; arraycounter++;}
00772     
00773     counter = 0;
00774     
00775     //(2) Energy systematics
00776     for(Int_t i=0;i<fNum2DSysts;i++){
00777       v2DSystematicSurfaces.at(i)->Fill(sin2_par,dm2_par,f2DMinSystematicShift[i]);
00778       v2DSystematicErrorSurfaces.at(i)->Fill(sin2_par,dm2_par,f2DMinErrorSystematicShift[i]);
00779     }
00780     for(Int_t i=0;i<fNum1DSysts;i++){
00781       v1DSystematicSurfaces.at(i)->Fill(sin2_par,dm2_par,f1DMinSystematicShift[i]);
00782       v1DSystematicErrorSurfaces.at(i)->Fill(sin2_par,dm2_par,f1DMinErrorSystematicShift[i]);
00783     }
00784 
00785     if(fFitDecoherence){
00786       hDecoherenceSurface->Fill(sin2_par,dm2_par,fMu2_par);
00787       hDecoherenceSurfaceError->Fill(sin2_par,dm2_par,fMu2_Error);
00788     }
00789     
00790     //-----------------------------------------------------------------------------------------------------------------------------------------------------
00791     
00792     
00793     //Format MSG output for this function
00794     //MsgFormat ifmt("%8.6f");
00795     
00796     cout << " sstt " << sin2_par << " dm2 " <<  dm2_par << "; ll = " << lnL << endl;
00797     if(fFitCCDet){cout << "cc_detector_scale (/sigma) = " << fMinScale[counter]/fCCDetSigma << endl; counter++;}
00798     if(fFitCCRock){cout << "cc_rock_scale (/sigma) = " << fMinScale[counter]/fCCRockSigma << endl; counter++;}
00799     if(fFitNC){cout << "nc_scale (/sigma) = " << fMinScale[counter]/fNCSigma << endl; counter++;}
00800     if(fFitNuMuBar){cout << "ws_scale (/sigma) = " << fMinScale[counter]/fNuMuBarSigma << endl; counter++;}
00801     
00802     for(Int_t i=0;i<fNum2DSysts;i++){
00803       cout << ((gInputs->f2DSystematicName).at(i)).Data() << "(/sigma) = " << f2DMinSystematicShift[i] << endl;
00804     }
00805     for(Int_t i=0;i<fNum1DSysts;i++){
00806       cout << ((gInputs->f1DSystematicName).at(i)).Data() << "(/sigma) = " << f1DMinSystematicShift[i] << endl;
00807     }
00808 
00809     if(fFitDecoherence) cout << "mu2 = " << fMu2_par << " eV^{2}" << endl;
00810     
00811     fFreeMinuit->Clear();
00812     delete fFreeMinuit;
00813     delete[] vstart;
00814     delete[] step;
00815     delete[] parName;
00816     delete[] lowphysbound;
00817     delete[] highphysbound;
00818     
00819     gMinuit = old; //Reset global minuit
00820     
00821   }//EO fit with systematics
00822   
00823   // ===============
00824   // Statistical Fit
00825   // ===============
00826   else if((fNumScalingFactors + fNum2DSysts + fNum1DSysts)==0&&!fFitDecoherence){
00827 
00828     MSG("GhostFramework",Msg::kDebug) << "Statistical Fit" << endl;
00829 
00830     //Format MSG output for this function
00831     //MsgFormat ifmt("%8.6f");
00832 
00833     //LnL calculation
00834     //Only need to consider central values
00835     Double_t temp_lik = 0.0;
00836   
00837     //Loop over samples
00838     for(Int_t i=0;i<fNumSamples;i++){
00839 
00840       vector<vector<ArrayTH1D*> > vStatPredicted2D;
00841       vector<ArrayTH1D*> vStatPredicted1D;
00842       for(Int_t j=0;j<vSamples.at(i)->fNDis;j++) vStatPredicted2D.push_back(vStatPredicted1D);
00843       
00844       Double_t* syst2d = 0; Double_t* syst1d = 0;
00845     
00846       vSamples.at(i)->GetSpectrum(fSin2_par,fDm2_par,0,0,0,0,syst2d,syst1d,vStatPredicted2D,false,0);
00847       
00848       //Loop over distributions for this sample and determine likelihood
00849       for(Int_t j=0;j<(vSamples.at(i))->fNDis;j++){
00850         
00851         //Choose which form of likelihood function to use: resolution shape or standard
00852         //(1) RESOLUTION LIKELIHOOD
00853         if(fResolutionLikelihood){
00854           
00855           //Doubles to hold OVERALL NORMALISATION TERM for this distribution (run) 
00856           Double_t PredNormTotal = 0.;
00857           Double_t DataNormTotal = 0.;
00858           
00859           //Loop over resolution bins
00860           for(Int_t res=0;res<(vSamples.at(i))->fNResBins;res++){
00861 
00862             TH1D* hPredicted = (vStatPredicted2D.at(j)).at(res)->ConvertTH1D("hPredicted");
00863             Double_t PredIntegral = hPredicted->Integral();
00864             TH1D* hData = (((gInputs->vFinalData3D.at(i)).at(j)).at(res))->ConvertTH1D("hData");
00865             Double_t DataIntegral = hData->Integral();
00866             PredNormTotal += PredIntegral;
00867             DataNormTotal += DataIntegral;
00868             
00869             //Loop over energy for this resolution bin
00870             Double_t PredBinContent = -1.0;
00871             Double_t DataBinContent = -1.0;
00872             
00873             for(Int_t bin=0; bin<(vSamples.at(i))->fNumFittingBins;bin++){
00874               
00875               if(PredIntegral>0) PredBinContent = (vStatPredicted2D.at(j)).at(res)->GetBinContent(bin+1)/PredIntegral;
00876               else cout << "Error: Empty quantile!" << endl;
00877               DataBinContent = (((gInputs->vFinalData3D.at(i)).at(j)).at(res))->GetBinContent(bin+1);
00878               
00879               if(PredBinContent>0){
00880                 temp_lik -= DataBinContent*log(PredBinContent);
00881               }
00882               else if(PredBinContent<0){
00883                 temp_lik += 1.0e10;
00884               }
00885               if(DataBinContent>0){
00886                 temp_lik += DataBinContent*log(DataBinContent) - DataBinContent;
00887               }
00888               
00889               
00890             }//EO Energy bins loop
00891             delete hPredicted;
00892             delete hData;
00893           }//EO Res Bins loop
00894           
00895           if(PredNormTotal>0) temp_lik += PredNormTotal - DataNormTotal*log(PredNormTotal);  //overall normalisation term
00896           
00897         }
00898  
00899         
00900         //(2) STANDARD LIKELIHOOD
00901         else{
00902 
00903           Double_t FD_pred_bin_content = -1.0;
00904           Double_t data_bin_content = -1.0; 
00905           
00906           for(Int_t res=0;res<(vSamples.at(i))->fNResBins;res++){
00907             for(Int_t bin=0; bin<(vSamples.at(i))->fNumFittingBins;bin++){
00908 
00909               FD_pred_bin_content = (vStatPredicted2D.at(j)).at(res)->GetBinContent(bin+1);
00910               data_bin_content = ((gInputs->vFinalData3D.at(i)).at(j)).at(res)->GetBinContent(bin+1);
00911 
00912               Double_t bin_lik = 0.0;
00913                       
00914               if(FD_pred_bin_content<1e-7) FD_pred_bin_content = 1e-7;
00915               if( FD_pred_bin_content>0.0 ){
00916                 temp_lik += 2*(FD_pred_bin_content-data_bin_content*log(FD_pred_bin_content));
00917                 bin_lik += 2*(FD_pred_bin_content-data_bin_content*log(FD_pred_bin_content));
00918                 if( data_bin_content>0.0 ){
00919                   temp_lik -= 2*(data_bin_content-data_bin_content*log(data_bin_content));
00920                   bin_lik -= 2*(data_bin_content-data_bin_content*log(data_bin_content));
00921                 }
00922               }
00923               if(FD_pred_bin_content==0){cout << "ERRORRRRRRRR" << endl;}
00924             }//EO loop over bins
00925           }//EO loop over resolution bins
00926           
00927         }//EO standard likelihood function
00928       }//EO loop over distributions
00929       
00930       for(UInt_t j=0;j<vStatPredicted2D.size();j++){
00931         for(UInt_t k=0;k<(vStatPredicted2D.at(j)).size();k++){
00932           delete (vStatPredicted2D.at(j)).at(k);
00933         }
00934       }
00935     }//EO loop over samples
00936     
00937     MSG("GhostFramework",Msg::kInfo) << " sstt " << sin2_par << " dm2 " <<  dm2_par << " ll " << temp_lik << endl;
00938     
00939     lnL = temp_lik;
00940     
00941     
00942     //-----------------------------------------------------------------------------------------------------------------------------------------------------
00943     //Fill appropriate histograms for this point in parameter space
00944     hOverallSurface->Fill(sin2_par,dm2_par,lnL);
00945     //-----------------------------------------------------------------------------------------------------------------------------------------------------
00946   }
00947   else{cout << "ERROR in setting number of systematics, exiting" << endl;return -1;}
00948   
00949   return lnL;
00950 }

Double_t GhostFramework::LogLikelihoodNuisance ( Double_t *  par  ) 

Definition at line 953 of file GhostFramework.cxx.

References MuELoss::e, fCCDetSigma, fCCRockSigma, fDm2_par, fFitCCDet, fFitCCRock, fFitDecoherence, fFitNC, fFitNuMuBar, fNCSigma, fNuisanceParameters, fNum1DSysts, fNum2DSysts, fNumSamples, fNumScalingFactors, fNuMuBarSigma, fResolutionLikelihood, fSin2_par, fSurfacePair, gInputs, Msg::kDebug, MSG, size, GhostInputs::vFinalData3D, and vSamples.

00954 {
00955 
00956   MSG("GhostFramework",Msg::kDebug) << " *-* GhostFramework::LogLikelihoodNuisance() *-* " << endl;
00957   Double_t temp_lik = 0.0;
00958 
00959   //Energy systematics
00960   Double_t* syst2d = 0; if(fNum2DSysts>0) syst2d = new Double_t[fNum2DSysts];
00961   Double_t* syst1d = 0; if(fNum1DSysts>0) syst1d = new Double_t[fNum1DSysts];
00962   
00963   for(Int_t j=0;j<fNum2DSysts;j++) syst2d[j] = par[fNumScalingFactors + j];
00964   for(Int_t j=0;j<fNum1DSysts;j++) syst1d[j] = par[fNumScalingFactors + fNum2DSysts +j];
00965   
00966   Double_t mu2_shift = 0.;
00967   if(fFitDecoherence){mu2_shift = par[fNumScalingFactors + fNum2DSysts + fNum1DSysts];}
00968   
00969   //Scaling systematics
00970   Int_t counter = 0;
00971 
00972   Double_t cc_detshift = 0.; Double_t cc_rockshift = 0.;
00973   Double_t nc_shift = 0.; Double_t numubar_shift = 0.;
00974   if(fFitCCDet){cc_detshift = par[counter]; counter++;}
00975   if(fFitCCRock){cc_rockshift = par[counter]; counter++;}
00976   if(fFitNC){nc_shift = par[counter]; counter++;}
00977   if(fFitNuMuBar){numubar_shift = par[counter]; counter++;}
00978   
00979   for(Int_t i=0;i<fNumSamples;i++){
00980     
00981     vector<vector<ArrayTH1D*> > vStatPredicted2D;
00982     vector<ArrayTH1D*> vStatPredicted1D;
00983     for(Int_t j=0;j<vSamples.at(i)->fNDis;j++) vStatPredicted2D.push_back(vStatPredicted1D);
00984     
00985     switch(fSurfacePair){
00986       case 1: 
00987       vSamples.at(i)->GetSpectrum(fSin2_par,mu2_shift,cc_detshift,cc_rockshift,nc_shift,numubar_shift,syst2d,syst1d,vStatPredicted2D,false,fDm2_par);
00988       break;
00989       case 2: 
00990       vSamples.at(i)->GetSpectrum(mu2_shift,fDm2_par,cc_detshift,cc_rockshift,nc_shift,numubar_shift,syst2d,syst1d,vStatPredicted2D,false,fSin2_par);
00991       break;
00992       default: 
00993       vSamples.at(i)->GetSpectrum(fSin2_par,fDm2_par,cc_detshift,cc_rockshift,nc_shift,numubar_shift,syst2d,syst1d,vStatPredicted2D,false,mu2_shift);
00994     }
00995 
00996     //Loop over distributions for this sample and determine likelihood
00997     for(Int_t j=0;j<(vSamples.at(i))->fNDis;j++){
00998       
00999       //Choose which form of likelihood function to use: resolution shape or standard
01000       //(1) RESOLUTION LIKELIHOOD
01001       if(fResolutionLikelihood){
01002         
01003         //Doubles to hold OVERALL NORMALISATION TERM for this distribution (run) 
01004         Double_t PredNormTotal = 0.;
01005         Double_t DataNormTotal = 0.;
01006         
01007         //Loop over resolution bins
01008         for(Int_t res=0;res<(vSamples.at(i))->fNResBins;res++){
01009           
01010           TH1D* hPredicted = (vStatPredicted2D.at(j)).at(res)->ConvertTH1D("hPredicted");
01011           Double_t PredIntegral = hPredicted->Integral();
01012           TH1D* hData = (((gInputs->vFinalData3D.at(i)).at(j)).at(res))->ConvertTH1D("hData");
01013           Double_t DataIntegral = hData->Integral();
01014           PredNormTotal += PredIntegral;
01015           DataNormTotal += DataIntegral;
01016             
01017           //Loop over energy for this resolution bin
01018           Double_t PredBinContent = -1.0;
01019           Double_t DataBinContent = -1.0;
01020           
01021           for(Int_t bin=0; bin<(vSamples.at(i))->fNumFittingBins;bin++){
01022             
01023             if(PredIntegral>0) PredBinContent = (vStatPredicted2D.at(j)).at(res)->GetBinContent(bin+1)/PredIntegral;
01024             else cout << "Error: Empty quantile!" << endl;
01025             DataBinContent = (((gInputs->vFinalData3D.at(i)).at(j)).at(res))->GetBinContent(bin+1);
01026             
01027             if(PredBinContent>0){
01028               temp_lik -= DataBinContent*log(PredBinContent);
01029             }
01030             else if(PredBinContent<0){
01031                 temp_lik += 1.0e10;
01032             }
01033             if(DataBinContent>0){
01034               temp_lik += DataBinContent*log(DataBinContent) - DataBinContent;
01035             }
01036             
01037             
01038           }//EO Energy bins loop
01039           delete hPredicted;
01040           delete hData;
01041         }//EO Res Bins loop
01042           
01043         if(PredNormTotal>0) temp_lik += PredNormTotal - DataNormTotal*log(PredNormTotal);  //overall normalisation term
01044         
01045       }
01046       
01047       
01048       //(2) STANDARD LIKELIHOOD
01049       else{
01050           
01051         Double_t FD_pred_bin_content = -1.0;
01052         Double_t data_bin_content = -1.0; 
01053         
01054         for(Int_t res=0;res<(vSamples.at(i))->fNResBins;res++){
01055           for(Int_t bin=0; bin<(vSamples.at(i))->fNumFittingBins;bin++){
01056             
01057             FD_pred_bin_content = (vStatPredicted2D.at(j)).at(res)->GetBinContent(bin+1);
01058             data_bin_content = ((gInputs->vFinalData3D.at(i)).at(j)).at(res)->GetBinContent(bin+1);
01059               
01060             //cout << "FD_pred_bin_content = " << FD_pred_bin_content << endl;
01061 
01062             if(FD_pred_bin_content<1e-7) FD_pred_bin_content = 1e-7;
01063 
01064             //cout << "FD_pred_bin_content after checking is too small = " << FD_pred_bin_content << endl;
01065 
01066               if( FD_pred_bin_content>0.0 ){
01067                 temp_lik += 2*(FD_pred_bin_content-data_bin_content*log(FD_pred_bin_content));
01068                 if( data_bin_content>0.0 ){
01069                   temp_lik -= 2*(data_bin_content-data_bin_content*log(data_bin_content));
01070                 }
01071               }
01072               if(FD_pred_bin_content==0){cout << "ERRORRRRRRRR" << endl;}
01073           }//EO loop over bins
01074         }//EO loop over resolution bins
01075         
01076       }//EO standard likelihood function
01077     }//EO loop over distributions
01078     
01079     for(UInt_t j=0;j<vStatPredicted2D.size();j++){
01080       for(UInt_t k=0;k<(vStatPredicted2D.at(j)).size();k++){
01081         delete (vStatPredicted2D.at(j)).at(k);
01082       }
01083     }
01084   }//EO loop over samples
01085   counter = 0;
01086   
01087   if(fNuisanceParameters){
01088     temp_lik += cc_detshift*cc_detshift/(fCCDetSigma*fCCDetSigma);
01089     temp_lik += cc_rockshift*cc_rockshift/(fCCRockSigma*fCCRockSigma);
01090     temp_lik += nc_shift*nc_shift/(fNCSigma*fNCSigma);
01091     temp_lik += numubar_shift*numubar_shift/(fNuMuBarSigma*fNuMuBarSigma);
01092     
01093     //2D systematics  
01094     for(Int_t j=0;j<fNum2DSysts;j++) temp_lik += syst2d[j]*syst2d[j];
01095     
01096     //1D systematics
01097     for(Int_t j=0;j<fNum1DSysts;j++) temp_lik += syst1d[j]*syst1d[j];
01098   }
01099   
01100   if(fNum2DSysts>0) delete[] syst2d;
01101   if(fNum1DSysts>0) delete[] syst1d;
01102   return temp_lik;
01103 }

void GhostFramework::OscFitGridSearch (  )  [private]

Definition at line 559 of file GhostFramework.cxx.

References fDm2BinWidth, fDm2Max, fDm2Min, fNumDm2Bins, fNumSin2Bins, fPlot, fSin2BinWidth, fSin2Max, fSin2Min, fTrackGlobal, hTrkChi2, Msg::kDebug, Msg::kInfo, LogLikelihood(), and MSG.

Referenced by RunFit().

00560 {
00561   MSG("GhostFramework",Msg::kInfo) << " *-* GhostFramework::OscFitGridSearch() *-* " << endl;
00562 
00563   //Double_t fTrackGlobalArray[] = {-1.66,-1.5,-1.4,-1.3,-1.2,-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.45,0.0};
00564   //hTrkChi2 = new TH2D("hTrkChi2","hTrkChi2",50,2e-3,2.8e-3,14,fTrackGlobalArray);
00565   
00566   
00567   //Loop over (this section of) parameter space
00568   for(Int_t sin=0;sin<fNumSin2Bins;sin++){
00569    
00570     Double_t dm2_temp = -1.0;
00571     Double_t sin2_temp = -1.0;
00572 
00573     sin2_temp = fSin2Min + (sin+0.50)*fSin2BinWidth;
00574     
00575     if(sin2_temp<fSin2Max){
00576 
00577       for(Int_t j=0;j<fNumDm2Bins;j++){
00578         
00579         dm2_temp = fDm2Min + (j+0.50)*fDm2BinWidth;
00580 
00581         if(dm2_temp<fDm2Max){
00582 
00583           //Calculate log likelihood      
00584           Double_t lnL;
00585           
00586           //If set plotting flag, set global variable for energy shift, then fix in likelihood function
00587           if(fPlot){
00588 
00589             //Loop over values
00590             //fTrackGlobal = fTrackGlobalArray[i];
00591             //cout << "fTrackGlobal = " << fTrackGlobal << endl;
00592             
00593             lnL = this->LogLikelihood(sin2_temp,dm2_temp);
00594             
00595             hTrkChi2->Fill(dm2_temp,fTrackGlobal,lnL);
00596           }
00597           if(!fPlot) lnL = this->LogLikelihood(sin2_temp,dm2_temp);
00598           
00599           //TCanvas* ctrk = new TCanvas();
00600           //ctrk->cd();
00601           //hTrkChi2->Draw();
00602           //this->Pause();
00603           
00604           MSG("GhostFramework",Msg::kDebug) << " *-* OscFitFramework::OscFitGridSearch() lnL = " << lnL << " *-* " << endl;
00605           
00606           Double_t percent_done = 0.;
00607           percent_done = 100.0000*(fNumDm2Bins*sin + j+1)/(fNumSin2Bins*fNumDm2Bins);
00608           MSG("GhostFramework",Msg::kInfo) << " *-* " << percent_done << "% completed *-* "  << endl;
00609           
00610         }//EO check on dm2 value
00611       }//EO loop over dm2
00612     }//EO check on sin2
00613   }//EO loop over sin2
00614   
00615   
00616   if(fPlot){
00617     TFile *f = new TFile("TrackShiftTestingPlots.root","update");
00618     f->cd();
00619     hTrkChi2->Write();
00620     f->Close();
00621     delete f;
00622   }
00623 }

void GhostFramework::Pause (  ) 

Definition at line 1218 of file GhostFramework.cxx.

References gSystem().

01219 {
01220   cout << "Press return to continue" << endl;
01221 
01222   int key = 0;
01223   while(1){
01224     gSystem->ProcessEvents();
01225     fcntl(STDIN_FILENO, F_SETFL, O_NONBLOCK);
01226     key = getchar();
01227     if(key=='\n') break;
01228     usleep(1000);
01229   }
01230 }

void GhostFramework::RebinRealData (  )  [private]

Definition at line 1193 of file GhostFramework.cxx.

References fNumSamples, gInputs, Msg::kDebug, MSG, size, GhostInputs::vFinalData3D, GhostInputs::vInputData3D, and vSamples.

Referenced by RunFit().

01194 {
01195   MSG("GhostFramework",Msg::kDebug) << " *-* GhostFramework::RebinRealData() *-*" << endl;
01196 
01197   for(Int_t i=0;i<fNumSamples;i++){ //Samples
01198 
01199     for(UInt_t j=0;j<((gInputs->vInputData3D).at(i)).size();j++){  //Distributions
01200       
01201       for(UInt_t k=0;k<((gInputs->vInputData3D.at(i)).at(j)).size();k++){  //Resolution bins
01202         
01203         TString name("hFinalData_"); name += i; name.Append("_"); name += j; name.Append("_"); name += k; 
01204         ((gInputs->vFinalData3D.at(i)).at(j)).push_back(new ArrayTH1D(name.Data(),name.Data(),(vSamples.at(i))->fNumFittingBins,0,(vSamples.at(i))->fMaxEnergy,(vSamples.at(i))->fArrayHistsFlag));
01205 
01206         for(Int_t bin=0;bin<(((gInputs->vInputData3D.at(i)).at(j)).at(k))->GetXaxis()->GetNbins();bin++){
01207           Double_t bin_content = ((gInputs->vInputData3D.at(i)).at(j)).at(k)->GetBinContent(bin+1);
01208           Double_t energy = ((gInputs->vInputData3D.at(i)).at(j)).at(k)->GetXaxis()->GetBinCenter(bin+1);
01209           ((gInputs->vFinalData3D.at(i)).at(j)).at(k)->Fill(energy,bin_content);
01210         }
01211       }
01212     }
01213   }
01214   MSG("GhostFramework",Msg::kDebug) << " *-* Rebinned Real Data Spectra *-* " << endl;
01215 }

void GhostFramework::RunFit (  ) 

Definition at line 194 of file GhostFramework.cxx.

References GhostInputs::f1DSystematicName, GhostInputs::f2DSystematicName, GhostInputs::fBestFitRootFilename, GhostInputs::fBestFitTextFilename, fCCDetSigma, fCCRockSigma, fDm2BinWidth, fDm2MaxOverall, fDm2MinOverall, fFitCCDet, fFitCCRock, fFitDecoherence, fFitNC, fFitNuMuBar, GhostInputs::fMCData, fMinIndexDm2, fMinIndexSin2, fMinLik, fMinLikDm2, fMinLikDm2Error, fMinLikSin2, fMinLikSin2Error, fNCSigma, fNum1DSysts, fNum2DSysts, fNumDm2Bins, fNumDm2BinsOverall, fNumMockDataExpts, fNumSamples, fNumSin2Bins, fNumSin2BinsOverall, fNuMuBarSigma, fOutputRootFile, GhostInputs::fRealData, fRunParallel, fSin2BinWidth, fSin2MaxOverall, fSin2MinOverall, gBestFit, GhostFakeData::GenerateSpectra(), gFakeData, gInputs, hCCDetNormalisation, hCCDetNormalisationError, hCCRockNormalisation, hCCRockNormalisationError, hDecoherenceSurface, hDecoherenceSurfaceError, hDeltaOverallSurface, hNCNormalisation, hNCNormalisationError, hNuMuBarNormalisation, hNuMuBarNormalisationError, hOverallSurface, Msg::kDebug, Msg::kInfo, MSG, OscFitGridSearch(), RebinRealData(), size, v1DSystematicErrorSurfaces, v1DSystematicSurfaces, v2DSystematicErrorSurfaces, v2DSystematicSurfaces, GhostInputs::vFinalData3D, and vSamples.

00195 {
00196   MSG("GhostFramework",Msg::kInfo) << " *-* GhostFramework::RunFit() *-* " << endl;
00197 
00198   gStyle->SetOptStat(0);
00199 
00200 
00201   //Rebin real data if using here, with binning schemes from each individual sample
00202   if(gInputs->fRealData){
00203     this->RebinRealData();
00204     MSG("GhostFramework",Msg::kInfo) << " *-* GhostFramework::RunFit() Rebinned Real Data*-* " << endl;
00205   }
00206   
00207   //LOOP OVER NUMBER OF MOCK DATA EXPERIMENTS
00208   for(int nexp=0;nexp<fNumMockDataExpts;nexp++){
00209 
00210     MSG("GhostFramework",Msg::kInfo) << " " << endl;
00211     MSG("GhostFramework",Msg::kInfo) << "\033[22;31m" << "*-* GhostFramework::RunFit() Experiment = " << nexp+1 << "/" << fNumMockDataExpts << " *-*" << "\033[0m" << endl;
00212     MSG("GhostFramework",Msg::kInfo) << " " << endl;
00213 
00214 
00215     //Make FakeData histograms for this experiment (if using real data, is rebinned when fit is initialised)
00216     if(gInputs->fMCData){
00217 
00218       //gInputs->FitPureDecoherence();      
00219       gFakeData->GenerateSpectra((gInputs->vFinalData3D));
00220       //gInputs->FitPureDecoherence(fFitDecoherence);      
00221       
00222       MSG("GhostFramework",Msg::kInfo) << " *-* GhostFramework::RunFit() Generated Fake Data Spectra for all samples *-* " << endl;
00223     }
00224     
00225     //Setup likelihood surfaces
00226     TString name;
00227   
00228     //Set bin width for 2D "overall" contour plots
00229     fDm2BinWidth = (fDm2MaxOverall - fDm2MinOverall)/fNumDm2BinsOverall;
00230     fSin2BinWidth = (fSin2MaxOverall - fSin2MinOverall)/fNumSin2BinsOverall;
00231   
00232     name.Resize(0); name.Append("hLogLikelihoodSurface_Exp_"); name += nexp;
00233     hOverallSurface = new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall);
00234     if(!fRunParallel){
00235       name.Resize(0); name.Append("hDeltaLogLikelihoodSurface_Exp_"); name += nexp;
00236       hDeltaOverallSurface = new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall);
00237     }
00238     
00239     //NOTE: ALL SYST SURFACES FILLED IN UNITS OF SIGMA
00240     //Scaling Surfaces
00241     if(fFitCCDet){
00242       name.Resize(0); name.Append("hCCDetScalingSurface_Exp_"); name += nexp;
00243       hCCDetNormalisation = new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall);
00244       name.Resize(0); name.Append("hCCDetScalingErrorSurface__Exp_"); name += nexp;
00245       hCCDetNormalisationError = new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall);
00246     }
00247     if(fFitCCRock){
00248       name.Resize(0); name.Append("hCCRockScalingSurface_Exp_"); name += nexp;
00249       hCCRockNormalisation = new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall);
00250       name.Resize(0); name.Append("hCCRockScalingErrorSurface__Exp_"); name += nexp;
00251       hCCRockNormalisationError = new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall);
00252     }
00253     if(fFitNC){
00254       name.Resize(0); name.Append("hNCScalingSurface_Exp_"); name += nexp;
00255       hNCNormalisation = new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall);
00256       name.Resize(0); name.Append("hNCScalingErrorSurface__Exp_"); name += nexp;
00257       hNCNormalisationError = new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall);
00258     }
00259     if(fFitNuMuBar){
00260       name.Resize(0); name.Append("hNuMuBarScalingSurface_Exp_"); name += nexp;
00261       hNuMuBarNormalisation = new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall);
00262       name.Resize(0); name.Append("hNuMuBarScalingErrorSurface__Exp_"); name += nexp;
00263       hNuMuBarNormalisationError = new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall);
00264     }
00265     
00266     //2D Systematics Surfaces
00267     for(Int_t i=0;i<fNum2DSysts;i++){
00268       name.Resize(0); name.Append(gInputs->f2DSystematicName.at(i)); name.Append("Surface_Exp_"); name += nexp;
00269       v2DSystematicSurfaces.push_back(new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall));
00270       name.Resize(0); name.Append(gInputs->f2DSystematicName.at(i)); name.Append("ErrorSurface_Exp_"); name += nexp;
00271       v2DSystematicErrorSurfaces.push_back(new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall));
00272     }
00273     //1D Systematics Surfaces
00274     for(Int_t i=0;i<fNum1DSysts;i++){
00275       name.Resize(0); name.Append(gInputs->f1DSystematicName.at(i)); name.Append("Surface_Exp_"); name += nexp;
00276       v1DSystematicSurfaces.push_back(new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall));
00277       name.Resize(0); name.Append(gInputs->f1DSystematicName.at(i)); name.Append("ErrorSurface_Exp_"); name += nexp;
00278       v1DSystematicErrorSurfaces.push_back(new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall));
00279     }    
00280     if(fFitDecoherence){
00281       name.Resize(0); name.Append("hDecoherenceSurface_Exp_"); name += nexp;
00282       hDecoherenceSurface = new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall);
00283       name.Resize(0); name.Append("hDecoherenceErrorSurface__Exp_"); name += nexp;
00284       hDecoherenceSurfaceError = new TH2D(name.Data(),name.Data(),fNumSin2BinsOverall,fSin2MinOverall,fSin2MaxOverall,fNumDm2BinsOverall,fDm2MinOverall,fDm2MaxOverall);
00285     }
00286 
00287     MSG("GhostFramework",Msg::kDebug) << " *-* GhostFramework::RunFit() Prepared Likelihood Surfaces Completed *-*" << endl;
00288     
00289     this->OscFitGridSearch(); 
00290       
00291     MSG("OscFitFramework",Msg::kInfo) << " *-* OscFitFramework::RunFit() GridSearch Completed *-*" << endl;
00292     
00293   
00294     //IF NOT PARALLEL, FINALISE SURFACES AND BEST FIT VALUES TO WRITE OUT
00295     if(!fRunParallel){
00296       
00297       //Calculate best fit point and delta loglikelihood surface
00298       for(Int_t i=0;i<fNumSin2BinsOverall;i++){
00299         for(Int_t j=0;j<fNumDm2BinsOverall;j++){
00300           
00301           Double_t tmp_lik = hOverallSurface->GetBinContent(i+1,j+1);
00302           if((i==0 && j==0) || fMinLik > tmp_lik){
00303             fMinLik = tmp_lik; fMinIndexSin2 = i; fMinIndexDm2 = j;
00304           }
00305         }
00306       }
00307       MSG("GhostFramework",Msg::kInfo) << " *-* GhostFramework::RunFit() Found Best Fit Point *-*" << endl;
00308       
00309       //Fill delta loglikelihood surface and get best fit values
00310       for(Int_t i=0;i<fNumSin2Bins;i++){
00311         for(Int_t j=0;j<fNumDm2Bins;j++){
00312           Double_t temp_delta = 0.0;
00313           temp_delta = (hOverallSurface->GetBinContent(i+1,j+1) - fMinLik);
00314           hDeltaOverallSurface->SetBinContent(i+1,j+1,temp_delta);
00315         }
00316       }
00317       
00318       MSG("GhostFramework",Msg::kInfo) << " *-* GhostFramework::RunFit() Made DeltaLogLikelihood Surface *-*" << endl;
00319       
00320       //Get Variables for text file
00321       fMinLikSin2 = hOverallSurface->GetXaxis()->GetBinCenter(fMinIndexSin2+1); fMinLikSin2Error = fSin2BinWidth/2.;
00322       fMinLikDm2 = hOverallSurface->GetYaxis()->GetBinCenter(fMinIndexDm2+1); fMinLikDm2Error = fDm2BinWidth/2.;
00323       
00324       if(nexp==0){
00325         ofstream ofile(gInputs->fBestFitTextFilename);
00326         if(!ofile.fail()){
00327           ofile << "Number of Experiments = " << fNumMockDataExpts << endl;
00328           ofile << endl;
00329         }
00330         ofile.close();
00331       }
00332       
00333       ofstream ofile(gInputs->fBestFitTextFilename,ios::app);
00334       MSG("GhostFramework",Msg::kDebug) << " *-* GhostFramework::RunFit() Opened text file *-* " << endl;
00335       
00336       if(!ofile.fail()){
00337         ofile << "EXPERIMENT = " << nexp << endl;
00338         ofile << "Best Fit Point (sin2,dm2) = (" << fMinLikSin2 << "," << fMinLikDm2 << ")" << endl;
00339         ofile << "Error on sin2 best fit point: " << fMinLikSin2Error << endl;
00340         ofile << "Error on dm2 best fit point: " << fMinLikDm2Error << endl;
00341         ofile << "Likelihood value at the best fit point = " << fMinLik << endl;
00342         ofile << endl;
00343                  
00344         if(fFitCCDet){ofile << "Best Fit CC Detector Scale (/sigma) = " << hCCDetNormalisation->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1) << endl;
00345         ofile << "Error on Best Fit CC Detector Scale (/sigma) = " << hCCDetNormalisationError->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1) << endl;}
00346         if(fFitCCRock){ofile << "Best Fit CC Rock Scale (/sigma) = " << hCCRockNormalisation->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1) << endl;
00347         ofile << "Error on Best Fit CC Rock Scale (/sigma) = " << hCCRockNormalisationError->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1) << endl;}
00348         if(fFitNC){ofile << "Best Fit NC Scale (/sigma) = " << hNCNormalisation->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1) << endl;
00349         ofile << "Error on Best Fit NC Scale (/sigma) = " << hNCNormalisationError->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1) << endl;}
00350         if(fFitNuMuBar){ofile << "Best Fit NuMuBar Scale (/sigma) = " << hNuMuBarNormalisation->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1) << endl;
00351         ofile << "Error on Best Fit NuMuBar Scale (/sigma) = " << hNuMuBarNormalisationError->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1) << endl;}
00352         
00353         for(Int_t i=0;i<fNum2DSysts;i++){
00354           ofile << "Best Fit " << (gInputs->f2DSystematicName).at(i) <<  " (/sigma)  = " << v2DSystematicSurfaces.at(i)->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1) << endl;
00355           ofile << "Error on " << (gInputs->f2DSystematicName).at(i) <<  " (/sigma)  = " << v2DSystematicErrorSurfaces.at(i)->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1) << endl;
00356         }
00357        
00358         for(Int_t i=0;i<fNum1DSysts;i++){
00359           ofile << "Best Fit " << (gInputs->f1DSystematicName).at(i) <<  " (/sigma)  = " << v1DSystematicSurfaces.at(i)->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1) << endl;
00360           ofile << "Error on " << (gInputs->f1DSystematicName).at(i) <<  " (/sigma)  = " << v1DSystematicErrorSurfaces.at(i)->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1) << endl;
00361         }
00362 
00363         if(fFitDecoherence){
00364           ofile << "Best Fit mu2 = " << hDecoherenceSurface->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1) << endl;
00365           ofile << "Error on Best Fit mu2 = " << hDecoherenceSurfaceError->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1) << endl;
00366         }
00367       }
00368       ofile.close();
00369       
00370       MSG("GhostFramework",Msg::kInfo) << " *-* GhostFramework::RunFit() Written Text File *-* " << endl; 
00371       
00372       //When running MDC
00373       ofstream bffile("BestFitPoints.txt",ios::app);
00374       if(!bffile.fail()){
00375         bffile << fMinLikSin2 << "," << fMinLikDm2 << ";" << endl;
00376       }
00377       bffile.close();
00378 
00379       ofstream dmfile("BestFitDm2.txt",ios::app);
00380       if(!dmfile.fail()){
00381         dmfile << fMinLikDm2 << endl;
00382       }
00383       dmfile.close();
00384 
00385       ofstream sinfile("BestFitSin2.txt",ios::app);
00386       if(!sinfile.fail()){
00387         sinfile << fMinLikSin2 << endl;
00388       }
00389       sinfile.close();
00390 
00391       //Make best fit graph
00392       Double_t sin2_bestfit_array[1] = {fMinLikSin2};
00393       Double_t dm2_bestfit_array[1] = {fMinLikDm2};
00394       
00395       gBestFit = new TGraph(1,sin2_bestfit_array,dm2_bestfit_array);
00396       gBestFit->SetMarkerColor(2);
00397       gBestFit->SetMarkerStyle(5);
00398       gBestFit->SetMarkerSize(5);
00399       TString tempstring("gBestFit_Exp");
00400       tempstring += nexp;                         
00401       gBestFit->SetName(tempstring.Data());
00402 
00403     }//EO non-parallel run
00404     
00405     //WRITE OUT ROOT FILE-----------------------------------------------------------------------------------------------------------
00406     if(nexp==0 && gInputs->fRealData) fOutputRootFile = new TFile(gInputs->fBestFitRootFilename,"recreate");
00407     else fOutputRootFile = new TFile(gInputs->fBestFitRootFilename,"update");
00408     fOutputRootFile->cd();
00409     hOverallSurface->Write();
00410     if(fFitCCDet){hCCDetNormalisation->Write();
00411                   hCCDetNormalisationError->Write();}
00412     if(fFitCCRock){hCCRockNormalisation->Write();
00413                   hCCRockNormalisationError->Write();}
00414     if(fFitNC){hNCNormalisation->Write();
00415                hNCNormalisationError->Write();}
00416     if(fFitNuMuBar){hNuMuBarNormalisation->Write();
00417                hNuMuBarNormalisationError->Write();}
00418 
00419     for(Int_t i=0;i<fNum2DSysts;i++){
00420       v2DSystematicSurfaces.at(i)->Write();
00421       v2DSystematicErrorSurfaces.at(i)->Write();
00422     }
00423     for(Int_t i=0;i<fNum1DSysts;i++){
00424       v1DSystematicSurfaces.at(i)->Write();
00425       v1DSystematicErrorSurfaces.at(i)->Write();
00426     }
00427     if(fFitDecoherence){
00428       hDecoherenceSurface->Write();
00429       hDecoherenceSurfaceError->Write();
00430     }
00431     if(!fRunParallel){
00432       gBestFit->Write();
00433       hDeltaOverallSurface->Write();
00434     }
00435     //Write out Fake/Real data spectra for comparison purposes
00436     if(!fRunParallel){
00437       for(UInt_t i=0;i<(gInputs->vFinalData3D).size();i++){
00438         for(UInt_t j=0;j<((gInputs->vFinalData3D).at(i)).size();j++){
00439           for(UInt_t k=0;k<(((gInputs->vFinalData3D).at(i)).at(j)).size();k++){
00440           
00441             TString name("hDataSpectrum_"); name += i; name.Append("_"); name += j; name.Append("_"); name += k;
00442             TH1D* htemp = (((gInputs->vFinalData3D).at(i)).at(j)).at(k)->ConvertTH1D(name.Data());
00443             htemp->Write();
00444             delete htemp;
00445           }
00446         }
00447       }
00448     }
00449     //Write out spectrum at best fit point for each sample
00450     for(Int_t i=0;i<fNumSamples;i++){
00451 
00452       TString label("BESTFIT");
00453       TString label_unosc("UNOSCILLATED");
00454 
00455       //Dummy 2D vector container
00456       vector<vector<ArrayTH1D*> > vDummy2D;
00457       vector<ArrayTH1D*> vDummy1D;
00458       for(Int_t j=0;j<vSamples.at(i)->fNDis;j++) vDummy2D.push_back(vDummy1D);
00459 
00460       //Get scales
00461       Double_t cc_detbestfitscale = 0.0;
00462       Double_t cc_rockbestfitscale = 0.0;
00463       Double_t nc_bestfitscale = 0.0;
00464       Double_t ws_bestfitscale = 0.0;
00465 
00466       if(fFitCCDet) cc_detbestfitscale = hCCDetNormalisation->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1)*fCCDetSigma;
00467       if(fFitCCRock) cc_rockbestfitscale = hCCRockNormalisation->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1)*fCCRockSigma;
00468       if(fFitNC) nc_bestfitscale = hNCNormalisation->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1)*fNCSigma;
00469       if(fFitNuMuBar) ws_bestfitscale = hNuMuBarNormalisation->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1)*fNuMuBarSigma;
00470       
00471       //Get energy systs
00472       Double_t* EnergySysts2D = new Double_t[fNum2DSysts];
00473       Double_t* EnergySysts1D = new Double_t[fNum1DSysts];
00474       for(Int_t j=0;j<fNum2DSysts;j++) EnergySysts2D[j] = v2DSystematicSurfaces.at(j)->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1);
00475       for(Int_t j=0;j<fNum1DSysts;j++) EnergySysts1D[j] = v1DSystematicSurfaces.at(j)->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1);
00476 
00477       Double_t mu2_bestfit = 0.0;
00478       if(fFitDecoherence) mu2_bestfit = hDecoherenceSurface->GetBinContent(fMinIndexSin2+1,fMinIndexDm2+1);
00479 
00480       MSG("GhostFramework",Msg::kInfo) << " *-* GhostFramework::RunFit() Writing out best fit spectra *-*" << endl;
00481       vSamples.at(i)->GetSpectrum(fMinLikSin2,fMinLikDm2,cc_detbestfitscale,cc_rockbestfitscale,nc_bestfitscale,ws_bestfitscale,EnergySysts2D,EnergySysts1D,vDummy2D,fOutputRootFile,label,mu2_bestfit);
00482 
00483       //Write out overall best fit spectrum for sample,dis,res bin
00484       for(Int_t j=0;j<vSamples.at(i)->fNDis;j++){
00485         for(Int_t k=0;k<vSamples.at(i)->fNResBins;k++){
00486           
00487           TString name; name.Append("hBestFitSpectrum_"); name += i; name.Append("_"); name += j; name.Append("_"); name += k;
00488           TH1D* hBestFit = (TH1D*)(vDummy2D.at(j)).at(k)->ConvertTH1D(name.Data());
00489           hBestFit->Write();
00490           delete hBestFit;
00491         }
00492       }
00493 
00494       //Get Unoscillated prediction for all contributions
00495       vector<vector<ArrayTH1D*> > vDummyUnosc2D;
00496       vector<ArrayTH1D*> vDummyUnosc1D;
00497       for(Int_t j=0;j<vSamples.at(i)->fNDis;j++) vDummyUnosc2D.push_back(vDummyUnosc1D);
00498 
00499       for(Int_t j=0;j<fNum2DSysts;j++) EnergySysts2D[j] = 0.;
00500       for(Int_t j=0;j<fNum1DSysts;j++) EnergySysts1D[j] = 0.;
00501       vSamples.at(i)->GetSpectrum(0.0,0.0,0.0,0.0,0.0,0.0,EnergySysts2D,EnergySysts1D,vDummyUnosc2D,fOutputRootFile,label_unosc);
00502     
00503       //Write out unoscillated spectrum for sample,dis,res bin
00504       for(Int_t j=0;j<vSamples.at(i)->fNDis;j++){
00505         for(Int_t k=0;k<vSamples.at(i)->fNResBins;k++){
00506           
00507           TString name; name.Append("hUnoscillatedSpectrum_"); name += i; name.Append("_"); name += j; name.Append("_"); name += k;
00508           TH1D* hUnoscillated = (TH1D*)(vDummyUnosc2D.at(j)).at(k)->ConvertTH1D(name.Data());
00509           hUnoscillated->Write();
00510           delete hUnoscillated;
00511         }
00512       }
00513       delete[] EnergySysts2D ;
00514       delete[] EnergySysts1D;
00515     }
00516     
00517     fOutputRootFile->Close();
00518      
00519     MSG("GhostFramework",Msg::kInfo) << " *-* GhostFramework::RunFit() Written Root File *-* " << endl; 
00520 
00521     //Remade for each new experiment
00522     delete gBestFit;
00523     delete hOverallSurface;
00524     if(!fRunParallel) delete hDeltaOverallSurface;
00525     if(fFitCCDet){
00526       delete hCCDetNormalisation;
00527       delete hCCDetNormalisationError;
00528     }
00529     if(fFitCCRock){
00530       delete hCCRockNormalisation;
00531       delete hCCRockNormalisationError;
00532     }
00533     if(fFitNC){
00534       delete hNCNormalisation;
00535       delete hNCNormalisationError;
00536     }
00537     if(fFitNuMuBar){
00538       delete hNuMuBarNormalisation;
00539       delete hNuMuBarNormalisationError;
00540     }
00541     for(Int_t i=0;i<fNum2DSysts;i++){
00542       delete v2DSystematicSurfaces.at(i);
00543       delete v2DSystematicErrorSurfaces.at(i);
00544     }
00545     for(Int_t i=0;i<fNum1DSysts;i++){
00546       delete v1DSystematicSurfaces.at(i);
00547       delete v1DSystematicErrorSurfaces.at(i);
00548     }
00549     if(fFitDecoherence){
00550       delete hDecoherenceSurface;
00551       delete hDecoherenceSurfaceError;
00552     }
00553     delete fOutputRootFile;
00554   }//EO loop over experiments
00555 }

void GhostFramework::SetChi2Plot ( Bool_t  plot  )  [inline]

Definition at line 43 of file GhostFramework.h.

References fPlot.

00043 { fPlot = plot; }

void GhostFramework::UseResolutionLikelihood ( int  reslik = 1  )  [inline]

Definition at line 34 of file GhostFramework.h.

References fResolutionLikelihood, and fStandardLikelihood.

void GhostFramework::UseStandardLikelihood ( int  stdlik = 1  )  [inline]

Definition at line 35 of file GhostFramework.h.

References fResolutionLikelihood, and fStandardLikelihood.


Member Data Documentation

Definition at line 113 of file GhostFramework.h.

Referenced by GhostFramework(), LogLikelihood(), and ~GhostFramework().

Definition at line 112 of file GhostFramework.h.

Referenced by GhostFramework(), LogLikelihood(), and ~GhostFramework().

Definition at line 110 of file GhostFramework.h.

Referenced by GhostFramework(), LogLikelihood(), and ~GhostFramework().

Definition at line 109 of file GhostFramework.h.

Referenced by GhostFramework(), LogLikelihood(), and ~GhostFramework().

Double_t GhostFramework::fCCDetScale [private]

Definition at line 92 of file GhostFramework.h.

Referenced by GhostFramework().

Double_t GhostFramework::fCCDetSigma [private]
Double_t GhostFramework::fCCRockScale [private]

Definition at line 93 of file GhostFramework.h.

Referenced by GhostFramework().

Double_t GhostFramework::fCCRockSigma [private]
Int_t GhostFramework::fDecoPower [private]

Definition at line 119 of file GhostFramework.h.

Referenced by GhostFramework().

Double_t GhostFramework::fDm2_par [private]

Definition at line 104 of file GhostFramework.h.

Referenced by LogLikelihood(), and LogLikelihoodNuisance().

Double_t GhostFramework::fDm2BinWidth [private]

Definition at line 69 of file GhostFramework.h.

Referenced by OscFitGridSearch(), and RunFit().

Double_t GhostFramework::fDm2Max [private]

Definition at line 84 of file GhostFramework.h.

Referenced by GhostFramework(), and OscFitGridSearch().

Double_t GhostFramework::fDm2MaxOverall [private]

Definition at line 74 of file GhostFramework.h.

Referenced by GhostFramework(), and RunFit().

Double_t GhostFramework::fDm2Min [private]

Definition at line 83 of file GhostFramework.h.

Referenced by GhostFramework(), and OscFitGridSearch().

Double_t GhostFramework::fDm2MinOverall [private]

Definition at line 73 of file GhostFramework.h.

Referenced by GhostFramework(), and RunFit().

Bool_t GhostFramework::fFitCCDet [private]
Bool_t GhostFramework::fFitCCRock [private]
Bool_t GhostFramework::fFitDecay [private]

Definition at line 61 of file GhostFramework.h.

Referenced by GhostFramework().

Bool_t GhostFramework::fFitNC [private]
Bool_t GhostFramework::fFitNuMuBar [private]
Bool_t GhostFramework::fFitScales [private]

Definition at line 65 of file GhostFramework.h.

Referenced by GhostFramework().

TMinuit* GhostFramework::fFreeMinuit [private]

Definition at line 121 of file GhostFramework.h.

Referenced by LogLikelihood().

TMinuit* GhostFramework::fInitialNMinuit [private]

Definition at line 120 of file GhostFramework.h.

Referenced by FindInitialNuisanceValues().

Definition at line 128 of file GhostFramework.h.

Referenced by RunFit().

Definition at line 129 of file GhostFramework.h.

Referenced by RunFit().

Double_t GhostFramework::fMinLik [private]

Definition at line 127 of file GhostFramework.h.

Referenced by RunFit().

Double_t GhostFramework::fMinLikDm2 [private]

Definition at line 132 of file GhostFramework.h.

Referenced by RunFit().

Double_t GhostFramework::fMinLikDm2Error [private]

Definition at line 133 of file GhostFramework.h.

Referenced by RunFit().

Double_t GhostFramework::fMinLikSin2 [private]

Definition at line 130 of file GhostFramework.h.

Referenced by RunFit().

Definition at line 131 of file GhostFramework.h.

Referenced by RunFit().

Double_t* GhostFramework::fMinScale [private]

Definition at line 106 of file GhostFramework.h.

Referenced by GhostFramework(), LogLikelihood(), and ~GhostFramework().

Double_t* GhostFramework::fMinScaleError [private]

Definition at line 107 of file GhostFramework.h.

Referenced by GhostFramework(), LogLikelihood(), and ~GhostFramework().

Definition at line 66 of file GhostFramework.h.

Referenced by GhostFramework(), and LogLikelihood().

Double_t GhostFramework::fMu2_Error [private]

Definition at line 117 of file GhostFramework.h.

Referenced by GhostFramework(), and LogLikelihood().

Double_t GhostFramework::fMu2_High [private]

Definition at line 115 of file GhostFramework.h.

Referenced by FindInitialNuisanceValues(), GhostFramework(), and LogLikelihood().

Double_t GhostFramework::fMu2_Initial [private]

Definition at line 114 of file GhostFramework.h.

Referenced by FindInitialNuisanceValues(), GhostFramework(), and LogLikelihood().

Double_t GhostFramework::fMu2_par [private]

Definition at line 116 of file GhostFramework.h.

Referenced by GhostFramework(), and LogLikelihood().

Double_t GhostFramework::fNCScale [private]

Definition at line 94 of file GhostFramework.h.

Referenced by GhostFramework().

Double_t GhostFramework::fNCSigma [private]

Definition at line 86 of file GhostFramework.h.

Referenced by GhostFramework(), and LogLikelihoodNuisance().

Int_t GhostFramework::fNum1DSysts [private]
Int_t GhostFramework::fNum2DSysts [private]
Int_t GhostFramework::fNumDm2Bins [private]

Definition at line 80 of file GhostFramework.h.

Referenced by GhostFramework(), OscFitGridSearch(), and RunFit().

Definition at line 75 of file GhostFramework.h.

Referenced by GhostFramework(), and RunFit().

Definition at line 67 of file GhostFramework.h.

Referenced by GhostFramework(), and RunFit().

Int_t GhostFramework::fNumSamples [private]

Definition at line 79 of file GhostFramework.h.

Referenced by GhostFramework(), OscFitGridSearch(), and RunFit().

Definition at line 76 of file GhostFramework.h.

Referenced by GhostFramework(), and RunFit().

Double_t GhostFramework::fNuMuBarScale [private]

Definition at line 95 of file GhostFramework.h.

Referenced by GhostFramework().

Double_t GhostFramework::fNuMuBarSigma [private]

Definition at line 136 of file GhostFramework.h.

Referenced by RunFit().

Bool_t GhostFramework::fPlot [private]

Definition at line 163 of file GhostFramework.h.

Referenced by LogLikelihood(), OscFitGridSearch(), and SetChi2Plot().

Bool_t GhostFramework::fRunParallel [private]

Definition at line 64 of file GhostFramework.h.

Referenced by GhostFramework(), and RunFit().

Double_t GhostFramework::fSin2_par [private]

Definition at line 103 of file GhostFramework.h.

Referenced by LogLikelihood(), and LogLikelihoodNuisance().

Double_t GhostFramework::fSin2BinWidth [private]

Definition at line 70 of file GhostFramework.h.

Referenced by OscFitGridSearch(), and RunFit().

Double_t GhostFramework::fSin2Max [private]

Definition at line 82 of file GhostFramework.h.

Referenced by GhostFramework(), and OscFitGridSearch().

Double_t GhostFramework::fSin2MaxOverall [private]

Definition at line 72 of file GhostFramework.h.

Referenced by GhostFramework(), and RunFit().

Double_t GhostFramework::fSin2Min [private]

Definition at line 81 of file GhostFramework.h.

Referenced by GhostFramework(), and OscFitGridSearch().

Double_t GhostFramework::fSin2MinOverall [private]

Definition at line 71 of file GhostFramework.h.

Referenced by GhostFramework(), and RunFit().

Definition at line 63 of file GhostFramework.h.

Referenced by UseResolutionLikelihood(), and UseStandardLikelihood().

Definition at line 118 of file GhostFramework.h.

Referenced by GhostFramework(), LogLikelihood(), and LogLikelihoodNuisance().

Double_t GhostFramework::fTrackGlobal [private]

Definition at line 164 of file GhostFramework.h.

Referenced by LogLikelihood(), and OscFitGridSearch().

TGraph* GhostFramework::gBestFit [private]

Definition at line 134 of file GhostFramework.h.

Referenced by RunFit().

Definition at line 57 of file GhostFramework.h.

Referenced by GhostFramework(), and RunFit().

Definition at line 145 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

Definition at line 146 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

Definition at line 147 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

Definition at line 148 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

Definition at line 153 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

Definition at line 154 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

Definition at line 143 of file GhostFramework.h.

Referenced by RunFit().

Double_t* GhostFramework::highphysbound [private]

Definition at line 126 of file GhostFramework.h.

Referenced by FindInitialNuisanceValues(), and LogLikelihood().

Definition at line 151 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

Definition at line 152 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

Definition at line 149 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

Definition at line 150 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

Definition at line 142 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

TH1D* GhostFramework::hShwChi2 [private]

Definition at line 165 of file GhostFramework.h.

TH2D* GhostFramework::hTrkChi2 [private]

Definition at line 166 of file GhostFramework.h.

Referenced by OscFitGridSearch().

Double_t* GhostFramework::lowphysbound [private]

Definition at line 125 of file GhostFramework.h.

Referenced by FindInitialNuisanceValues(), and LogLikelihood().

string* GhostFramework::parName [private]

Definition at line 124 of file GhostFramework.h.

Referenced by FindInitialNuisanceValues(), and LogLikelihood().

Double_t* GhostFramework::step [private]

Definition at line 123 of file GhostFramework.h.

Referenced by FindInitialNuisanceValues(), and LogLikelihood().

Definition at line 159 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

Definition at line 158 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

Definition at line 157 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

Definition at line 156 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

Double_t* GhostFramework::vstart [private]

Definition at line 122 of file GhostFramework.h.

Referenced by FindInitialNuisanceValues(), and LogLikelihood().


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

Generated on 14 Dec 2017 for loon by  doxygen 1.6.1