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

NuMMRunTransition.cxx

Go to the documentation of this file.
00001 #include "TCanvas.h"
00002 #include "TFile.h"
00003 #include "TH1D.h"
00004 #include "TGraph.h"
00005 
00006 #include "MessageService/MsgService.h"
00007 #include "NtupleUtils/NuHistInterpolator.h"
00008 #include "NtupleUtils/NuMatrixSpectrum.h"
00009 #include "NtupleUtils/NuMMHelperCPT.h"
00010 #include "NtupleUtils/NuMMParameters.h"
00011 #include "NtupleUtils/NuMMRunTransition.h"
00012 
00013 ClassImp(NuMMRunTransition)
00014 
00015 CVSID("$Id: NuMMRunTransition.cxx,v 1.18 2009/10/15 16:40:55 ahimmel Exp $");
00016 
00017 using namespace Sample;
00018 
00019 //____________________________________________________________________72
00020 NuMMRunTransition::NuMMRunTransition(NuMMHelperCPT* helper,
00021                                      NuMatrixSpectrum* ndNuData,
00022                                      NuMatrixSpectrum* ndBarData,
00023                                      double pot)
00024   : NuMMRunNuBar(),
00025     fPredictNeutrinos(false),
00026     fPredictAntiNeutrinos(true)
00027 {
00028   TString name = "hError";
00029   name += (counter++);
00030   TH1D hBlank(name, "Something Went Wrong", 10, 0, 30);
00031   
00032   fndNuData = ndNuData;
00033   fndBarData = ndBarData;
00034   fHelper = helper;
00035   
00036   ffdNuData = new NuMatrixSpectrum(hBlank, pot);
00037   ffdBarData = new NuMatrixSpectrum(hBlank, pot);
00038   
00039   
00040   // Blank out all the Precache variables
00041   PreCalcnuPrediction   = PreCalcbarPrediction      = 0;
00042   PreCalcnuAppeared     = PreCalcbarAppeared        = 0;
00043   PreCalcpotentialNuTaus= PreCalcpotentialTauBars   = 0;
00044   PreCalcnuNCBackground = PreCalcbarNCBackground    = 0;
00045   PreCalcwrongSignNuMus = PreCalcwrongSignNuBars    = 0;  
00046   
00047   PreCalcnuPredictionReco   = PreCalcbarPredictionReco      = 0;
00048   PreCalcnuAppearedReco     = PreCalcbarAppearedReco        = 0;
00049   PreCalcpotentialNuTausReco= PreCalcpotentialTauBarsReco   = 0;
00050   PreCalcwrongSignNuMusReco = PreCalcwrongSignNuBarsReco    = 0;  
00051   
00052 }
00053 
00054 //____________________________________________________________________72
00055 NuMMRunTransition::NuMMRunTransition(NuMMHelperCPT* helper,
00056                                      NuMatrixSpectrum* ndNuData,
00057                                      NuMatrixSpectrum* ndBarData,
00058                                      NuMatrixSpectrum* fdNuData,
00059                                      NuMatrixSpectrum* fdBarData)
00060   : NuMMRunNuBar(),
00061     fPredictNeutrinos(false),
00062     fPredictAntiNeutrinos(true)
00063 { 
00064   fndNuData = ndNuData;
00065   fndBarData = ndBarData;
00066   ffdNuData = fdNuData;
00067   ffdBarData = fdBarData;
00068   fHelper = helper;
00069    
00070   
00071   // Blank out all the Precache variables
00072   PreCalcnuPrediction   = PreCalcbarPrediction      = 0;
00073   PreCalcnuAppeared     = PreCalcbarAppeared        = 0;
00074   PreCalcpotentialNuTaus= PreCalcpotentialTauBars   = 0;
00075   PreCalcnuNCBackground = PreCalcbarNCBackground    = 0;
00076   PreCalcwrongSignNuMus = PreCalcwrongSignNuBars    = 0;  
00077 
00078   PreCalcnuPredictionReco   = PreCalcbarPredictionReco      = 0;
00079   PreCalcnuAppearedReco     = PreCalcbarAppearedReco        = 0;
00080   PreCalcpotentialNuTausReco= PreCalcpotentialTauBarsReco   = 0;
00081   PreCalcwrongSignNuMusReco = PreCalcwrongSignNuBarsReco    = 0;  
00082   
00083 }
00084 
00085 //____________________________________________________________________72
00086 
00087 NuMMRunTransition::~NuMMRunTransition()
00088 {
00089   ResetCache();
00090 }
00091 
00092 //___________________________________________________________________
00093 
00094 void NuMMRunTransition::ResetCache()
00095 {
00096   // delete the precalculated MM spectra we previously created
00097   if (PreCalcnuPrediction)    delete PreCalcnuPrediction;
00098   if (PreCalcbarPrediction)   delete PreCalcbarPrediction;
00099   
00100   if (PreCalcnuAppeared)      delete PreCalcnuAppeared;
00101   if (PreCalcbarAppeared)     delete PreCalcbarAppeared; 
00102   
00103   if (PreCalcpotentialNuTaus) delete PreCalcpotentialNuTaus; 
00104   if (PreCalcpotentialTauBars)delete PreCalcpotentialTauBars; 
00105   
00106   if (PreCalcnuNCBackground)  delete PreCalcnuNCBackground; 
00107   if (PreCalcbarNCBackground) delete PreCalcbarNCBackground; 
00108   
00109   if (PreCalcwrongSignNuMus)  delete PreCalcwrongSignNuMus; 
00110   if (PreCalcwrongSignNuBars) delete PreCalcwrongSignNuBars; 
00111   
00112   if (PreCalcnuPredictionReco)    delete PreCalcnuPredictionReco;
00113   if (PreCalcbarPredictionReco)   delete PreCalcbarPredictionReco;
00114   
00115   if (PreCalcnuAppearedReco)      delete PreCalcnuAppearedReco;
00116   if (PreCalcbarAppearedReco)     delete PreCalcbarAppearedReco;
00117   
00118   if (PreCalcpotentialNuTausReco) delete PreCalcpotentialNuTausReco; 
00119   if (PreCalcpotentialTauBarsReco)delete PreCalcpotentialTauBarsReco; 
00120   
00121   if (PreCalcwrongSignNuMusReco)  delete PreCalcwrongSignNuMusReco; 
00122   if (PreCalcwrongSignNuBarsReco) delete PreCalcwrongSignNuBarsReco; 
00123   
00124   PreCalcnuPrediction   = PreCalcbarPrediction      = 0;
00125   PreCalcnuAppeared     = PreCalcbarAppeared        = 0;
00126   PreCalcpotentialNuTaus= PreCalcpotentialTauBars   = 0;
00127   PreCalcnuNCBackground = PreCalcbarNCBackground    = 0;
00128   PreCalcwrongSignNuMus = PreCalcwrongSignNuBars    = 0;  
00129   
00130   PreCalcnuPredictionReco   = PreCalcbarPredictionReco      = 0;
00131   PreCalcnuAppearedReco     = PreCalcbarAppearedReco        = 0;
00132   PreCalcpotentialNuTausReco= PreCalcpotentialTauBarsReco   = 0;
00133   PreCalcwrongSignNuMusReco = PreCalcwrongSignNuBarsReco    = 0; 
00134   
00135   fCached = false;
00136 }
00137 
00138 //___________________________________________________________________
00139 
00140 void NuMMRunTransition::CacheExtrapolation(const NuMMParameters& pars)
00141 {
00142   // Only cache once.
00143   if (fCached) return;
00144   
00145   ResetCache();
00146   
00147   if (!fQuietMode)
00148     MSG("NuMMRunTransitions",Msg::kInfo) << "Pre-calculating static extrapolation variables..." << endl;
00149   
00150   // Caching variables, to reduce the workload of each step of the fitter
00151   // Initialise these the way they would be before
00152   PreCalcnuPrediction = new NuMatrixSpectrum(*fndNuData);
00153   PreCalcbarPrediction = new NuMatrixSpectrum(*fndBarData);
00154   
00155   // Be lazy and just alias the pointers to references so can keep
00156   // the code the same
00157   NuMatrixSpectrum &nuPrediction = *PreCalcnuPrediction;
00158   NuMatrixSpectrum &barPrediction = *PreCalcbarPrediction;
00159 
00160   
00161   // Now we have prepared the extrapolation variables, calculate them
00162   // Get the neutrinos to the FD
00163   nuPrediction.Multiply(fHelper->NDNuPurity());
00164   nuPrediction.RecoToTrue(fHelper->NDNuRecoVsTrue());
00165   nuPrediction.Divide(fHelper->NDNuEfficiency());
00166   nuPrediction.Divide(fHelper->NuXSecGraph());
00167   nuPrediction.Divide(fndNuData->PoT());
00168   nuPrediction.Divide(fNDFidMass);
00169   nuPrediction.ExtrapolateNDToFD(fHelper->NuBeamMatrix());
00170   nuPrediction.Multiply(fFDFidMass);
00171   nuPrediction.Multiply(ffdNuData->PoT());
00172   // Reset the spectrum to the new PoT
00173   nuPrediction.ResetPoT(ffdNuData->PoT());
00174   
00175   //Get the antineutrinos to the FD
00176   barPrediction.Multiply(fHelper->NDBarPurity());
00177   barPrediction.RecoToTrue(fHelper->NDBarRecoVsTrue());
00178   barPrediction.Divide(fHelper->NDBarEfficiency());
00179   barPrediction.Divide(fHelper->BarXSecGraph());
00180   barPrediction.Divide(fndBarData->PoT());
00181   barPrediction.Divide(fNDFidMass);
00182   barPrediction.ExtrapolateNDToFD(fHelper->BarBeamMatrix());
00183   barPrediction.Multiply(fFDFidMass);
00184   barPrediction.Multiply(ffdBarData->PoT());
00185   // Reset the spectrum to the new PoT
00186   barPrediction.ResetPoT(ffdBarData->PoT());
00187   
00188   PreCalcpotentialNuTaus = new NuMatrixSpectrum(nuPrediction);
00189   PreCalcpotentialTauBars = new NuMatrixSpectrum(barPrediction);  
00190   // Create and link to the Tau components, before we include the
00191   // cross section in the predictions (tau Xsec is different)
00192   NuMatrixSpectrum &potentialNuTaus = *PreCalcpotentialNuTaus;
00193   NuMatrixSpectrum &potentialTauBars = *PreCalcpotentialTauBars;
00194   
00195   if (fHelper->DoTaus()) {
00196     //Get the taus (ignoring any wrong-sign taubars)
00197     potentialNuTaus.Multiply(fHelper->XSecGraphNuTaus());
00198     potentialNuTaus.Multiply(fHelper->FDNuTauEfficiency());
00199     potentialNuTaus.InverseOscillate(pars.Dm2(),pars.Sn2());
00200   
00201     //Get the taubars (ignoring any wrong-sign taus)
00202     potentialTauBars.Multiply(fHelper->XSecGraphTauBars());
00203     potentialTauBars.Multiply(fHelper->FDTauBarEfficiency());
00204     potentialTauBars.InverseOscillate(pars.Dm2Bar(),pars.Sn2Bar());
00205   }
00206   else {
00207     PreCalcpotentialNuTaus->Multiply(0.0);
00208     PreCalcpotentialTauBars->Multiply(0.0);
00209   }
00210   
00211   // Create and link to the appeared components, before we include the
00212   // cross section in the predictions (Xsec is different)
00213   PreCalcnuAppeared = new NuMatrixSpectrum(barPrediction);
00214   PreCalcbarAppeared = new NuMatrixSpectrum(nuPrediction);
00215   NuMatrixSpectrum &appearedNuMus = *PreCalcnuAppeared;
00216   NuMatrixSpectrum &appearedNuBars = *PreCalcbarAppeared;
00217   appearedNuMus.SetName("Appeared_NuMus");
00218   appearedNuBars.SetName("Appeared_NuMuBars");
00219   
00220   // Apply the appeared cross-sectinos, efficiencies, true to reco
00221   appearedNuMus.Multiply(fHelper->NuXSecGraph());
00222   appearedNuMus.InverseOscillate(pars.Dm2Bar(),pars.Sn2Bar());
00223   appearedNuMus.Multiply(fHelper->FDNuEfficiency());
00224 
00225   appearedNuBars.Multiply(fHelper->BarXSecGraph());
00226   appearedNuBars.InverseOscillate(pars.Dm2(),pars.Sn2());
00227   appearedNuBars.Multiply(fHelper->FDBarEfficiency());
00228   
00229   
00230   // Now do the main spectrum, xsecs
00231   nuPrediction.Multiply(fHelper->NuXSecGraph());
00232   barPrediction.Multiply(fHelper->BarXSecGraph());
00233   
00234   
00235   // Now build the NC Spectrum, after the cross section calculation!
00236   PreCalcnuNCBackground = new NuMatrixSpectrum(*PreCalcnuPrediction);;
00237   PreCalcbarNCBackground = new NuMatrixSpectrum(*PreCalcbarPrediction);
00238   NuMatrixSpectrum &nuNCBackground = *PreCalcnuNCBackground;
00239   NuMatrixSpectrum &barNCBackground = *PreCalcbarNCBackground;  
00240 
00241   
00242   //Get the neutrino NC background
00243   nuNCBackground.Multiply(fHelper->FDNuEfficiency());
00244   nuNCBackground.TrueToReco(fHelper->FDNuRecoVsTrue());
00245   nuNCBackground.Divide(fHelper->FDNuPurity());
00246   nuNCBackground.Multiply(fHelper->FDNuNCContamination());
00247   
00248   //Get the antineutrino NC background
00249   barNCBackground.Multiply(fHelper->FDBarEfficiency());
00250   barNCBackground.TrueToReco(fHelper->FDBarRecoVsTrue());
00251   barNCBackground.Divide(fHelper->FDBarPurity());
00252   barNCBackground.Multiply(fHelper->FDBarNCContamination());
00253   
00254   
00255   
00256   // Now build the wrong sign spectrum
00257   PreCalcwrongSignNuMus = new NuMatrixSpectrum(*PreCalcnuPrediction);
00258   PreCalcwrongSignNuBars = new NuMatrixSpectrum(*PreCalcbarPrediction);
00259   NuMatrixSpectrum &wrongSignNuMus = *PreCalcwrongSignNuMus;
00260   NuMatrixSpectrum &wrongSignNuBars = *PreCalcwrongSignNuBars;
00261 
00262   // Do the Neutrino background for the antineutrino spectrum
00263   wrongSignNuMus.Multiply(fHelper->FDWrongSignNuEfficiency());
00264   wrongSignNuMus.Oscillate(pars.Dm2(),pars.Sn2());
00265   
00266   // Do the Antineutrino background for the neutrino spectrum
00267   wrongSignNuBars.Multiply(fHelper->FDWrongSignBarEfficiency());
00268   wrongSignNuBars.Oscillate(pars.Dm2Bar(),pars.Sn2Bar());
00269   
00270   
00271   // Now apply oscillations, efficeicnies to the main spectrum
00272   nuPrediction.Oscillate(pars.Dm2(),pars.Sn2());
00273   nuPrediction.Multiply(fHelper->FDNuEfficiency());
00274   
00275   barPrediction.Oscillate(pars.Dm2Bar(),pars.Sn2Bar());
00276   barPrediction.Multiply(fHelper->FDBarEfficiency());
00277   
00278 
00279   // If we have precalculated the wrong sign oscillations,
00280   // then we need to recalculate them next time we are given
00281   // the oscilaltion parameters. Reset the flags used to
00282   // control this
00283   //fDoneWsNeutrinos = fDoneWsAntineutrinos = false;
00284 
00285   nuPrediction.SetName("NuMus");
00286   barPrediction.SetName("NuMuBars");
00287 
00288   potentialNuTaus.SetName("NuTaus");
00289   potentialTauBars.SetName("NuTauBars");
00290   
00291   wrongSignNuMus.SetName("WrongSign_NuMus");
00292   wrongSignNuBars.SetName("WrongSign_NuMuBars");
00293   
00294   nuNCBackground.SetName("NC_NuMus");
00295   barNCBackground.SetName("NC_NuMuBars");
00296   
00297   // Now perform all the Reco To True's for the Reco prediction
00298   PreCalcnuPredictionReco = new NuMatrixSpectrum(*PreCalcnuPrediction);
00299   PreCalcbarPredictionReco = new NuMatrixSpectrum(*PreCalcbarPrediction);
00300   PreCalcnuPredictionReco->TrueToReco(fHelper->FDNuRecoVsTrue());
00301   PreCalcbarPredictionReco->TrueToReco(fHelper->FDBarRecoVsTrue());
00302   
00303   if (fHelper->DoTaus()) {
00304     PreCalcpotentialNuTausReco = new NuMatrixSpectrum(*PreCalcpotentialNuTaus);
00305     PreCalcpotentialTauBarsReco = new NuMatrixSpectrum(*PreCalcpotentialTauBars);
00306     PreCalcpotentialNuTausReco->TrueToReco(fHelper->FDNuTauRecoVsTrue());
00307     PreCalcpotentialTauBarsReco->TrueToReco(fHelper->FDTauBarRecoVsTrue());    
00308   }
00309   
00310   PreCalcnuAppearedReco = new NuMatrixSpectrum(*PreCalcnuAppeared);
00311   PreCalcbarAppearedReco = new NuMatrixSpectrum(*PreCalcbarAppeared);
00312   PreCalcnuAppearedReco->TrueToReco(fHelper->FDNuRecoVsTrue());
00313   PreCalcbarAppearedReco->TrueToReco(fHelper->FDBarRecoVsTrue());
00314   
00315   PreCalcwrongSignNuMusReco = new NuMatrixSpectrum(*PreCalcwrongSignNuMus);
00316   PreCalcwrongSignNuBarsReco = new NuMatrixSpectrum(*PreCalcwrongSignNuBars);
00317   PreCalcwrongSignNuMusReco->TrueToReco(fHelper->FDWrongSignNuRecoVsTrue());    
00318   PreCalcwrongSignNuBarsReco->TrueToReco(fHelper->FDWrongSignBarRecoVsTrue());
00319   
00320   fCached = true;
00321     
00322   // Disable this output for now - this process should be fast anyway!
00323   if (!fQuietMode)
00324     cout << "Extrapolation variables cached..." << endl;
00325 }
00326 
00327 
00328 //____________________________________________________________________72
00329 const std::vector<TH1D>
00330 NuMMRunTransition::WriteFDPredHistos(const NuMMParameters& pars) const 
00331 {
00332   //Set up a vector to push the histograms into.
00333   vector<TH1D> vHistos;
00334 
00335   //Put the ND data in the vector:
00336   TH1D hNDNuData(*(fndNuData->Spectrum()));
00337   hNDNuData.SetName("ndDataNQ");
00338   hNDNuData.SetTitle("ndDataNQ");
00339   vHistos.push_back(*(new TH1D(hNDNuData)));
00340   TH1D hNDBarData(*(fndBarData->Spectrum()));
00341   hNDBarData.SetName("ndDataPQ");
00342   hNDBarData.SetTitle("ndDataPQ");
00343   vHistos.push_back(*(new TH1D(hNDBarData)));
00344 
00345 //cout << "sn2: " << pars.Sn2() << "; dm2: " << pars.Dm2() 
00346 //     << "; sn2bar: " << pars.Sn2Bar() << "; dm2bar: " << pars.Dm2Bar()
00347 //     << "; Prob: " << pars.TransitionProb()
00348 //     << endl;
00349   //Get the neutrinos to the FD
00350   NuMatrixSpectrum nuPrediction(*fndNuData);
00351   nuPrediction.Multiply(fHelper->NDNuPurity());
00352   nuPrediction.RecoToTrue(fHelper->NDNuRecoVsTrue());
00353   nuPrediction.Divide(fHelper->NDNuEfficiency());
00354   nuPrediction.Divide(fHelper->NuXSecGraph());
00355   nuPrediction.Divide(fndNuData->PoT());
00356   nuPrediction.Divide(fNDFidMass);
00357   nuPrediction.ExtrapolateNDToFD(fHelper->NuBeamMatrix());
00358   nuPrediction.Multiply(fFDFidMass);
00359   nuPrediction.Multiply(ffdNuData->PoT());
00360 
00361   //Get the antineutrinos to the FD
00362   NuMatrixSpectrum barPrediction(*fndBarData);
00363   barPrediction.Multiply(fHelper->NDBarPurity());
00364   barPrediction.RecoToTrue(fHelper->NDBarRecoVsTrue());
00365   barPrediction.Divide(fHelper->NDBarEfficiency());
00366   barPrediction.Divide(fHelper->BarXSecGraph());
00367   barPrediction.Divide(fndBarData->PoT());
00368   barPrediction.Divide(fNDFidMass);
00369   barPrediction.ExtrapolateNDToFD(fHelper->BarBeamMatrix());
00370   barPrediction.Multiply(fFDFidMass);
00371   barPrediction.Multiply(ffdBarData->PoT());
00372 
00373   NuMatrixSpectrum potentialNuTaus(nuPrediction);
00374   NuMatrixSpectrum potentialTauBars(barPrediction);
00375   if (fHelper->DoTaus())  {
00376     //Get the taus (ignoring any wrong-sign taubars)
00377     potentialNuTaus.Multiply(fHelper->XSecGraphNuTaus());
00378     potentialNuTaus.Multiply(fHelper->FDNuTauEfficiency());
00379     potentialNuTaus.InverseOscillate(pars.Dm2(),pars.Sn2());
00380     potentialNuTaus.TrueToReco(fHelper->FDNuTauRecoVsTrue());
00381     //Get the taubars (ignoring any wrong-sign taus)
00382     potentialTauBars.Multiply(fHelper->XSecGraphTauBars());
00383     potentialTauBars.Multiply(fHelper->FDTauBarEfficiency());
00384     potentialTauBars.InverseOscillate(pars.Dm2Bar(),pars.Sn2Bar());
00385     potentialTauBars.TrueToReco(fHelper->FDTauBarRecoVsTrue());
00386   }
00387   else {
00388     potentialNuTaus.Multiply(0.0);
00389     potentialTauBars.Multiply(0.0);
00390   }
00391   
00392   //Get the transitioned neutrinos (they're now antineutrinos)
00393   NuMatrixSpectrum appearedNuBars(nuPrediction);
00394   TH1D hRawAppearPQ(*(appearedNuBars.Spectrum()));
00395   hRawAppearPQ.SetName("hRawAppearPQ");
00396   vHistos.push_back(hRawAppearPQ); 
00397   
00398   appearedNuBars.Multiply(fHelper->BarXSecGraph());
00399   TH1D hTrueEffAppearPQ(*(appearedNuBars.Spectrum()));
00400   hTrueEffAppearPQ.SetName("hTrueEffAppearPQ");
00401   vHistos.push_back(hTrueEffAppearPQ); 
00402   
00403   appearedNuBars.Multiply(fHelper->FDBarEfficiency());
00404   TH1D hTrueAppearPQ(*(appearedNuBars.Spectrum()));
00405   hTrueAppearPQ.SetName("hTrueAppearPQ");
00406   vHistos.push_back(hTrueAppearPQ); 
00407   
00408   appearedNuBars.InverseOscillate(pars.Dm2(),pars.Sn2());
00409   appearedNuBars.Multiply(pars.TransitionProb());
00410   TH1D hTrueTransitionPQ(*(appearedNuBars.Spectrum()));
00411   hTrueTransitionPQ.SetName("hTrueTransitionPQ");
00412   vHistos.push_back(hTrueTransitionPQ); 
00413     
00414   appearedNuBars.TrueToReco(fHelper->FDBarRecoVsTrue());
00415   TH1D hRecoAppearPQ(*(appearedNuBars.Spectrum()));
00416   hRecoAppearPQ.SetName("hRecoAppearPQ");
00417   vHistos.push_back(hRecoAppearPQ);
00418   
00419   //Get the transitioned antineutrinos (they're now neutrinos)
00420   NuMatrixSpectrum appearedNuMus(barPrediction);
00421   TH1D hRawAppearNQ(*(appearedNuMus.Spectrum()));
00422   hRawAppearNQ.SetName("hRawAppearNQ");
00423   vHistos.push_back(hRawAppearNQ); 
00424   
00425   appearedNuMus.Multiply(fHelper->NuXSecGraph());
00426   TH1D hTrueEffAppearNQ(*(appearedNuMus.Spectrum()));
00427   hTrueEffAppearNQ.SetName("hTrueEffAppearNQ");
00428   vHistos.push_back(hTrueEffAppearNQ); 
00429   
00430   appearedNuMus.Multiply(fHelper->FDNuEfficiency());
00431   TH1D hTrueAppearNQ(*(appearedNuMus.Spectrum()));
00432   hTrueAppearNQ.SetName("hTrueAppearNQ");
00433   vHistos.push_back(hTrueAppearNQ); 
00434   
00435   appearedNuMus.InverseOscillate(pars.Dm2(),pars.Sn2());
00436   appearedNuMus.Multiply(pars.TransitionProb());
00437   TH1D hTrueTransitionNQ(*(appearedNuMus.Spectrum()));
00438   hTrueTransitionNQ.SetName("hTrueTransitionNQ");
00439   vHistos.push_back(hTrueTransitionNQ); 
00440   
00441   appearedNuMus.TrueToReco(fHelper->FDNuRecoVsTrue());
00442   TH1D hRecoAppearNQ(*(appearedNuMus.Spectrum()));
00443   hRecoAppearNQ.SetName("hRecoAppearNQ");
00444   vHistos.push_back(hRecoAppearNQ);
00445   
00446 
00447   //Cross sections
00448   nuPrediction.Multiply(fHelper->NuXSecGraph());
00449   barPrediction.Multiply(fHelper->BarXSecGraph());
00450   
00451   // True NuMus (flux x sigma)
00452   TH1D hTrueNuMus(*(nuPrediction.Spectrum()));
00453   hTrueNuMus.SetName("hTrueNuMus");
00454   vHistos.push_back(hTrueNuMus); 
00455   
00456   // True NuBars (flux x sigma)
00457   TH1D hTrueNuBars(*(barPrediction.Spectrum()));
00458   hTrueNuBars.SetName("hTrueNuBars");
00459   vHistos.push_back(hTrueNuBars); 
00460 
00461   //Get the neutrino NC background
00462   NuMatrixSpectrum nuNCBackground(nuPrediction);
00463   nuNCBackground.Multiply(fHelper->FDNuEfficiency());
00464   nuNCBackground.TrueToReco(fHelper->FDNuRecoVsTrue());
00465   nuNCBackground.Divide(fHelper->FDNuPurity());
00466   nuNCBackground.Multiply(fHelper->FDNuNCContamination());
00467   
00468   //Get the antineutrino NC background
00469   NuMatrixSpectrum barNCBackground(barPrediction);
00470   barNCBackground.Multiply(fHelper->FDBarEfficiency());
00471   barNCBackground.TrueToReco(fHelper->FDBarRecoVsTrue());
00472   barNCBackground.Divide(fHelper->FDBarPurity());
00473   barNCBackground.Multiply(fHelper->FDBarNCContamination());
00474   
00475   //Get the wrong-sign neutrino background for the antineutrino
00476   //prediction
00477   NuMatrixSpectrum wrongSignNuMus(nuPrediction);
00478   wrongSignNuMus.Multiply(fHelper->FDWrongSignNuEfficiency());
00479   wrongSignNuMus.Oscillate(pars.Dm2(),pars.Sn2());
00480   
00481   //Get the wrong-sign antineutrino background for the neutrino
00482   //prediction
00483   NuMatrixSpectrum wrongSignNuBars(barPrediction);
00484   wrongSignNuBars.Multiply(fHelper->FDWrongSignBarEfficiency());
00485   wrongSignNuBars.Oscillate(pars.Dm2Bar(),pars.Sn2Bar());  
00486 
00487 
00488   TH1D hFDWrongSignPQ(*(wrongSignNuBars.Spectrum()));
00489   hFDWrongSignPQ.SetName("fdWrongSignPQ");
00490   vHistos.push_back(hFDWrongSignPQ);
00491 
00492   TH1D hFDWrongSignNQ(*(wrongSignNuMus.Spectrum()));
00493   hFDWrongSignNQ.SetName("fdWrongSignNQ");
00494   vHistos.push_back(hFDWrongSignNQ);
00495 
00496   // Do the truetoreco seperately so that we can write out the histos first
00497   wrongSignNuBars.TrueToReco(fHelper->FDWrongSignBarRecoVsTrue());
00498   wrongSignNuMus.TrueToReco(fHelper->FDWrongSignNuRecoVsTrue()); 
00499 
00500   //Oscillations
00501   nuPrediction.Oscillate(pars.Dm2(),pars.Sn2());
00502   barPrediction.Oscillate(pars.Dm2Bar(),pars.Sn2Bar());
00503   //Efficiencies
00504   nuPrediction.Multiply(fHelper->FDNuEfficiency());
00505   barPrediction.Multiply(fHelper->FDBarEfficiency());
00506 
00507   //True to reco
00508   nuPrediction.TrueToReco(fHelper->FDNuRecoVsTrue());
00509   barPrediction.TrueToReco(fHelper->FDBarRecoVsTrue());
00510   
00511   // Dump out the prediction spectrum before adding in the backgrounds
00512   TH1D hFDNoBackNQ(*(nuPrediction.Spectrum()));
00513   hFDNoBackNQ.SetName("fdBasePredictionNQ");
00514   vHistos.push_back(hFDNoBackNQ);
00515   
00516   TH1D hFDNoBackPQ(*(barPrediction.Spectrum()));
00517   hFDNoBackPQ.SetName("fdBasePredictionPQ");
00518   vHistos.push_back(hFDNoBackPQ);
00519   
00520   //Add in backgrounds
00521   nuPrediction.Add(wrongSignNuBars);
00522   barPrediction.Add(wrongSignNuMus);
00523   nuPrediction.Add(nuNCBackground);
00524   if (fHelper->DoTaus()) nuPrediction.Add(potentialNuTaus);
00525   barPrediction.Add(barNCBackground);
00526   if (fHelper->DoTaus()) barPrediction.Add(potentialTauBars);
00527 
00528   
00529   TH1D hFDPredictionNoTransNQ(*(nuPrediction.Spectrum()));
00530   hFDPredictionNoTransNQ.SetName("fdPredictionNoTransNQ");
00531   vHistos.push_back(hFDPredictionNoTransNQ);
00532   
00533   TH1D hFDPredictionNoTransPQ(*(barPrediction.Spectrum()));
00534   hFDPredictionNoTransPQ.SetName("fdPredictionNoTransPQ");
00535   vHistos.push_back(hFDPredictionNoTransPQ);
00536   
00537   
00538   
00539   //Add in transitions
00540   nuPrediction.Add(appearedNuMus);
00541   barPrediction.Add(appearedNuBars);
00542   
00543   //Put the componants of the predictions in the vector
00544 
00545   TH1D hFDPredictionNQ(*(nuPrediction.Spectrum()));
00546   hFDPredictionNQ.SetName("fdPredictionNQ");
00547   vHistos.push_back(hFDPredictionNQ);
00548 
00549   TH1D hFDPredictionPQ(*(barPrediction.Spectrum()));
00550   hFDPredictionPQ.SetName("fdPredictionPQ");
00551   vHistos.push_back(hFDPredictionPQ);
00552 
00553   TH1D hFDNCNQ(*(nuNCBackground.Spectrum()));
00554   hFDNCNQ.SetName("fdNCNQ");
00555   vHistos.push_back(hFDNCNQ);
00556 
00557   TH1D hFDNCPQ(*(barNCBackground.Spectrum()));
00558   hFDNCPQ.SetName("fdNCPQ");
00559   vHistos.push_back(hFDNCPQ);
00560   
00561   if (fHelper->DoTaus()) {
00562     TH1D hFDTausNQ(*(potentialNuTaus.Spectrum()));
00563     hFDTausNQ.SetName("fdTausNQ");
00564     vHistos.push_back(hFDTausNQ);
00565     
00566     TH1D hFDTausPQ(*(potentialTauBars.Spectrum()));
00567     hFDTausPQ.SetName("fdTausPQ");
00568     vHistos.push_back(hFDTausPQ);
00569   }
00570   
00571   TH1D hFDDataNQ(*(ffdNuData->Spectrum()));
00572   hFDDataNQ.SetName("fdDataNQ");
00573   vHistos.push_back(hFDDataNQ);
00574 
00575   TH1D hFDDataPQ(*(ffdBarData->Spectrum()));
00576   hFDDataPQ.SetName("fdDataPQ");
00577   vHistos.push_back(hFDDataPQ);
00578 
00579   return vHistos;
00580 }
00581 //____________________________________________________________________72
00582 const pair<NuMatrixSpectrum,NuMatrixSpectrum>
00583 NuMMRunTransition::MakeFDPred(const NuMMParameters& pars)
00584 {
00585   if (!fQuietMode) {
00586     cout << "sn2: " << pars.Sn2() << "; dm2: " << pars.Dm2() 
00587      << "; sn2bar: " << pars.Sn2Bar() << "; dm2bar: " << pars.Dm2Bar()
00588      << "; alpha: " << pars.TransitionProb()
00589      << endl;
00590   }
00591   
00592   CacheExtrapolation(pars);
00593   
00595   // For all neutrinos
00596   
00597   // Make copies of the prepared extrapolation variables, so that
00598   // we can alter them locally with the oscillation parameters
00599   NuMatrixSpectrum nuPrediction(*PreCalcnuPredictionReco);
00600   NuMatrixSpectrum barPrediction(*PreCalcbarPredictionReco);
00601   
00602   NuMatrixSpectrum potentialNuTaus(*PreCalcpotentialNuTausReco);
00603   NuMatrixSpectrum potentialTauBars(*PreCalcpotentialTauBarsReco);
00604   
00605   NuMatrixSpectrum appearedNuMus(*PreCalcnuAppearedReco);
00606   NuMatrixSpectrum appearedNuBars(*PreCalcbarAppearedReco);
00607   
00608   NuMatrixSpectrum wrongSignNuMus(*PreCalcwrongSignNuMusReco);
00609   NuMatrixSpectrum wrongSignNuBars(*PreCalcwrongSignNuBarsReco);
00610   
00611   NuMatrixSpectrum nuNCBackground(*PreCalcnuNCBackground);
00612   NuMatrixSpectrum barNCBackground(*PreCalcbarNCBackground);  
00613   
00614   
00616   // Processes for neutrinos
00617   
00618   if (fPredictNeutrinos) {  
00619     // Apply alpha
00620     if (fHelper->DoTaus()) potentialNuTaus.Multiply(1.-pars.TransitionProb());
00621     appearedNuMus.Multiply(pars.TransitionProb());
00622     
00623     //Add everything together
00624     nuPrediction.Add(appearedNuMus);
00625     if (fHelper->DoTaus()) nuPrediction.Add(potentialNuTaus);    
00626     nuPrediction.Add(wrongSignNuBars);
00627     nuPrediction.Add(nuNCBackground);
00628   } else {
00629     // Scale the histogram to zero, so as not to confuse people
00630     // if they look at the spectrum with neutrinos turned off
00631     nuPrediction.Spectrum()->Scale(0);
00632   }
00633   
00635   // Processes for Antineutrinos
00636   
00637   if (fPredictAntiNeutrinos) {
00638     if (fHelper->DoTaus()) potentialTauBars.Multiply(1.-pars.TransitionProb());
00639     appearedNuBars.Multiply(pars.TransitionProb());
00640     
00641     //Add everything together
00642     barPrediction.Add(appearedNuBars);
00643     if (fHelper->DoTaus()) barPrediction.Add(potentialTauBars);
00644     barPrediction.Add(wrongSignNuMus);
00645     barPrediction.Add(barNCBackground);
00646   } else {
00647     // Scale the histogram to zero, so as not to confuse people
00648     barPrediction.Spectrum()->Scale(0);
00649   }
00650   
00651   pair<NuMatrixSpectrum,NuMatrixSpectrum> predictions(nuPrediction,barPrediction);
00652   return predictions;
00653 }
00654 
00655 //____________________________________________________________________72
00656 NuMatrixSpectrum* NuMMRunTransition::MakeFDPredNuMu(const NuMMParameters& pars)
00657 {
00658   if (!fQuietMode) {
00659     cout << "sn2: " << pars.Sn2() << "; dm2: " << pars.Dm2() 
00660     << "; sn2bar: " << pars.Sn2Bar() << "; dm2bar: " << pars.Dm2Bar()
00661     << "; alpha: " << pars.TransitionProb()
00662     << endl;
00663   }
00664   
00665   CacheExtrapolation(pars);
00666   
00668   // For all neutrinos
00669   
00670   // Make copies of the prepared extrapolation variables, so that
00671   // we can alter them locally with the oscillation parameters
00672   NuMatrixSpectrum *nuPrediction = new NuMatrixSpectrum(*PreCalcnuPredictionReco);
00673   
00674   NuMatrixSpectrum potentialNuTaus(*PreCalcpotentialNuTausReco);  
00675   NuMatrixSpectrum appearedNuMus(*PreCalcnuAppearedReco);  
00676   NuMatrixSpectrum wrongSignNuBars(*PreCalcwrongSignNuBarsReco);  
00677   NuMatrixSpectrum nuNCBackground(*PreCalcnuNCBackground);
00678   
00679   
00681   // Processes for neutrinos
00682   
00683   if (fPredictNeutrinos) {  
00684     // Apply alpha
00685     if (fHelper->DoTaus()) potentialNuTaus.Multiply(1.-pars.TransitionProb());
00686     appearedNuMus.Multiply(pars.TransitionProb());
00687     
00688     //Add everything together
00689     nuPrediction->Add(appearedNuMus);
00690     if (fHelper->DoTaus()) nuPrediction->Add(potentialNuTaus);    
00691     nuPrediction->Add(wrongSignNuBars);
00692     nuPrediction->Add(nuNCBackground);
00693   } else {
00694     // Scale the histogram to zero, so as not to confuse people
00695     // if they look at the spectrum with neutrinos turned off
00696     nuPrediction->Spectrum()->Scale(0);
00697   }
00698   
00699   return nuPrediction;
00700 }
00701 
00702 
00703 //____________________________________________________________________72
00704 NuMatrixSpectrum* NuMMRunTransition::MakeFDPredNuBar(const NuMMParameters& pars)
00705 {
00706   if (!fQuietMode) {
00707     cout << "sn2: " << pars.Sn2() << "; dm2: " << pars.Dm2() 
00708     << "; sn2bar: " << pars.Sn2Bar() << "; dm2bar: " << pars.Dm2Bar()
00709     << "; alpha: " << pars.TransitionProb()
00710     << endl;
00711   }
00712   
00713   CacheExtrapolation(pars);
00714   
00716   // For all neutrinos
00717   
00718   // Make copies of the prepared extrapolation variables, so that
00719   // we can alter them locally with the oscillation parameters
00720   NuMatrixSpectrum *barPrediction = new NuMatrixSpectrum(*PreCalcbarPredictionReco);
00721   
00722   NuMatrixSpectrum potentialTauBars(*PreCalcpotentialTauBarsReco);
00723   NuMatrixSpectrum appearedNuBars(*PreCalcbarAppearedReco);  
00724   NuMatrixSpectrum wrongSignNuMus(*PreCalcwrongSignNuMusReco);  
00725   NuMatrixSpectrum barNCBackground(*PreCalcbarNCBackground);  
00726   
00728   // Processes for Antineutrinos
00729   
00730   if (fPredictAntiNeutrinos) {
00731     if (fHelper->DoTaus()) potentialTauBars.Multiply(1.-pars.TransitionProb());
00732     appearedNuBars.Multiply(pars.TransitionProb());
00733     
00734     //Add everything together
00735     barPrediction->Add(appearedNuBars);
00736     if (fHelper->DoTaus()) barPrediction->Add(potentialTauBars);
00737     barPrediction->Add(wrongSignNuMus);
00738     barPrediction->Add(barNCBackground);
00739   } else {
00740     // Scale the histogram to zero, so as not to confuse people
00741     barPrediction->Spectrum()->Scale(0);
00742   }
00743 
00744   return barPrediction;
00745 }
00746 
00747 
00748 //____________________________________________________________________72
00749 
00750 const NuMatrixSpectrum NuMMRunTransition::TrueComponents(const NuMMParameters& pars, Sample_t s)
00751 {
00752   if (!fQuietMode) {
00753     cout << "sn2: " << pars.Sn2() << "; dm2: " << pars.Dm2() 
00754     << "; sn2bar: " << pars.Sn2Bar() << "; dm2bar: " << pars.Dm2Bar()
00755     << "; alpha: " << pars.TransitionProb()
00756     << endl;
00757   }
00758   
00759   // Request NQ spectrum, but no NuBar Predictions OR
00760   // Request PQ spectrum, but no NuMu Predictions
00761   if ( (s < 0 && !fPredictNeutrinos) || 
00762        (s > 0 && !fPredictAntiNeutrinos) ) {
00763     NuMatrixSpectrum blank(*PreCalcnuPrediction);
00764     blank.SetName("Blank");
00765     blank.Spectrum()->Scale(0);
00766     return blank;
00767   }
00768   
00769   CacheExtrapolation(pars);
00770   
00771 
00772   if (s == kSignalPQ) {
00773     if (!fQuietMode) cout << "Getting kSignalPQ" << endl;
00774     NuMatrixSpectrum barPrediction(*PreCalcbarPrediction);
00775     if (!fQuietMode) cout << "Returning kSignalPQ" << endl;
00776     return barPrediction;
00777   }
00778   else if (s == kSignalNQ) {
00779     if (!fQuietMode) cout << "Getting kSignalNQ" << endl;
00780     NuMatrixSpectrum nuPrediction(*PreCalcnuPrediction);  
00781     if (!fQuietMode) cout << "Returning kSignalNQ" << endl;
00782     return nuPrediction;
00783   }
00784   else if (s == kWrongSignPQ) {
00785     if (!fQuietMode) cout << "Getting kWrongSignPQ" << endl;
00786     NuMatrixSpectrum wrongSignNuMus(*PreCalcwrongSignNuMus);
00787     if (!fQuietMode) cout << "Returning kWrongSignPQ" << endl;
00788     return wrongSignNuMus;
00789   }
00790   else if (s == kWrongSignNQ) {
00791     if (!fQuietMode) cout << "Getting kWrongSignNQ" << endl;
00792     NuMatrixSpectrum wrongSignNuBars(*PreCalcwrongSignNuBars);
00793     if (!fQuietMode) cout << "Returning kWrongSignNQ" << endl;
00794     return wrongSignNuBars;
00795   }  
00796   else if (s == kNCPQ) {
00797     if (!fQuietMode) cout << "Getting kNCPQ" << endl;
00798     NuMatrixSpectrum barNCBackground(*PreCalcbarNCBackground);
00799     if (!fQuietMode) cout << "Returning kNCPQ" << endl;
00800     return barNCBackground;
00801   }
00802   else if (s == kNCNQ) {
00803     if (!fQuietMode) cout << "Getting kNCNQ" << endl;
00804     NuMatrixSpectrum nuNCBackground(*PreCalcnuNCBackground);
00805     if (!fQuietMode) cout << "Returning kNCNQ" << endl;
00806     return nuNCBackground;
00807   }
00808   else if (s == kTauPQ) {
00809     if (!fQuietMode) cout << "Getting kTauPQ" << endl;
00810     NuMatrixSpectrum potentialTauBars(*PreCalcpotentialTauBars);
00811     potentialTauBars.Multiply(1.-pars.TransitionProb());
00812     if (!fQuietMode) cout << "Returning kTauPQ" << endl;
00813     return potentialTauBars;
00814   }
00815   else if (s == kTauNQ) {
00816     if (!fQuietMode) cout << "Getting kTauNQ" << endl;
00817     NuMatrixSpectrum potentialNuTaus(*PreCalcpotentialNuTaus);
00818     potentialNuTaus.Multiply(1.-pars.TransitionProb());
00819     if (!fQuietMode) cout << "Returning kTauNQ" << endl;
00820     return potentialNuTaus;
00821   }
00822   else if (s == kAppearedPQ) {
00823     if (!fQuietMode) cout << "Getting kAppearedPQ" << endl;
00824     NuMatrixSpectrum appearedNuBars(*PreCalcbarAppeared);
00825     appearedNuBars.Multiply(pars.TransitionProb());
00826     if (!fQuietMode) cout << "Returning kAppearedPQ" << endl;
00827     return appearedNuBars;
00828   }
00829   else if (s == kAppearedNQ) {
00830     if (!fQuietMode) cout << "Getting kAppearedNQ" << endl;
00831     NuMatrixSpectrum appearedNuMus(*PreCalcnuAppeared);
00832     appearedNuMus.Multiply(pars.TransitionProb());
00833     if (!fQuietMode) cout << "Returning kAppearedNQ" << endl;
00834     return appearedNuMus;
00835   }
00836   else {
00837     MSG("NuMMRunTransitions",Msg::kWarning) << "Requested unknown prediction #" << s << ", returning blank." << endl;
00838     NuMatrixSpectrum blank(*PreCalcnuPrediction);
00839     blank.SetName("Blank");
00840     blank.Spectrum()->Scale(0);
00841     return blank;
00842   }
00843 }
00844 
00845 //____________________________________________________________________
00846 
00847 const Double_t NuMMRunTransition::ComparePredWithData(const NuMMParameters& pars)
00848 {
00849   // Generate a prediction
00850   const pair<NuMatrixSpectrum,NuMatrixSpectrum> predictions =
00851     this->MakeFDPred(pars);
00852   Double_t like = 0;
00853   // Now, only count a likelihood for components we want - this is just
00854   // a constant offset, but should help speed things up a little.
00855   if (fPredictNeutrinos) {
00856     NuMatrixSpectrum* nuPrediction = MakeFDPredNuMu(pars);
00857     like += this->StatsLikelihood(nuPrediction->Spectrum(),
00858                                   ffdNuData->Spectrum());
00859     delete nuPrediction;
00860   }
00861   if (fPredictAntiNeutrinos) {
00862     NuMatrixSpectrum* barPrediction = MakeFDPredNuBar(pars);
00863     like += this->StatsLikelihood(barPrediction->Spectrum(),
00864                                   ffdBarData->Spectrum());
00865     delete barPrediction;
00866   }
00867   return like;
00868 }
Generated on Mon Feb 15 11:07:15 2010 for loon by  doxygen 1.3.9.1