00001
00002
00003
00004
00005 #include "MakeSterilePred_Interpolated.h"
00006 #include "OscProb/NuOscProbCalc.h"
00007 #include "NtupleUtils/NuXMLConfig.h"
00008 #include "NtupleUtils/NuMMParameters.h"
00009
00010 #include <cassert>
00011 #include <iostream>
00012 #include <TFile.h>
00013 #include <TH1D.h>
00014 #include <TH2D.h>
00015 #include <TAxis.h>
00016 #include <TGraph.h>
00017
00018 using namespace std;
00019
00020 ClassImp(MakeSterilePred_Interpolated)
00021
00022
00023 MakeSterilePred_Interpolated::MakeSterilePred_Interpolated() : oscAvgOverBinNC(0), oscAvgOverBinNuMu(0),
00024 oscIterOverBinNC(0), oscIterOverBinNuMu(0),
00025 ND_MC(0), ND_Nue(0), ND_NuTau(0), ND_data(0),
00026 FD_MC(0), FD_Nue(0), FD_NuTau(0), FD_data(0),
00027 POT(0.0), got_histos(false), fOscBinCenter(false)
00028 {
00029 ND_baseline = 1.04*Munits::km;
00030 FD_baseline = 735.00*Munits::km;
00031
00032 fParMap[kShwScale] = NuMMParameters::kShwScale;
00033 fParMap[kTrkScale] = NuMMParameters::kTrkScale;
00034 fParMap[kNcBack] = NuMMParameters::kNcBack;
00035 fParMap[kAbsShwNC] = NuMMParameters::kAbsShwNC;
00036 fParMap[kRelShwNC] = NuMMParameters::kRelShwNC;
00037 fParMap[kCCBack] = NuMMParameters::kCCBack;
00038 fParMap[kFDClean] = NuMMParameters::kFDClean;
00039 fParMap[kNDClean] = NuMMParameters::kNDClean;
00040 fParMap[kFDCosmic] = NuMMParameters::kFDCosmic;
00041
00042 fFitSysts = false;
00043
00044 }
00045
00046
00047 void MakeSterilePred_Interpolated::cleanup()
00048 {
00049 if(oscAvgOverBinNC) { delete oscAvgOverBinNC; oscAvgOverBinNC = 0; }
00050 if(oscAvgOverBinNuMu) { delete oscAvgOverBinNuMu; oscAvgOverBinNuMu = 0; }
00051 if(oscIterOverBinNC) { delete oscIterOverBinNC; oscIterOverBinNC = 0; }
00052 if(oscIterOverBinNuMu) { delete oscIterOverBinNuMu; oscIterOverBinNuMu = 0; }
00053 }
00054
00055
00056 MakeSterilePred_Interpolated::~MakeSterilePred_Interpolated()
00057 {
00058 cleanup();
00059 }
00060
00061
00062 void MakeSterilePred_Interpolated::SetNDMC(TFile *MC_file, TFile *Nue_file, TFile *NuTau_file)
00063 {
00064 ND_MC = MC_file;
00065 ND_Nue = Nue_file;
00066 ND_NuTau = NuTau_file;
00067 }
00068
00069
00070 void MakeSterilePred_Interpolated::SetNDData(TFile *Data_file)
00071 {
00072 ND_data = Data_file;
00073 }
00074
00075
00076 void MakeSterilePred_Interpolated::SetNDBaseline(double baseline)
00077 {
00078 ND_baseline = baseline;
00079 }
00080
00081
00082 void MakeSterilePred_Interpolated::SetFDMC(TFile *MC_file, TFile *Nue_file, TFile *NuTau_file)
00083 {
00084 FD_MC = MC_file;
00085 FD_Nue = Nue_file;
00086 FD_NuTau = NuTau_file;
00087 }
00088
00089
00090 void MakeSterilePred_Interpolated::SetFDData(TFile *Data_file)
00091 {
00092 FD_data = Data_file;
00093 }
00094
00095
00096 void MakeSterilePred_Interpolated::SetFDBaseline(double baseline)
00097 {
00098 FD_baseline = baseline;
00099 }
00100
00101
00102 void MakeSterilePred_Interpolated::SetPOT(double POT_)
00103 {
00104 POT = POT_;
00105 }
00106
00107
00108 void MakeSterilePred_Interpolated::GetHistos()
00109 {
00110 if(got_histos) return;
00111
00112
00113 assert(FD_data);
00114 assert(ND_data);
00115 assert(ND_MC);
00116 assert(ND_Nue);
00117 assert(ND_NuTau);
00118 assert(FD_MC);
00119 assert(FD_Nue);
00120 assert(FD_NuTau);
00121
00122 FD_data->GetObject("hRecoEnergyNC", FD_dataNC); assert(FD_dataNC);
00123 FD_data->GetObject("hRecoEnergyCCAll", FD_dataCC); assert(FD_dataCC);
00124
00125 ND_data->GetObject("hRecoEnergyNC", ND_dataNC); assert(ND_dataNC);
00126 ND_data->GetObject("hRecoEnergyCCAll", ND_dataCC); assert(ND_dataCC);
00127
00128 ND_MC->GetObject("hRecoToTrueNCSelectedTrueNC", NDNC_TrueNC); assert(NDNC_TrueNC);
00129 ND_MC->GetObject("hRecoToTrueNCSelectedTrueCCNuMu", NDNC_NuMu); assert(NDNC_NuMu);
00130 ND_MC->GetObject("hRecoToTrueNCSelectedBeamNue", NDNC_BeamNue); assert(NDNC_BeamNue);
00131 ND_Nue->GetObject("hRecoToTrueNCSelectedAppearNue", NDNC_AppNue); assert(NDNC_AppNue);
00132 ND_NuTau->GetObject("hRecoToTrueNCSelectedAppearNuTau", NDNC_AppNuTau); assert(NDNC_AppNuTau);
00133
00134 ND_MC->GetObject("hRecoToTrueCCSelectedTrueNC", NDCC_TrueNC); assert(NDCC_TrueNC);
00135 ND_MC->GetObject("hRecoToTrueCCSelectedTrueCCNuMu", NDCC_NuMu); assert(NDCC_NuMu);
00136 ND_MC->GetObject("hRecoToTrueCCSelectedBeamNue", NDCC_BeamNue); assert(NDCC_BeamNue);
00137 ND_Nue->GetObject("hRecoToTrueCCSelectedAppearNue", NDCC_AppNue); assert(NDCC_AppNue);
00138 ND_NuTau->GetObject("hRecoToTrueCCSelectedAppearNuTau", NDCC_AppNuTau); assert(NDCC_AppNuTau);
00139
00140 FD_MC->GetObject("hRecoToTrueNCSelectedTrueNC", FDNC_TrueNC); assert(FDNC_TrueNC);
00141 FD_MC->GetObject("hRecoToTrueNCSelectedTrueCCNuMu", FDNC_NuMu); assert(FDNC_NuMu);
00142 FD_MC->GetObject("hRecoToTrueNCSelectedBeamNue", FDNC_BeamNue); assert(FDNC_BeamNue);
00143 FD_Nue->GetObject("hRecoToTrueNCSelectedAppearNue", FDNC_AppNue); assert(FDNC_AppNue);
00144 FD_NuTau->GetObject("hRecoToTrueNCSelectedAppearNuTau", FDNC_AppNuTau); assert(FDNC_AppNuTau);
00145
00146 FD_MC->GetObject("hRecoToTrueCCSelectedTrueNC", FDCC_TrueNC); assert(FDCC_TrueNC);
00147 FD_MC->GetObject("hRecoToTrueCCSelectedTrueCCNuMu", FDCC_NuMu); assert(FDCC_NuMu);
00148 FD_MC->GetObject("hRecoToTrueCCSelectedBeamNue", FDCC_BeamNue); assert(FDCC_BeamNue);
00149 FD_Nue->GetObject("hRecoToTrueCCSelectedAppearNue", FDCC_AppNue); assert(FDCC_AppNue);
00150 FD_NuTau->GetObject("hRecoToTrueCCSelectedAppearNuTau", FDCC_AppNuTau); assert(FDCC_AppNuTau);
00151
00152 ND_data->GetObject("hTotalPot", ND_data_POT); assert(ND_data_POT);
00153 ND_MC->GetObject("hTotalPot", ND_MC_POT); assert(ND_MC_POT);
00154 ND_Nue->GetObject("hTotalPot", ND_Nue_POT); assert(ND_Nue_POT);
00155 ND_NuTau->GetObject("hTotalPot", ND_NuTau_POT); assert(ND_NuTau_POT);
00156 FD_MC->GetObject("hTotalPot", FD_MC_POT); assert(FD_MC_POT);
00157 FD_Nue->GetObject("hTotalPot", FD_Nue_POT); assert(FD_Nue_POT);
00158 FD_NuTau->GetObject("hTotalPot", FD_NuTau_POT); assert(FD_NuTau_POT);
00159
00160
00161 double dataPOT = ND_data_POT->Integral();
00162 ND_dataNC->Scale(POT/dataPOT);
00163 ND_dataCC->Scale(POT/dataPOT);
00164
00165 double ndmcPOT = ND_MC_POT->Integral();
00166 NDNC_TrueNC->Scale(POT/ndmcPOT);
00167 NDNC_NuMu->Scale(POT/ndmcPOT);
00168 NDNC_BeamNue->Scale(POT/ndmcPOT);
00169 NDCC_TrueNC->Scale(POT/ndmcPOT);
00170 NDCC_NuMu->Scale(POT/ndmcPOT);
00171 NDCC_BeamNue->Scale(POT/ndmcPOT);
00172
00173 double ndnuePOT = ND_Nue_POT->Integral();
00174 NDNC_AppNue->Scale(POT/ndnuePOT);
00175 NDCC_AppNue->Scale(POT/ndnuePOT);
00176
00177 double ndnutauPOT = ND_NuTau_POT->Integral();
00178 NDNC_AppNuTau->Scale(POT/ndnutauPOT);
00179 NDCC_AppNuTau->Scale(POT/ndnutauPOT);
00180
00181 double fdmcPOT = FD_MC_POT->Integral();
00182 FDNC_TrueNC->Scale(POT/fdmcPOT);
00183 FDNC_NuMu->Scale(POT/fdmcPOT);
00184 FDNC_BeamNue->Scale(POT/fdmcPOT);
00185 FDCC_TrueNC->Scale(POT/fdmcPOT);
00186 FDCC_NuMu->Scale(POT/fdmcPOT);
00187 FDCC_BeamNue->Scale(POT/fdmcPOT);
00188
00189 double fdnuePOT = FD_Nue_POT->Integral();
00190 FDNC_AppNue->Scale(POT/fdnuePOT);
00191 FDCC_AppNue->Scale(POT/fdnuePOT);
00192
00193 double fdnutauPOT = FD_NuTau_POT->Integral();
00194 FDNC_AppNuTau->Scale(POT/fdnutauPOT);
00195 FDCC_AppNuTau->Scale(POT/fdnutauPOT);
00196
00197 FDNCTrueE_TrueNC = new NuMatrix1D(*FDNC_TrueNC->ProjectionX(),0);
00198 FDNCTrueE_NuMu = new NuMatrix1D(*FDNC_NuMu->ProjectionX(),0);
00199 FDNCTrueE_BeamNue = new NuMatrix1D(*FDNC_BeamNue->ProjectionX(),0);
00200 FDNCTrueE_AppNue = new NuMatrix1D(*FDNC_AppNue->ProjectionX(),0);
00201 FDNCTrueE_AppNuTau = new NuMatrix1D(*FDNC_AppNuTau->ProjectionX(),0);
00202
00203 FDCCTrueE_TrueNC = new NuMatrix1D(*FDCC_TrueNC->ProjectionX(),0);
00204 FDCCTrueE_NuMu = new NuMatrix1D(*FDCC_NuMu->ProjectionX(),0);
00205 FDCCTrueE_BeamNue = new NuMatrix1D(*FDCC_BeamNue->ProjectionX(),0);
00206 FDCCTrueE_AppNue = new NuMatrix1D(*FDCC_AppNue->ProjectionX(),0);
00207 FDCCTrueE_AppNuTau = new NuMatrix1D(*FDCC_AppNuTau->ProjectionX(),0);
00208
00209 NDNCTrueE_TrueNC = new NuMatrix1D(*NDNC_TrueNC->ProjectionX(),0);
00210 NDNCTrueE_NuMu = new NuMatrix1D(*NDNC_NuMu->ProjectionX(),0);
00211 NDNCTrueE_BeamNue = new NuMatrix1D(*NDNC_BeamNue->ProjectionX(),0);
00212 NDNCTrueE_AppNue = new NuMatrix1D(*NDNC_AppNue->ProjectionX(),0);
00213 NDNCTrueE_AppNuTau = new NuMatrix1D(*NDNC_AppNuTau->ProjectionX(),0);
00214
00215 NDCCTrueE_TrueNC = new NuMatrix1D(*NDCC_TrueNC->ProjectionX(),0);
00216 NDCCTrueE_NuMu = new NuMatrix1D(*NDCC_NuMu->ProjectionX(),0);
00217 NDCCTrueE_BeamNue = new NuMatrix1D(*NDCC_BeamNue->ProjectionX(),0);
00218 NDCCTrueE_AppNue = new NuMatrix1D(*NDCC_AppNue->ProjectionX(),0);
00219 NDCCTrueE_AppNuTau = new NuMatrix1D(*NDCC_AppNuTau->ProjectionX(),0);
00220
00221 got_histos = true;
00222
00223 cout << "got all histos" << endl;
00224 }
00225
00226
00227 Double_t MakeSterilePred_Interpolated::ComparePredWithData(const NuMMParameters &pars)
00228 {
00229 GetHistos();
00230
00231 TH1D *FD_predNC = FDPrediction(pars, NC);
00232 TH1D *FD_predCC = FDPrediction(pars, CC);
00233
00234 Double_t likelihood = 0.0;
00235
00236 likelihood += StatsLikelihood(FD_predNC, FD_dataNC);
00237 likelihood += StatsLikelihood(FD_predCC, FD_dataCC);
00238
00239 delete FD_predNC;
00240 delete FD_predCC;
00241
00242 return likelihood;
00243 }
00244
00245
00246 std::vector<TH1D> MakeSterilePred_Interpolated::WriteFDPredHistos(const NuMMParameters &) const
00247 {
00248 return vector<TH1D>();
00249 }
00250
00251
00252 TH1D *MakeSterilePred_Interpolated::FDPrediction(const NuMMParameters &pars, whichPred Pred)
00253 {
00254 GetHistos();
00255 TH1D *ND_dataE = 0;
00256 switch(Pred)
00257 {
00258 case NC:
00259 ND_dataE = ND_dataNC;
00260 break;
00261 case CC:
00262 ND_dataE = ND_dataCC;
00263 break;
00264 default:
00265 cout << "wrong pred type" << endl;
00266 break;
00267 }
00268
00269 TH1D *FDOsc_MC = MCOscillated(pars, Pred, FD);
00270 TH1D *NDOsc_MC = MCOscillated(pars, Pred, ND);
00271
00272 FDOsc_MC->Divide(NDOsc_MC);
00273 FDOsc_MC->Multiply(ND_dataE);
00274
00275 return FDOsc_MC;
00276 }
00277
00278
00279 TH1D *MakeSterilePred_Interpolated::FDPrediction(NuXMLConfig *xml_config, whichPred Pred)
00280 {
00281 NuMMParameters pars = XMLtoParameters(xml_config);
00282 return FDPrediction(pars, Pred);
00283 }
00284
00285
00286 TH1D *MakeSterilePred_Interpolated::MCOscillated(const NuMMParameters &pars, whichPred Pred, whichDet Det)
00287 {
00288
00289 GetHistos();
00290
00291 cleanup();
00292
00293 double baseline = 0;
00294
00295 TH2D *RecoVsTrue[5];
00296 NuMatrix1D *unoscTrueE[5];
00297 switch(Det*1000+Pred)
00298 {
00299 case ND*1000+NC:
00300 RecoVsTrue[0] = NDNC_TrueNC;
00301 RecoVsTrue[1] = NDNC_NuMu;
00302 RecoVsTrue[2] = NDNC_BeamNue;
00303 RecoVsTrue[3] = NDNC_AppNue;
00304 RecoVsTrue[4] = NDNC_AppNuTau;
00305 unoscTrueE[0] = NDNCTrueE_TrueNC;
00306 unoscTrueE[1] = NDNCTrueE_NuMu;
00307 unoscTrueE[2] = NDNCTrueE_BeamNue;
00308 unoscTrueE[3] = NDNCTrueE_AppNue;
00309 unoscTrueE[4] = NDNCTrueE_AppNuTau;
00310 baseline = ND_baseline;
00311 break;
00312 case ND*1000+CC:
00313 RecoVsTrue[0] = NDCC_TrueNC;
00314 RecoVsTrue[1] = NDCC_NuMu;
00315 RecoVsTrue[2] = NDCC_BeamNue;
00316 RecoVsTrue[3] = NDCC_AppNue;
00317 RecoVsTrue[4] = NDCC_AppNuTau;
00318 unoscTrueE[0] = NDCCTrueE_TrueNC;
00319 unoscTrueE[1] = NDCCTrueE_NuMu;
00320 unoscTrueE[2] = NDCCTrueE_BeamNue;
00321 unoscTrueE[3] = NDCCTrueE_AppNue;
00322 unoscTrueE[4] = NDCCTrueE_AppNuTau;
00323 baseline = ND_baseline;
00324 break;
00325 case FD*1000+NC:
00326 RecoVsTrue[0] = FDNC_TrueNC;
00327 RecoVsTrue[1] = FDNC_NuMu;
00328 RecoVsTrue[2] = FDNC_BeamNue;
00329 RecoVsTrue[3] = FDNC_AppNue;
00330 RecoVsTrue[4] = FDNC_AppNuTau;
00331 unoscTrueE[0] = FDNCTrueE_TrueNC;
00332 unoscTrueE[1] = FDNCTrueE_NuMu;
00333 unoscTrueE[2] = FDNCTrueE_BeamNue;
00334 unoscTrueE[3] = FDNCTrueE_AppNue;
00335 unoscTrueE[4] = FDNCTrueE_AppNuTau;
00336 baseline = FD_baseline;
00337 break;
00338 case FD*1000+CC:
00339 RecoVsTrue[0] = FDCC_TrueNC;
00340 RecoVsTrue[1] = FDCC_NuMu;
00341 RecoVsTrue[2] = FDCC_BeamNue;
00342 RecoVsTrue[3] = FDCC_AppNue;
00343 RecoVsTrue[4] = FDCC_AppNuTau;
00344 unoscTrueE[0] = FDCCTrueE_TrueNC;
00345 unoscTrueE[1] = FDCCTrueE_NuMu;
00346 unoscTrueE[2] = FDCCTrueE_BeamNue;
00347 unoscTrueE[3] = FDCCTrueE_AppNue;
00348 unoscTrueE[4] = FDCCTrueE_AppNuTau;
00349 baseline = FD_baseline;
00350 break;
00351 default:
00352 cout << "what detector and what pred type???" << endl;
00353 break;
00354 }
00355
00356
00357
00358 TH1D *binTemplate = RecoVsTrue[0]->ProjectionY();
00359 binTemplate->Reset();
00360 TH1D *oscRecoE = (TH1D*) binTemplate->Clone("oscRecoE");
00361
00362 oscAvgOverBinNC = (TH1D*) binTemplate->Clone("oscAvgOverBinNC");
00363 oscAvgOverBinNuMu = (TH1D*) binTemplate->Clone("oscAvgOverBinNuMu");
00364 oscIterOverBinNC = new TGraph();
00365 oscIterOverBinNuMu = new TGraph();
00366
00367 for(int kcomp=0; kcomp<5; kcomp++){
00368
00369
00370 double xmin = unoscTrueE[kcomp]->GetXmin();
00371 double xmax = unoscTrueE[kcomp]->GetXmax();
00372
00373
00374
00375
00376
00377 TAxis *Yaxis = RecoVsTrue[kcomp]->GetYaxis();
00378 TAxis *Xaxis = RecoVsTrue[kcomp]->GetXaxis();
00379 for(int x = 1; x <= Xaxis->GetNbins(); x++)
00380 {
00381 double oscWeight = 0.0;
00382 double norm = 0.0;
00383
00384 double bup = Xaxis->GetBinUpEdge(x);
00385 double blow = Xaxis->GetBinLowEdge(x);
00386
00387 if( baseline > 0 && unoscTrueE[kcomp]->Spectrum()->GetBinContent(x) > 0 )
00388 {
00389 if(fOscBinCenter==false){
00390
00391
00392 if(blow < xmin ) blow = xmin;
00393 if(bup > xmax ) bup = xmax;
00394
00395
00396 double bctr = (bup + blow) / 2;
00397 double bwdt = (bup - blow);
00398
00399
00400 double freq = 2 * 1.267 * TMath::Max(pars.Dm243(),pars.Dm2()) * (baseline/Munits::km) / pow(bctr,2);
00401
00402
00403
00404 int n_div = TMath::Ceil( 0.37 * pow(0.05/100,-0.4) * pow(bwdt/bctr,0.8) );
00405
00406
00407 double arg = freq * bwdt/n_div / 2;
00408 double sample = bwdt/n_div / (2 * sqrt(3));
00409 if(arg>0) sample = TMath::ACos( TMath::Sin(arg)/arg ) / freq;
00410
00411 for(int i = 0; i < 2*n_div; i++){
00412
00413
00414 double E = blow + ((i-i%2)/2+0.5)*bwdt/n_div + (2*(i%2)-1)*sample;
00415
00416
00417 double interp = unoscTrueE[kcomp]->GetEventDensity(E);
00418 oscWeight += interp * OscProb(E, pars, baseline, kcomp);
00419 norm += interp;
00420
00421 }
00422
00423
00424 if(norm>0) oscWeight /= norm;
00425
00426 }
00427
00428 else{
00429 double E = Xaxis->GetBinCenter(x);
00430 oscWeight = OscProb(E, pars, baseline, kcomp);
00431 }
00432
00433 }
00434 else
00435 {
00436 if( kcomp!=1 && kcomp!=2 ) oscWeight = 0.0;
00437 else oscWeight = 1.0;
00438 }
00439
00440 if(kcomp==0) oscAvgOverBinNC->SetBinContent(x, oscWeight);
00441 else if(kcomp==1) oscAvgOverBinNuMu->SetBinContent(x, oscWeight);
00442
00443
00444 for(int y = 1; y <= Yaxis->GetNbins(); y++){
00445 if(kcomp==0) oscRecoE->AddBinContent(y, RecoVsTrue[kcomp]->GetBinContent(x,y)*(1.0-oscWeight));
00446 else oscRecoE->AddBinContent(y, RecoVsTrue[kcomp]->GetBinContent(x,y)*oscWeight);
00447 }
00448
00449 }
00450
00451 }
00452
00453 if( fFitSysts ){
00454
00455
00456 if(Det==FD){
00457 if(Pred==NC) oscRecoE->Scale(pars.GetParameterByIndex(NuMMParameters::kNormNC));
00458 else if(Pred==CC) oscRecoE->Scale(pars.GetParameterByIndex(NuMMParameters::kNorm));
00459 }
00460
00461
00462 for(int kPar=0; kPar<kNumSyst; kPar++){
00463
00464
00465 if(kPar==kShwScale || kPar==kTrkScale || kPar==kNcBack){ if(Pred==NC) continue; }
00466 else{
00467 if(Pred==CC) continue;
00468 else{
00469 if(kPar!=kAbsShwNC && kPar!=kCCBack && kPar!=kNDClean && Det==ND) continue;
00470 if(kPar==kNDClean && Det==FD) continue;
00471 }
00472 }
00473
00474
00475 double parvalue = pars.GetParameterByIndex(fParMap.at(kPar));
00476 if(kPar==kNcBack) parvalue = (parvalue-1.0)/0.2;
00477
00478
00479 if(parvalue==0) continue;
00480
00481
00482 TH1D *ScaledSyst;
00483 if(parvalue<0)
00484 ScaledSyst = (TH1D*)fSystMinus[Det][kPar]->Clone();
00485 else
00486 ScaledSyst = (TH1D*)fSystPlus[Det][kPar]->Clone();
00487
00488 ScaledSyst->Scale(TMath::Abs(parvalue));
00489 ScaledSyst->Add(fhOne);
00490 oscRecoE->Multiply(ScaledSyst);
00491
00492 }
00493
00494 }
00495
00496 return oscRecoE;
00497 }
00498
00499
00500 TH1D *MakeSterilePred_Interpolated::MCOscillated(NuXMLConfig *xml_config, whichPred Pred, whichDet Det)
00501 {
00502 NuMMParameters pars = XMLtoParameters(xml_config);
00503 return MCOscillated(pars, Pred, Det);
00504 }
00505
00506
00507 NuMMParameters MakeSterilePred_Interpolated::XMLtoParameters(NuXMLConfig *xml_config)
00508 {
00509 NuMMParameters pars;
00510
00511 pars.Dm2( xml_config->DM2Nu());
00512 pars.Theta23( xml_config->Theta23());
00513 pars.Dm221( xml_config->DM221());
00514 pars.Dm243( xml_config->DM243());
00515 pars.Delta1( xml_config->Delta1());
00516 pars.Delta2( xml_config->Delta2());
00517 pars.Delta3( xml_config->Delta3());
00518 pars.Theta12( xml_config->Theta12());
00519 pars.Theta13( xml_config->Theta13());
00520 pars.Theta14( xml_config->Theta14());
00521 pars.Theta24( xml_config->Theta24());
00522 pars.Theta34( xml_config->Theta34());
00523 return pars;
00524 }
00525
00526
00527 void MakeSterilePred_Interpolated::SetOscBinCenter(bool IsOscBinCenter)
00528 {
00529 fOscBinCenter = IsOscBinCenter;
00530 }
00531
00532
00533 void MakeSterilePred_Interpolated::SetSystShifts(TFile *systFile)
00534 {
00535
00536 string Pred[2] = {"CC","NC"};
00537 string Det[2] = {"ND","FD"};
00538 string SystName[kNumSyst];
00539 SystName[kShwScale] = "absHadCC";
00540 SystName[kTrkScale] = "absLep";
00541 SystName[kNcBack] = "ncBkg";
00542 SystName[kAbsShwNC] = "absHad";
00543 SystName[kRelShwNC] = "relativeHadronic";
00544 SystName[kCCBack] = "ccBkg";
00545 SystName[kFDClean] = "fdCleaning";
00546 SystName[kNDClean] = "ndCleaning";
00547 SystName[kFDCosmic] = "fdCosmics";
00548
00549 assert(FD_data);
00550 TH1D *binTemplate;
00551 FD_data->GetObject("hRecoEnergyNC",binTemplate);
00552
00553 assert(systFile);
00554
00555 for(int kDet=0; kDet<2; kDet++){
00556 for(int kPar=0; kPar<kNumSyst; kPar++){
00557
00558 int kPred = 1;
00559 if(kPar==kShwScale || kPar==kTrkScale || kPar==kNcBack) kPred=0;
00560
00561 systFile->GetObject((Det[kDet]+Pred[kPred]+SystName[kPar]+"Plus").c_str(),fSystPlus[kDet][kPar]);
00562 systFile->GetObject((Det[kDet]+Pred[kPred]+SystName[kPar]+"Minus").c_str(),fSystMinus[kDet][kPar]);
00563
00564 fSystPlus[kDet][kPar] = RebinSystShift(fSystPlus[kDet][kPar],binTemplate);
00565 fSystMinus[kDet][kPar] = RebinSystShift(fSystMinus[kDet][kPar],binTemplate);
00566
00567 }}
00568
00569 fhOne = (TH1D*)binTemplate->Clone();
00570 for(int i=0; i<=fhOne->GetNbinsX()+1; i++){
00571 fhOne->SetBinContent(i,1);
00572 fhOne->SetBinError(i,0);
00573 }
00574
00575 fFitSysts = true;
00576
00577 }
00578
00579
00580 TH1D* MakeSterilePred_Interpolated::RebinSystShift(TH1D*& SystHist,TH1D* binTemplate) const
00581 {
00582
00583 TH1D *outhist = (TH1D*)binTemplate->Clone();
00584
00585 for(int i=0;i<=outhist->GetNbinsX()+1;i++){
00586 double xval = outhist->GetBinCenter(i);
00587 double yval = SystHist->Interpolate(xval);
00588 outhist->SetBinContent(i,yval);
00589 outhist->SetBinError(i,0);
00590 }
00591
00592 return outhist;
00593
00594 }
00595
00596
00597 Double_t MakeSterilePred_Interpolated::OscProb(double E, const NuMMParameters &pars, double baseline, int component)
00598 {
00599
00600 double prob;
00601
00602 switch(component){
00603
00604 case 1: prob = NuOscProbCalc::FourFlavourDisappearanceWeight(E, pars.Dm2(), pars.Theta23(), pars.Dm221(),
00605 pars.Dm243(), pars.Theta12(), pars.Theta13(), pars.Theta14(),
00606 pars.Theta24(), pars.Theta34(), pars.Delta1(), pars.Delta2(),
00607 pars.Delta3(), baseline); break;
00608
00609 case 2: prob = NuOscProbCalc::FourFlavourNuESurvivalProbability(E, pars.Dm2(), pars.Theta23(), pars.Dm221(),
00610 pars.Dm243(), pars.Theta12(), pars.Theta13(), pars.Theta14(),
00611 pars.Theta24(), pars.Theta34(), pars.Delta1(), pars.Delta2(),
00612 pars.Delta3(), baseline); break;
00613
00614 case 3: prob = NuOscProbCalc::FourFlavourNuMuToNuEProbability(E, pars.Dm2(), pars.Theta23(), pars.Dm221(),
00615 pars.Dm243(), pars.Theta12(), pars.Theta13(), pars.Theta14(),
00616 pars.Theta24(), pars.Theta34(), pars.Delta1(), pars.Delta2(),
00617 pars.Delta3(), baseline); break;
00618
00619 case 4: prob = NuOscProbCalc::FourFlavourNuMuToNuTauProbability(E, pars.Dm2(), pars.Theta23(), pars.Dm221(),
00620 pars.Dm243(), pars.Theta12(), pars.Theta13(), pars.Theta14(),
00621 pars.Theta24(), pars.Theta34(), pars.Delta1(), pars.Delta2(),
00622 pars.Delta3(), baseline); break;
00623
00624 default: prob = NuOscProbCalc::FourFlavourNuMuToNuSProbability(E, pars.Dm2(), pars.Theta23(), pars.Dm221(),
00625 pars.Dm243(), pars.Theta12(), pars.Theta13(), pars.Theta14(),
00626 pars.Theta24(), pars.Theta34(), pars.Delta1(), pars.Delta2(),
00627 pars.Delta3(), baseline); break;
00628
00629 }
00630
00631 return prob;
00632
00633 }
00634