Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages
NtupleUtils

NuMMRunCC2010.cxx

Go to the documentation of this file.
00001 #include "TH1.h"
00002 
00003 #include "MessageService/MsgService.h"
00004 #include "NtupleUtils/NuHistInterpolator.h"
00005 #include "NtupleUtils/NuMMHelperPRL.h"
00006 #include "NtupleUtils/NuMMHelperCPT.h"
00007 #include "NtupleUtils/NuMMParameters.h"
00008 #include "NtupleUtils/NuMMRunCC2010.h"
00009 
00010 ClassImp(NuMMRunCC2010)
00011 
00012 CVSID("$Id: NuMMRunCC2010.cxx,v 1.10 2010/02/03 14:41:01 bckhouse Exp $");
00013 
00014 //____________________________________________________________________72
00015 NuMMRunCC2010::NuMMRunCC2010(NuMMHelperPRL* helper,
00016                              NuMatrixSpectrum* ndCCData,
00017                              NuMatrixSpectrum* fdCCData)
00018   : NuMMRun(),
00019     fNDCCData(ndCCData),
00020     fFDCCData(fdCCData),
00021     fHelper(helper),
00022     fHelperCPT(0)
00023 {
00024   //Independent of fit parameters for a stats-only fit
00025   fFDFlux = this->CalculateFDFlux();
00026 }
00027 
00028 
00029 //____________________________________________________________________72
00030 NuMMRunCC2010::NuMMRunCC2010(NuMMHelperCPT* helper,
00031                              NuMatrixSpectrum* ndCCData,
00032                              NuMatrixSpectrum* ndCCDataBar,
00033                              NuMatrixSpectrum* fdCCData,
00034                              NuMatrixSpectrum* fdCCDataBar)
00035   : NuMMRun(),
00036     fNDCCData(ndCCData),
00037     fFDCCData(fdCCData),
00038     fNDCCDataBar(ndCCDataBar),
00039     fFDCCDataBar(fdCCDataBar),
00040     fHelper(0),
00041     fHelperCPT(helper)
00042 {
00043   fFDFlux = this->CalculateFDFlux();
00044 }
00045 
00046 //____________________________________________________________________72
00047 const std::vector<TH1D> NuMMRunCC2010
00048 ::WriteFDPredHistos(const NuMMParameters& pars) const
00049 {
00050   vector<TH1D> vHistos;
00051 
00052   pair<NuMatrixSpectrum,NuMatrixSpectrum> predictions = fFDFlux;  
00053   NuMatrixSpectrum nuPrediction = fFDFlux.first;
00054   NuMatrixSpectrum barPrediction = fFDFlux.first;
00055   
00056   if(fHelper){
00057     TH1D hNDData(*(fNDCCData->Spectrum()));
00058     hNDData.SetName("ndData");
00059     hNDData.SetTitle("ndData");
00060     vHistos.push_back(*(new TH1D(hNDData)));
00061   
00062     NuMatrixSpectrum potentialTaus(nuPrediction);
00063     potentialTaus.Multiply(fHelper->XSecGraphTaus());
00064     potentialTaus.Multiply(fHelper->FDTauEfficiency());
00065     
00066     nuPrediction.Multiply(fHelper->XSecGraph());
00067     nuPrediction.Multiply(fHelper->FDEfficiency());
00068     
00069     NuMatrixSpectrum unoscTrueSpectrum(nuPrediction);
00070     unoscTrueSpectrum.TrueToReco(fHelper->FDRecoVsTrue());
00071     NuMatrixSpectrum ncBackground(unoscTrueSpectrum);
00072     ncBackground.Divide(fHelper->FDPurity());
00073     ncBackground.Subtract(unoscTrueSpectrum);
00074 
00075     OscillateNuMu(nuPrediction, pars.Dm2(), pars.Sn2());
00076     InverseOscillateNuTau(potentialTaus, pars.Dm2(), pars.Sn2());
00077     
00078     //Output the prediction in true energy
00079     TH1D hTrueFDPrediction(*(nuPrediction.Spectrum()));
00080     hTrueFDPrediction.SetName("fdTruePrediction");
00081     vHistos.push_back(hTrueFDPrediction);
00082     
00083     potentialTaus.TrueToReco(fHelper->FDTauRecoVsTrue());  
00084     nuPrediction.TrueToReco(fHelper->FDRecoVsTrue());
00085     
00086     //Get the prediction in reco energy
00087     TH1D hBaseFDPrediction(*(nuPrediction.Spectrum()));
00088     hBaseFDPrediction.SetName("fdBasePrediction");
00089     vHistos.push_back(hBaseFDPrediction);
00090     
00091     nuPrediction.Add(ncBackground);
00092     if (0==fDisappearanceModel){
00093       nuPrediction.Add(potentialTaus);
00094     }
00095     
00096     //Get the prediction with backgrounds
00097     TH1D hFDPrediction(*(nuPrediction.Spectrum()));  
00098     hFDPrediction.SetName("fdPrediction");
00099     vHistos.push_back(hFDPrediction);
00100     
00101     TH1D hFDTaus(*(potentialTaus.Spectrum()));
00102     hFDTaus.SetName("fdTaus");
00103     vHistos.push_back(hFDTaus);
00104     
00105     TH1D hFDNCBackground(*(ncBackground.Spectrum()));
00106     hFDNCBackground.SetName("fdNCBackground");
00107     vHistos.push_back(hFDNCBackground);
00108     
00109     TH1D hfdData(*(fFDCCData->Spectrum()));
00110     hfdData.SetName("fdData");
00111     vHistos.push_back(hfdData);
00112   }
00113   
00114   if(fHelperCPT){
00115     //Put the ND data in the vector:
00116     TH1D hNDNuData(*(fNDCCData->Spectrum()));
00117     hNDNuData.SetName("ndDataNQ");
00118     hNDNuData.SetTitle("ndDataNQ");
00119     vHistos.push_back(*(new TH1D(hNDNuData)));
00120     TH1D hNDBarData(*(fNDCCDataBar->Spectrum()));
00121     hNDBarData.SetName("ndDataPQ");
00122     hNDBarData.SetTitle("ndDataPQ");
00123     vHistos.push_back(*(new TH1D(hNDBarData)));
00124     
00125     cout << "sn2: " << pars.Sn2() << "; dm2: " << pars.Dm2() 
00126          << "; sn2bar: " << pars.Sn2Bar() << "; dm2bar: " << pars.Dm2Bar()
00127          << endl
00128          << "; norm: " << pars.Normalisation() 
00129          << "; NCBack: " << pars.NCBackgroundScale()
00130          << "; ShwEn: " << pars.ShwEnScale() << endl;
00131     
00132     //Get the neutrinos to the FD
00133     NuMatrixSpectrum nuPrediction(*fNDCCData);
00134     nuPrediction.Multiply(fHelperCPT->NDNuPurity());
00135     nuPrediction.RecoToTrue(fHelperCPT->NDNuRecoVsTrue());
00136     nuPrediction.Divide(fHelperCPT->NDNuEfficiency());
00137     nuPrediction.Divide(fHelperCPT->NuXSecGraph());
00138     nuPrediction.Divide(fNDCCData->PoT());
00139     nuPrediction.Divide(fNDFidMass);
00140     nuPrediction.ExtrapolateNDToFD(fHelperCPT->NuBeamMatrix());
00141     nuPrediction.Multiply(fFDFidMass);
00142     nuPrediction.Multiply(fFDCCData->PoT());
00143     nuPrediction.ResetPoT(fFDCCData->PoT());
00144     
00145     //Get the antineutrinos to the FD
00146     NuMatrixSpectrum barPrediction(*fNDCCDataBar);
00147     barPrediction.Multiply(fHelperCPT->NDBarPurity());
00148     barPrediction.RecoToTrue(fHelperCPT->NDBarRecoVsTrue());
00149     barPrediction.Divide(fHelperCPT->NDBarEfficiency());
00150     barPrediction.Divide(fHelperCPT->BarXSecGraph());
00151     barPrediction.Divide(fNDCCDataBar->PoT());
00152     barPrediction.Divide(fNDFidMass);
00153     barPrediction.ExtrapolateNDToFD(fHelperCPT->BarBeamMatrix());
00154     barPrediction.Multiply(fFDFidMass);
00155     barPrediction.Multiply(fFDCCDataBar->PoT());
00156     barPrediction.ResetPoT(fFDCCDataBar->PoT());
00157 
00158     //Get the taus (ignoring any wrong-sign taubars)
00159     NuMatrixSpectrum potentialNuTaus(nuPrediction);
00160     potentialNuTaus.Multiply(fHelperCPT->XSecGraphNuTaus());
00161     potentialNuTaus.Multiply(fHelperCPT->FDNuTauEfficiency());
00162     InverseOscillateNuTau(potentialNuTaus, pars.Dm2(), pars.Sn2());
00163     //Get the taubars (ignoring any wrong-sign taus)
00164     NuMatrixSpectrum potentialTauBars(barPrediction);
00165     potentialTauBars.Multiply(fHelperCPT->XSecGraphTauBars());
00166     potentialTauBars.Multiply(fHelperCPT->FDTauBarEfficiency());
00167     InverseOscillateNuTau(potentialTauBars, pars.Dm2(), pars.Sn2());
00168     
00169     // Write out the true tau spectra
00170     TH1D hFDTausNuMu(*(potentialNuTaus.Spectrum()));
00171     hFDTausNuMu.SetName("fdTausNuMu");
00172     vHistos.push_back(hFDTausNuMu);
00173     TH1D hFDTausNuBar(*(potentialTauBars.Spectrum()));
00174     hFDTausNuBar.SetName("fdTausNuBar");
00175     vHistos.push_back(hFDTausNuBar);
00176     
00177     // Do the tau true to reco speerately so that we can write out the spectra
00178     potentialNuTaus.TrueToReco(fHelperCPT->FDNuTauRecoVsTrue());
00179     potentialTauBars.TrueToReco(fHelperCPT->FDTauBarRecoVsTrue());
00180     
00181     //Cross sections
00182     nuPrediction.Multiply(fHelperCPT->NuXSecGraph());
00183     barPrediction.Multiply(fHelperCPT->BarXSecGraph());
00184     
00185     //Get the neutrino NC background
00186     NuMatrixSpectrum nuNCBackground(nuPrediction);
00187     nuNCBackground.Multiply(fHelperCPT->FDNuEfficiency());
00188     nuNCBackground.TrueToReco(fHelperCPT->FDNuRecoVsTrue());
00189     nuNCBackground.Divide(fHelperCPT->FDNuPurity());
00190     nuNCBackground.Multiply(fHelperCPT->FDNuNCContamination());
00191     
00192     //Get the antineutrino NC background
00193     NuMatrixSpectrum barNCBackground(barPrediction);
00194     barNCBackground.Multiply(fHelperCPT->FDBarEfficiency());
00195     barNCBackground.TrueToReco(fHelperCPT->FDBarRecoVsTrue());
00196     barNCBackground.Divide(fHelperCPT->FDBarPurity());
00197     barNCBackground.Multiply(fHelperCPT->FDBarNCContamination());
00198     
00199     //Get the wrong-sign neutrino background for the antineutrino
00200     //prediction
00201     NuMatrixSpectrum wrongSignNuMus(nuPrediction);
00202     wrongSignNuMus.Multiply(fHelperCPT->FDWrongSignNuEfficiency());
00203     OscillateNuMu(wrongSignNuMus, pars.Dm2(), pars.Sn2());
00204     
00205     //Get the wrong-sign antineutrino background for the neutrino
00206     //prediction
00207     NuMatrixSpectrum wrongSignNuBars(barPrediction);
00208     wrongSignNuBars.Multiply(fHelperCPT->FDWrongSignBarEfficiency());
00209     OscillateNuMu(wrongSignNuBars, pars.Dm2(), pars.Sn2());
00210     
00211     TH1D hFDWrongSignNQ(*(wrongSignNuBars.Spectrum()));
00212     hFDWrongSignNQ.SetName("fdWrongSignNQ");
00213     vHistos.push_back(hFDWrongSignNQ);
00214     
00215     TH1D hFDWrongSignPQ(*(wrongSignNuMus.Spectrum()));
00216     hFDWrongSignPQ.SetName("fdWrongSignPQ");
00217     vHistos.push_back(hFDWrongSignPQ);
00218     
00219     // Do the truetoreco seperately so that we can write out the histos first
00220     wrongSignNuBars.TrueToReco(fHelperCPT->FDWrongSignBarRecoVsTrue());
00221     wrongSignNuMus.TrueToReco(fHelperCPT->FDWrongSignNuRecoVsTrue());
00222     
00223     //Oscillations
00224     OscillateNuMu(nuPrediction, pars.Dm2(), pars.Sn2());
00225     OscillateNuMu(barPrediction, pars.Dm2(), pars.Sn2());
00226     //Efficiencies
00227     nuPrediction.Multiply(fHelperCPT->FDNuEfficiency());
00228     barPrediction.Multiply(fHelperCPT->FDBarEfficiency());
00229     
00230     // Write out the true spectra for the prediction components
00231     TH1D hFDNoBackTrueNQ(*(nuPrediction.Spectrum()));
00232     hFDNoBackTrueNQ.SetName("fdBasePredictionNuMus");
00233     vHistos.push_back(hFDNoBackTrueNQ);
00234     
00235     TH1D hFDNoBackTruePQ(*(barPrediction.Spectrum()));
00236     hFDNoBackTruePQ.SetName("fdBasePredictionNuBars");
00237     vHistos.push_back(hFDNoBackTruePQ);
00238     
00239     //True to reco
00240     nuPrediction.TrueToReco(fHelperCPT->FDNuRecoVsTrue());
00241     barPrediction.TrueToReco(fHelperCPT->FDBarRecoVsTrue());
00242     
00243     // Dump out the prediction spectrum before adding in the backgrounds
00244     TH1D hFDNoBackNQ(*(nuPrediction.Spectrum()));
00245     hFDNoBackNQ.SetName("fdBasePredictionNQ");
00246     vHistos.push_back(hFDNoBackNQ);
00247     
00248     TH1D hFDNoBackPQ(*(barPrediction.Spectrum()));
00249     hFDNoBackPQ.SetName("fdBasePredictionPQ");
00250     vHistos.push_back(hFDNoBackPQ);
00251     
00252     //Add in backgrounds
00253     nuPrediction.Add(wrongSignNuBars);
00254     barPrediction.Add(wrongSignNuMus);
00255     nuPrediction.Add(nuNCBackground);
00256     nuPrediction.Add(potentialNuTaus);
00257     barPrediction.Add(barNCBackground);
00258     barPrediction.Add(potentialTauBars);
00259     
00260     //Put the components of the predictions in the vector
00261     
00262     TH1D hFDPredictionNQ(*(nuPrediction.Spectrum()));
00263     hFDPredictionNQ.SetName("fdPredictionNQ");
00264     vHistos.push_back(hFDPredictionNQ);
00265     
00266     TH1D hFDPredictionPQ(*(barPrediction.Spectrum()));
00267     hFDPredictionPQ.SetName("fdPredictionPQ");
00268     vHistos.push_back(hFDPredictionPQ);
00269     
00270     TH1D hFDNCNQ(*(nuNCBackground.Spectrum()));
00271     hFDNCNQ.SetName("fdNCNQ");
00272     vHistos.push_back(hFDNCNQ);
00273     
00274     TH1D hFDNCPQ(*(barNCBackground.Spectrum()));
00275     hFDNCPQ.SetName("fdNCPQ");
00276     vHistos.push_back(hFDNCPQ);
00277     
00278     TH1D hFDTausNQ(*(potentialNuTaus.Spectrum()));
00279     hFDTausNQ.SetName("fdTausNQ");
00280     vHistos.push_back(hFDTausNQ);
00281     
00282     TH1D hFDTausPQ(*(potentialTauBars.Spectrum()));
00283     hFDTausPQ.SetName("fdTausPQ");
00284     vHistos.push_back(hFDTausPQ);
00285     
00286     TH1D hFDDataNQ(*(fFDCCData->Spectrum()));
00287     hFDDataNQ.SetName("fdDataNQ");
00288     vHistos.push_back(hFDDataNQ);
00289     
00290     TH1D hFDDataPQ(*(fFDCCDataBar->Spectrum()));
00291     hFDDataPQ.SetName("fdDataPQ");
00292     vHistos.push_back(hFDDataPQ);
00293   }  
00294   else{
00295     MSG("NuMMRunCC2010",Msg::kError)<<"No Helper Set!"<<endl;
00296   }
00297   return vHistos;
00298 }
00299 
00300 //____________________________________________________________________72
00301 const pair<NuMatrixSpectrum,NuMatrixSpectrum>
00302 NuMMRunCC2010::CalculateFDFlux() const
00303 {
00304   if(fHelper){
00305     NuMatrixSpectrum prediction(*fNDCCData);
00306     prediction.Multiply(fHelper->NDPurity());
00307     prediction.RecoToTrue(fHelper->NDRecoVsTrue());
00308     prediction.Divide(fHelper->NDEfficiency());
00309     prediction.Divide(fHelper->XSecGraph());
00310     prediction.Divide(fNDCCData->PoT());
00311     prediction.Divide(fNDFidMass);
00312     prediction.ExtrapolateNDToFD(fHelper->BeamMatrix());
00313     prediction.Multiply(fFDFidMass);
00314     prediction.Multiply(fFDCCData->PoT());
00315     
00316     NuMatrixSpectrum dummy;
00317     pair<NuMatrixSpectrum,NuMatrixSpectrum>
00318       predictions(prediction,dummy);
00319     return predictions;
00320 
00321   }
00322   else if(fHelperCPT){
00323     
00324     //Get the neutrinos to the FD
00325     NuMatrixSpectrum nuPrediction(*fNDCCData);
00326     nuPrediction.Multiply(fHelperCPT->NDNuPurity());
00327     nuPrediction.RecoToTrue(fHelperCPT->NDNuRecoVsTrue());
00328     nuPrediction.Divide(fHelperCPT->NDNuEfficiency());
00329     nuPrediction.Divide(fHelperCPT->NuXSecGraph());
00330     nuPrediction.Divide(fNDCCData->PoT());
00331     nuPrediction.Divide(fNDFidMass);
00332     nuPrediction.ExtrapolateNDToFD(fHelperCPT->NuBeamMatrix());
00333     nuPrediction.Multiply(fFDFidMass);
00334     nuPrediction.Multiply(fFDCCData->PoT());
00335     nuPrediction.ResetPoT(fFDCCData->PoT());
00336 
00337     //Get the antineutrinos to the FD
00338     NuMatrixSpectrum barPrediction(*fNDCCDataBar);
00339     barPrediction.Multiply(fHelperCPT->NDBarPurity());
00340     barPrediction.RecoToTrue(fHelperCPT->NDBarRecoVsTrue());
00341     barPrediction.Divide(fHelperCPT->NDBarEfficiency());
00342     barPrediction.Divide(fHelperCPT->BarXSecGraph());
00343     barPrediction.Divide(fNDCCDataBar->PoT());
00344     barPrediction.Divide(fNDFidMass);
00345     barPrediction.ExtrapolateNDToFD(fHelperCPT->BarBeamMatrix());
00346     barPrediction.Multiply(fFDFidMass);
00347     barPrediction.Multiply(fFDCCDataBar->PoT());
00348     barPrediction.ResetPoT(fFDCCDataBar->PoT());
00349     
00350     pair<NuMatrixSpectrum,NuMatrixSpectrum>
00351       predictions(nuPrediction,barPrediction);
00352     return predictions;
00353   }
00354   else{
00355     //Woah, how did you get here without a helper?
00356     MSG("NuMMRunCC2010",Msg::kError) <<"Helpers not set!"<<endl;
00357     NuMatrixSpectrum dummy;
00358     pair<NuMatrixSpectrum,NuMatrixSpectrum>
00359       predictions(dummy,dummy);
00360     return predictions;
00361   }
00362 }
00363 //____________________________________________________________________72
00364 const pair<NuMatrixSpectrum,NuMatrixSpectrum>
00365 NuMMRunCC2010::MakeFDPred(const NuMMParameters& pars)
00366 {
00367   if (!fQuietMode) {
00368     cout << "sn2: " << pars.Sn2() << "; dm2: " << pars.Dm2() 
00369          << "; sn2bar: " << pars.Sn2Bar() << "; dm2bar: " << pars.Dm2Bar()
00370          << endl
00371          << "; norm: " << pars.Normalisation() 
00372          << "; NCBack: " << pars.NCBackgroundScale()
00373          << "; ShwEn: " << pars.ShwEnScale() << endl;
00374   }
00375   
00376   const pair<NuMatrixSpectrum,NuMatrixSpectrum> predictions =
00377     fFDFlux;
00378   NuMatrixSpectrum nuPrediction = predictions.first;
00379   NuMatrixSpectrum barPrediction = predictions.second;
00380   
00381   if(fHelperCPT){
00382     //Get the taus (ignoring any wrong-sign taubars)
00383     NuMatrixSpectrum potentialNuTaus(nuPrediction);
00384     potentialNuTaus.Multiply(fHelperCPT->XSecGraphNuTaus());
00385     potentialNuTaus.Multiply(fHelperCPT->FDNuTauEfficiency());
00386     InverseOscillateNuTau(potentialNuTaus, pars.Dm2(), pars.Sn2());
00387     potentialNuTaus.TrueToReco(fHelperCPT->FDNuTauRecoVsTrue());
00388     //Get the taubars (ignoring any wrong-sign taus)
00389     NuMatrixSpectrum potentialTauBars(barPrediction);
00390     potentialTauBars.Multiply(fHelperCPT->XSecGraphTauBars());
00391     potentialTauBars.Multiply(fHelperCPT->FDTauBarEfficiency());
00392     InverseOscillateNuTau(potentialTauBars, pars.Dm2(), pars.Sn2());
00393     potentialTauBars.TrueToReco(fHelperCPT->FDTauBarRecoVsTrue());
00394     
00395     //Cross sections
00396     nuPrediction.Multiply(fHelperCPT->NuXSecGraph());
00397     barPrediction.Multiply(fHelperCPT->BarXSecGraph());
00398     
00399     //Get the neutrino NC background
00400     NuMatrixSpectrum nuNCBackground(nuPrediction);
00401     nuNCBackground.Multiply(fHelperCPT->FDNuEfficiency());
00402     nuNCBackground.TrueToReco(fHelperCPT->FDNuRecoVsTrue());
00403     nuNCBackground.Divide(fHelperCPT->FDNuPurity());
00404     nuNCBackground.Multiply(fHelperCPT->FDNuNCContamination());
00405     
00406     //Get the antineutrino NC background
00407     NuMatrixSpectrum barNCBackground(barPrediction);
00408     barNCBackground.Multiply(fHelperCPT->FDBarEfficiency());
00409     barNCBackground.TrueToReco(fHelperCPT->FDBarRecoVsTrue());
00410     barNCBackground.Divide(fHelperCPT->FDBarPurity());
00411     barNCBackground.Multiply(fHelperCPT->FDBarNCContamination());
00412     
00413     //Get the wrong-sign neutrino background for the antineutrino
00414     //prediction
00415     NuMatrixSpectrum wrongSignNuMus(nuPrediction);
00416     wrongSignNuMus.Multiply(fHelperCPT->FDWrongSignNuEfficiency());
00417     OscillateNuMu(wrongSignNuMus, pars.Dm2(), pars.Sn2());
00418     wrongSignNuMus.TrueToReco(fHelperCPT->FDWrongSignNuRecoVsTrue());
00419     
00420     //Get the wrong-sign antineutrino background for the neutrino
00421     //prediction
00422     NuMatrixSpectrum wrongSignNuBars(barPrediction);
00423     wrongSignNuBars.Multiply(fHelperCPT->FDWrongSignBarEfficiency());
00424     OscillateNuMu(wrongSignNuBars, pars.Dm2(), pars.Sn2());
00425     wrongSignNuBars.TrueToReco(fHelperCPT->FDWrongSignBarRecoVsTrue());
00426     
00427     //Oscillations - both sets oscillate together
00428     OscillateNuMu(nuPrediction, pars.Dm2(), pars.Sn2());
00429     OscillateNuMu(barPrediction, pars.Dm2(), pars.Sn2());
00430     //Efficiencies
00431     nuPrediction.Multiply(fHelperCPT->FDNuEfficiency());
00432     barPrediction.Multiply(fHelperCPT->FDBarEfficiency());
00433     
00434     //True to reco
00435     nuPrediction.TrueToReco(fHelperCPT->FDNuRecoVsTrue());
00436     barPrediction.TrueToReco(fHelperCPT->FDBarRecoVsTrue());
00437     
00438     //Add in backgrounds
00439     nuPrediction.Add(wrongSignNuBars);
00440     barPrediction.Add(wrongSignNuMus);
00441     nuPrediction.Add(nuNCBackground);
00442     nuPrediction.Add(potentialNuTaus);
00443     barPrediction.Add(barNCBackground);
00444     barPrediction.Add(potentialTauBars);
00445     
00446     pair<NuMatrixSpectrum,NuMatrixSpectrum>
00447       ccPredictions(nuPrediction,barPrediction);
00448     return ccPredictions;
00449   }
00450 
00451   else if(fHelper){
00452     NuMatrixSpectrum potentialTaus(nuPrediction);
00453     potentialTaus.Multiply(fHelper->XSecGraphTaus());
00454     potentialTaus.Multiply(fHelper->FDTauEfficiency());
00455     
00456     nuPrediction.Multiply(fHelper->XSecGraph());
00457     nuPrediction.Multiply(fHelper->FDEfficiency());
00458     
00459     NuMatrixSpectrum unoscTrueSpectrum(nuPrediction);
00460     unoscTrueSpectrum.TrueToReco(fHelper->FDRecoVsTrue());
00461     NuMatrixSpectrum ncBackground(unoscTrueSpectrum);
00462     ncBackground.Divide(fHelper->FDPurity());
00463     ncBackground.Subtract(unoscTrueSpectrum);
00464     
00465     OscillateNuMu(nuPrediction, pars.Dm2(), pars.Sn2());
00466     InverseOscillateNuTau(potentialTaus, pars.Dm2(), pars.Sn2());
00467     
00468     potentialTaus.TrueToReco(fHelper->FDTauRecoVsTrue());
00469     nuPrediction.TrueToReco(fHelper->FDRecoVsTrue());
00470     nuPrediction.Add(ncBackground);
00471    
00472     nuPrediction.Add(potentialTaus);
00473     
00474     NuMatrixSpectrum dummy;
00475     pair<NuMatrixSpectrum,NuMatrixSpectrum>
00476       ccPredictions(nuPrediction,dummy);
00477     return ccPredictions;
00478   }
00479   else{
00480     //Woah, how did you get here without a helper?
00481     MSG("NuMMRunCC2010",Msg::kError)<<"Helper not set!"<<endl;
00482     NuMatrixSpectrum dummy;
00483     pair<NuMatrixSpectrum,NuMatrixSpectrum>
00484       predictions(dummy,dummy);
00485     return predictions;
00486   }
00487 }
00488 
00489 
00490 //____________________________________________________________________72
00491 const Double_t NuMMRunCC2010
00492 ::ComparePredWithData(const NuMMParameters& pars)
00493 {
00494   const pair<NuMatrixSpectrum,NuMatrixSpectrum> predictions =
00495     this->MakeFDPred(pars);
00496   Double_t like = 0.0;
00497   
00498   like += this->StatsLikelihood(predictions.first.Spectrum(),
00499                                 fFDCCData->Spectrum());
00500   
00501   if (fHelperCPT) {
00502     like += this->StatsLikelihood(predictions.second.Spectrum(),
00503                                   fFDCCDataBar->Spectrum());
00504   }
00505   return like;
00506 }
00507 
00508 
00509 //____________________________________________________________________72
00510 void NuMMRunCC2010::OscillateNuMu(NuMatrixSpectrum& numu,
00511                                   double dm2, double sn2) const
00512 {
00513   switch(fDisappearanceModel){
00514   case 0:
00515     numu.OscillateLinearInterp(dm2, sn2);
00516     return;
00517   case 1:
00518     numu.DecayCC(dm2, sn2);
00519     return;
00520   case 2:
00521     numu.Decohere(dm2, sn2);
00522     return;
00523   default:
00524     assert(0 && "Badly configured disappearance model");
00525   }
00526 }
00527 
00528 
00529 //____________________________________________________________________72
00530 void NuMMRunCC2010::InverseOscillateNuTau(NuMatrixSpectrum& nutau,
00531                                           double dm2, double sn2) const
00532 {
00533   if(fDisappearanceModel == 0)
00534     nutau.InverseOscillateLinearInterp(dm2, sn2);
00535 }
Generated on Mon Feb 15 11:07:15 2010 for loon by  doxygen 1.3.9.1