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, GhostInputs::fCCDetScale, fCCDetScale, GhostInputs::fCCDetSigma, fCCDetSigma, GhostInputs::fCCRockScale, fCCRockScale, GhostInputs::fCCRockSigma, fCCRockSigma, GhostInputs::fDecoPower, fDecoPower, GhostInputs::fDm2Max, fDm2Max, GhostInputs::fDm2MaxOverall, fDm2MaxOverall, GhostInputs::fDm2Min, fDm2Min, GhostInputs::fDm2MinOverall, fDm2MinOverall, GhostInputs::fFitCCDet, fFitCCDet, GhostInputs::fFitCCRock, fFitCCRock, GhostInputs::fFitDecay, fFitDecay, GhostInputs::fFitDecoherence, fFitDecoherence, GhostInputs::fFitNC, fFitNC, GhostInputs::fFitNuMuBar, fFitNuMuBar, GhostInputs::fFitScales, fFitScales, fInitialScalingFactors, fMinScale, fMinScaleError, GhostInputs::fMinuitInitialValues, fMinuitInitialValues, fMu2_Error, fMu2_High, fMu2_Initial, fMu2_par, GhostInputs::fNCScale, fNCScale, GhostInputs::fNCSigma, fNCSigma, GhostInputs::fNuisanceParameters, fNuisanceParameters, GhostInputs::fNum1DSysts, fNum1DSysts, GhostInputs::fNum2DSysts, fNum2DSysts, GhostInputs::fNumDm2Bins, fNumDm2Bins, GhostInputs::fNumDm2BinsOverall, fNumDm2BinsOverall, GhostInputs::fNumMockDataExpts, fNumMockDataExpts, GhostInputs::fNumSamples, fNumSamples, GhostInputs::fNumScalingFactors, fNumScalingFactors, GhostInputs::fNumSin2Bins, fNumSin2Bins, GhostInputs::fNumSin2BinsOverall, fNumSin2BinsOverall, GhostInputs::fNuMuBarScale, fNuMuBarScale, fNuMuBarSigma, GhostInputs::fRunParallel, fRunParallel, GhostInputs::fSin2Max, fSin2Max, GhostInputs::fSin2MaxOverall, fSin2MaxOverall, GhostInputs::fSin2Min, fSin2Min, GhostInputs::fSin2MinOverall, 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, BeamData::error(), 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, parName, step, and vstart.

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, 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, parName, 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, 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, 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(), 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

Double_t* GhostFramework::f1DInitialSystematicShift [private]

Definition at line 111 of file GhostFramework.h.

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

Double_t* GhostFramework::f1DMinErrorSystematicShift [private]

Definition at line 113 of file GhostFramework.h.

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

Double_t* GhostFramework::f1DMinSystematicShift [private]

Definition at line 112 of file GhostFramework.h.

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

Double_t* GhostFramework::f2DInitialSystematicShift [private]

Definition at line 108 of file GhostFramework.h.

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

Double_t* GhostFramework::f2DMinErrorSystematicShift [private]

Definition at line 110 of file GhostFramework.h.

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

Double_t* GhostFramework::f2DMinSystematicShift [private]

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]

Definition at line 88 of file GhostFramework.h.

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

Double_t GhostFramework::fCCRockScale [private]

Definition at line 93 of file GhostFramework.h.

Referenced by GhostFramework().

Double_t GhostFramework::fCCRockSigma [private]

Definition at line 89 of file GhostFramework.h.

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

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]

Definition at line 96 of file GhostFramework.h.

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

Bool_t GhostFramework::fFitCCRock [private]

Definition at line 97 of file GhostFramework.h.

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

Bool_t GhostFramework::fFitDecay [private]

Definition at line 61 of file GhostFramework.h.

Referenced by GhostFramework().

Bool_t GhostFramework::fFitDecoherence [private]

Definition at line 60 of file GhostFramework.h.

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

Bool_t GhostFramework::fFitNC [private]

Definition at line 98 of file GhostFramework.h.

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

Bool_t GhostFramework::fFitNuMuBar [private]

Definition at line 99 of file GhostFramework.h.

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

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().

Double_t* GhostFramework::fInitialScalingFactors [private]

Definition at line 105 of file GhostFramework.h.

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

Int_t GhostFramework::fMinIndexDm2 [private]

Definition at line 128 of file GhostFramework.h.

Referenced by RunFit().

Int_t GhostFramework::fMinIndexSin2 [private]

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().

Double_t GhostFramework::fMinLikSin2Error [private]

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().

Bool_t GhostFramework::fMinuitInitialValues [private]

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 90 of file GhostFramework.h.

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

Bool_t GhostFramework::fNuisanceParameters [private]

Definition at line 86 of file GhostFramework.h.

Referenced by GhostFramework(), and LogLikelihoodNuisance().

Int_t GhostFramework::fNum1DSysts [private]

Definition at line 78 of file GhostFramework.h.

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

Int_t GhostFramework::fNum2DSysts [private]

Definition at line 77 of file GhostFramework.h.

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

Int_t GhostFramework::fNumDm2Bins [private]

Definition at line 80 of file GhostFramework.h.

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

Int_t GhostFramework::fNumDm2BinsOverall [private]

Definition at line 75 of file GhostFramework.h.

Referenced by GhostFramework(), and RunFit().

Int_t GhostFramework::fNumMockDataExpts [private]

Definition at line 67 of file GhostFramework.h.

Referenced by GhostFramework(), and RunFit().

Int_t GhostFramework::fNumSamples [private]

Definition at line 68 of file GhostFramework.h.

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

Int_t GhostFramework::fNumScalingFactors [private]

Definition at line 85 of file GhostFramework.h.

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

Int_t GhostFramework::fNumSin2Bins [private]

Definition at line 79 of file GhostFramework.h.

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

Int_t GhostFramework::fNumSin2BinsOverall [private]

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 91 of file GhostFramework.h.

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

TFile* GhostFramework::fOutputRootFile [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::fResolutionLikelihood [private]

Definition at line 62 of file GhostFramework.h.

Referenced by LogLikelihood(), LogLikelihoodNuisance(), UseResolutionLikelihood(), and UseStandardLikelihood().

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().

Bool_t GhostFramework::fStandardLikelihood [private]

Definition at line 63 of file GhostFramework.h.

Referenced by UseResolutionLikelihood(), and UseStandardLikelihood().

Int_t GhostFramework::fSurfacePair [private]

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().

GhostFakeData* GhostFramework::gFakeData [private]

Definition at line 57 of file GhostFramework.h.

Referenced by GhostFramework(), and RunFit().

GhostInputs* GhostFramework::gInputs [private]

Definition at line 56 of file GhostFramework.h.

Referenced by FindInitialNuisanceValues(), GhostFramework(), LogLikelihood(), LogLikelihoodNuisance(), RebinRealData(), and RunFit().

TH2D* GhostFramework::hCCDetNormalisation [private]

Definition at line 145 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

TH2D* GhostFramework::hCCDetNormalisationError [private]

Definition at line 146 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

TH2D* GhostFramework::hCCRockNormalisation [private]

Definition at line 147 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

TH2D* GhostFramework::hCCRockNormalisationError [private]

Definition at line 148 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

TH2D* GhostFramework::hDecoherenceSurface [private]

Definition at line 153 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

TH2D* GhostFramework::hDecoherenceSurfaceError [private]

Definition at line 154 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

TH2D* GhostFramework::hDeltaOverallSurface [private]

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().

TH2D* GhostFramework::hNCNormalisation [private]

Definition at line 151 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

TH2D* GhostFramework::hNCNormalisationError [private]

Definition at line 152 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

TH2D* GhostFramework::hNuMuBarNormalisation [private]

Definition at line 149 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

TH2D* GhostFramework::hNuMuBarNormalisationError [private]

Definition at line 150 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

TH2D* GhostFramework::hOverallSurface [private]

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().

vector<TH2D*> GhostFramework::v1DSystematicErrorSurfaces [private]

Definition at line 159 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

vector<TH2D*> GhostFramework::v1DSystematicSurfaces [private]

Definition at line 158 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

vector<TH2D*> GhostFramework::v2DSystematicErrorSurfaces [private]

Definition at line 157 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

vector<TH2D*> GhostFramework::v2DSystematicSurfaces [private]

Definition at line 156 of file GhostFramework.h.

Referenced by LogLikelihood(), and RunFit().

vector<GhostSample*> GhostFramework::vSamples [private]

Definition at line 139 of file GhostFramework.h.

Referenced by AddSample(), GhostFramework(), LogLikelihood(), LogLikelihoodNuisance(), RebinRealData(), 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 Mon Sep 1 00:51:38 2014 for loon by  doxygen 1.4.7