loon: MakeSterilePred_Interpolated.cxx Source File

00001 // a class that takes reco vs. true matrices and oscillates the reco energy based on an oscillation weight
00002 // outputs a total oscillated reco spectrum for certain oscillation parameters and baseline
00003 // August 17, 2012: by Alena Devan
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         // GetObjects
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         // POT normalize
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 &/*pars*/) 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         // clone the TrueNC projection plot so the new oscillated spectra have the same binning
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           // Get the min and max energies from the true energy interpolation
00370           double xmin = unoscTrueE[kcomp]->GetXmin();
00371           double xmax = unoscTrueE[kcomp]->GetXmax();
00372 
00373           // loop over every true energy bin in the reco vs. true matrices, then loop over every reco energy in that bin
00374           // to calculate an oscillation weight for that reco energy based on the true energy. Get the 4 flavor oscillations
00375           // from OscProb/NuOscProbCalc. Fill a new histogram for each of the 5 types of events with the now oscillated 
00376           // reco energy. 
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                           // Remove empty intervals from boundary bins
00392                           if(blow < xmin ) blow = xmin;
00393                           if(bup > xmax ) bup = xmax;
00394 
00395                           // Calculate effective bin center and width
00396                           double bctr = (bup + blow) / 2;
00397                           double bwdt = (bup - blow);
00398 
00399                           // Calculate max frequency of oscillations
00400                           double freq = 2 * 1.267 * TMath::Max(pars.Dm243(),pars.Dm2()) * (baseline/Munits::km) / pow(bctr,2);
00401 
00402                           // Set number of subdivisions to get approx. 0.05% precision
00403                           // This is usually an overestimate.
00404                           int n_div = TMath::Ceil( 0.37 * pow(0.05/100,-0.4) * pow(bwdt/bctr,0.8) );
00405 
00406                           // Find optimal sampling length within subdivisions
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                             // Sample in center +/- the sampling length
00414                             double E = blow + ((i-i%2)/2+0.5)*bwdt/n_div + (2*(i%2)-1)*sample;
00415 
00416                             // Oscillate weighted by event density
00417                             double interp = unoscTrueE[kcomp]->GetEventDensity(E);
00418                             oscWeight += interp * OscProb(E, pars, baseline, kcomp);
00419                             norm += interp;
00420                         
00421                           }// Loop over bin subdivisions
00422 
00423                           // Divide by sum of weights to get average
00424                           if(norm>0) oscWeight /= norm;
00425 
00426                         }// If averaging over bin
00427 
00428                         else{
00429                           double E = Xaxis->GetBinCenter(x);
00430                           oscWeight = OscProb(E, pars, baseline, kcomp);
00431                         }// If oscillate at bin center
00432 
00433                 }// If there is something to oscillate
00434                 else
00435                 {
00436                         if( kcomp!=1 && kcomp!=2 ) oscWeight = 0.0; // Appearance defaults 0.0
00437                         else                       oscWeight = 1.0; // Disappearance defaults 1.0
00438                 }// If there is nothing to oscillate
00439 
00440                 if(kcomp==0)      oscAvgOverBinNC->SetBinContent(x, oscWeight);
00441                 else if(kcomp==1) oscAvgOverBinNuMu->SetBinContent(x, oscWeight);
00442                                                                 
00443                 // using the oscillation weight, fill a 1d histogram for each type of event with the oscillated reco energy
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                 }//Loop over Reco Energy 
00448 
00449           }//Loop over True Energy
00450 
00451         }//Loop over components
00452 
00453         if( fFitSysts ){
00454 
00455           // Apply normalization systematics
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           // Apply energy dependent systematics
00462           for(int kPar=0; kPar<kNumSyst; kPar++){
00463 
00464             // Only apply shifts to spectra that do change
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             // Get the shift value in sigmas
00475             double parvalue = pars.GetParameterByIndex(fParMap.at(kPar));
00476             if(kPar==kNcBack) parvalue = (parvalue-1.0)/0.2;
00477 
00478             // Skip parameters that do nothing
00479             if(parvalue==0) continue;
00480 
00481             // Get the shift histogram by sign
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)); // Scale to shift value
00489             ScaledSyst->Add(fhOne);                  // Subtract one
00490             oscRecoE->Multiply(ScaledSyst);         // Apply systematic
00491 
00492           }// Loop over systematics
00493 
00494         }// Apply systematics if present
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         // set pars to xml_config
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 // vi:noet:ts=4:sw=4:sts=0