loon: NuMatrixSpectrum.h Source File

00001 #ifndef NUMATRIXSPECTRUM_H
00002 #define NUMATRIXSPECTRUM_H
00003 
00004 #include <string>
00005 
00006 #include "TArrayD.h"
00007 #include "TH1D.h"
00008 
00009 #include "NtupleUtils/NuMatrix1D.h"
00010 #include "NtupleUtils/NuUtilities.h" // For NuBinningScheme
00011 
00012 class TGraph;
00013 class TH2D;
00014 
00015 class NuMatrixSpectrum : public NuMatrix1D
00016 {
00017  public:
00018   // *** Constructors ***
00019   NuMatrixSpectrum();
00020   explicit NuMatrixSpectrum(const TH1D& spectrum);
00021   NuMatrixSpectrum(const TH1D& spectrum,
00022                             const Double_t pot);
00024   NuMatrixSpectrum(const Double_t pot);
00026   NuMatrixSpectrum(const std::string& filename, const std::string& histname,
00027                    const std::string& POTHistname = "hTotalPot");
00028   NuMatrixSpectrum(TFile* specfile, const std::string& histname,
00029                    const std::string& POTHistname = "hTotalPot");
00030   virtual ~NuMatrixSpectrum();
00031   
00032   const char* GetName() const;
00033   void SetName(const char* _name);
00034 
00035   // NuMatrix1D has maths operations in terms of concrete objects, whereas we
00036   // have them in terms of pointers. Fun times converting...
00037   // The "using" declarations are neeeded to prevent hiding NuMatrix1D versions
00038   // *** Multiply ***
00039   using NuMatrix1D::Multiply;
00040   virtual void Multiply(const TH1* c){Multiply(*c);}
00041   virtual void Multiply(const NuMatrixSpectrum* c){Multiply(c->Spectrum());}
00042   virtual void Multiply(const TGraph* c){Multiply(*c);}
00043   // This one is actually unique to NuMatrixSpectrum
00044   virtual void Multiply(TF1* c);
00045 
00046   // *** Divide ***
00047   using NuMatrix1D::Divide;
00048   // We propogate Poisson errors here
00049   virtual void Divide(const TH1* correction, Option_t* option = "");
00050   virtual void Divide(const NuMatrixSpectrum* c, Option_t* option = "")
00051     {Divide(c->Spectrum(), option);}
00052   virtual void Divide(const NuMatrixSpectrum c, Option_t* option = "")
00053     {Divide(c.Spectrum(), option);}
00054   virtual void Divide(const TGraph* correction){Divide(*correction);}
00055 
00056   // *** Purity ***
00057   virtual void MultiplyPurity(const TH1D* corr, const Double_t corrShift);
00058   virtual void DividePurity(const TH1D* corr, const Double_t corrShift);
00059 
00060   // *** Reco <-> True
00061   using NuMatrix1D::RecoToTrue;
00062   virtual void RecoToTrue(const TH2D* correction){RecoToTrue(*correction);}
00063   using NuMatrix1D::TrueToReco;
00064   virtual void TrueToReco(const TH2D* correction){TrueToReco(*correction);}
00065 
00066   // *** Errors ***
00067   virtual void GenerateSystErrors();
00068   virtual void GenerateErrors(Double_t scale = 1);
00069   virtual void RemoveErrors();
00070 
00071   // ** Presentation ***
00072   virtual void BinWidthNormalize(double width = -1);
00073   virtual void NuBarCompressX();
00074   virtual void NCCompressX();
00075   virtual void CC2010CompressX();
00076   virtual void CC2014CompressX();
00077   virtual void CompressTopBins();
00078   static TF1* GetAxisFunc();
00079   virtual void DrawNuBarAxes(int lsize = -1);
00080   virtual void DrawNCAxes(int lsize = -1);
00081   virtual void DrawCC2010Axes(bool labels = true);
00082   virtual void DrawCC2014Axes(bool labels = true);
00083   virtual void CutLastBin();
00084 
00085   virtual void BlessedND(bool rebin = true, bool scale = true);
00086   virtual void BlessedFD(bool scale = true);
00087   virtual void BlessedForNC(void);
00088   //  virtual void BlessedNDNC(void);
00089   virtual void BlessedFDRatio(bool scale = true);
00090   virtual void SetFontAndSize(int textFont, int textSize);
00091   virtual void BlessedFD2010();
00092   virtual void BlessedFD2014();
00093   virtual void BlessedFD7e20(bool scale = true);
00094 
00095   virtual const TH2D MakeJessTemplate(const TH2D* recoVtrue) const;
00096  protected:
00097   virtual Double_t PoiErr(double y, bool up);
00098   virtual void CalcPoisson();
00099  public:
00100   virtual void DrawPoisson(Option_t* opt = "");
00101   static double compress;
00102   static double xlimit;
00103 
00104 
00105   // Create some pass-through methods to make life easier
00106   virtual void SetLineColor(Color_t color){fSpectrum->SetLineColor(color);}
00107   virtual void SetLineWidth(Width_t width){fSpectrum->SetLineWidth(width);}
00108   virtual void SetLineStyle(Style_t style){fSpectrum->SetLineStyle(style);}
00109   virtual void SetMarkerColor(Color_t col){fSpectrum->SetMarkerColor(col);}
00110   virtual void SetMarkerSize(Size_t size) {fSpectrum->SetMarkerSize(size);}
00111   virtual void SetMarkerStyle(Style_t sty){fSpectrum->SetMarkerStyle(sty);}
00112   virtual void SetFillColor(Color_t color){fSpectrum->SetFillColor(color);}
00113   virtual void SetFillStyle(Style_t style){fSpectrum->SetFillStyle(style);}
00114 
00115   virtual void SetTitle(const char* title){fSpectrum->SetTitle(title);}
00116   virtual void CenterTitles(Bool_t cent = true)
00117   {
00118     fSpectrum->GetXaxis()->CenterTitle(cent);
00119     fSpectrum->GetYaxis()->CenterTitle(cent);
00120   }
00121   virtual void SetRangeUser(double low, double high){fSpectrum->GetXaxis()->SetRangeUser(low, high);}
00122   virtual TAxis* GetXaxis(){return fSpectrum->GetXaxis();}
00123   virtual TAxis* GetYaxis(){return fSpectrum->GetYaxis();}
00124 
00125   virtual Double_t Integral(Option_t* op = ""){return fSpectrum->Integral(op);}
00126   virtual Double_t Integral(Int_t binx1, Int_t binx2, Option_t* option = "")
00127   {
00128     return fSpectrum->Integral(binx1, binx2, option);
00129   }
00130   virtual Double_t IntegralVals(Double_t low, Double_t high, Option_t* option = "");
00131 
00132   virtual Int_t GetNbinsX() const {return fSpectrum->GetNbinsX();}
00133   virtual Double_t GetEntries() const {return fSpectrum->GetEntries();}
00134 
00135   // *** Extrapolate ***
00136   using NuMatrix1D::ExtrapolateNDToFD;
00137   virtual void ExtrapolateNDToFD(const TH2D* corr){ExtrapolateNDToFD(*corr);}
00138 
00139   // *** Oscillations ***
00140   virtual void Oscillate(Double_t dm2, Double_t sn2);
00141   virtual void Oscillate_3flavour(Double_t dm232, Double_t t23, Double_t dm221, Double_t t12, Double_t t13, Double_t delta, int inu);
00142   virtual void OscillateSimpleAverage(Double_t dm2, Double_t sn2);
00143   virtual void OscillatePreAveraged(Double_t dm2, Double_t sn2);
00144   virtual void OscillateLinearInterp(Double_t dm2, Double_t sn2);
00145   virtual void OscillateHybrid(Double_t dm2, Double_t sn2);
00146   virtual void OscillateHybrid_3flavour(Double_t dm232, Double_t t23, Double_t dm221, Double_t t12, Double_t t13, Double_t delta, int inu);
00147   virtual void OscillateBasic_3flavour(Double_t dm232, Double_t t23, Double_t dm221, Double_t t12, Double_t t13, Double_t delta, int inu);
00148   virtual void OscillateNSI(Double_t dm2, Double_t sn2,
00149                             Double_t epsilon, Double_t sign);
00150 
00151   // for large extra dimension (LED)
00152   void OscillateLED(int initialFlavor, int finalFlavor, int cp_sign,
00153                     double th12, double th13, double th23,
00154                     double deltacp,
00155                     double sdm, double ldm,
00156                     double m0, double a);
00157   
00158  //for transitions SME
00159   virtual void OscillateMu2Mu(Double_t dm2, Double_t sn2, Double_t gmu, Double_t gtau, Double_t cmumu, Double_t ctautau); 
00160   virtual void OscillateMu2MuBar(Double_t dm2, Double_t sn2, Double_t gmu, Double_t gtau, Double_t cmumu, Double_t ctautau); 
00161   virtual void OscillateMu2Tau(Double_t dm2, Double_t sn2, Double_t gmu, Double_t gtau, Double_t cmumu, Double_t ctautau); 
00162   virtual void OscillateMu2TauBar(Double_t dm2, Double_t sn2, Double_t gmu, Double_t gtau, Double_t cmumu, Double_t ctautau); 
00163 
00164   void OscillateFourFlavour(const double dm2, const double theta23,
00165                             const double dm221, const double dm243,
00166                             const double theta12,
00167                             const double theta13, const double theta14,
00168                             const double theta24, const double theta34,
00169                             const double delta1, const double delta2,
00170                             const double delta3,
00171                             const double baseline);
00172   void OscillateFourFlavourAverage(const double dm2, const double theta23,
00173                                    const double dm221, const double dm243,
00174                                    const double theta12,
00175                                    const double theta13, const double theta14,
00176                                    const double theta24, const double theta34,
00177                                    const double delta1, const double delta2,
00178                                    const double delta3,
00179                                    const double baseline);
00180   void UnoscillateND(const double dm2, const double theta23,
00181                      const double dm221, const double dm243,
00182                      const double theta12,
00183                      const double theta13, const double theta14,
00184                      const double theta24, const double theta34,
00185                      const double delta1, const double delta2,
00186                      const double delta3,
00187                      const double baseline);
00188   void UnoscillateNDAverage(const double dm2, const double theta23,
00189                             const double dm221, const double dm243,
00190                             const double theta12,
00191                             const double theta13, const double theta14,
00192                             const double theta24, const double theta34,
00193                             const double delta1, const double delta2,
00194                             const double delta3,
00195                             const double baseline);
00196   
00197   virtual void InverseOscillate(Double_t dm2, Double_t sn2);
00198   virtual void InverseOscillateSimpleAverage(Double_t dm2, Double_t sn2);
00199   virtual void InverseOscillatePreAveraged(Double_t dm2, Double_t sn2);
00200   virtual void InverseOscillateLinearInterp(Double_t dm2, Double_t sn2);
00201   virtual void InverseOscillateHybrid(Double_t dm2, Double_t sn2);
00202   virtual void InverseOscillateNSI(Double_t dm2, Double_t sn2,
00203                                    Double_t epsilon, Double_t sign);
00204 
00205   // *** Rebin ***
00206   virtual void RebinToArray(std::vector<Double_t>& bins);
00207   virtual void RebinToArray(int nbins, const double binedges[999]);
00208   virtual void RebinToScheme(int scheme);
00209   virtual void RebinToScheme(NuBinningScheme::NuBinningScheme_t scheme);
00210   virtual void RebinForFit();
00211   virtual void RebinForPlotAsFit();
00212   virtual void RebinToHalfGeV(){RebinToGeV(0.5);}
00213   virtual void RebinForPlots();
00215   virtual void RebinTo4GeV(){RebinToGeV(4);}
00220   virtual void RebinToGeV(Double_t binwidth);
00221   virtual void RebinToNMB(){RebinToScheme(NuBinningScheme::kDisplayFD);}
00222 
00223   virtual void ResetSpectrum(TH1D& spectrum);
00224 
00225   NuMatrixSpectrum* Fluctuate();
00226   void FluctuateMe();
00227   
00229   virtual NuMatrixSpectrum *Copy(void) const { return new NuMatrixSpectrum(*this); };
00230 
00231   // TODO CJB - The NuMatrix1D version of this returns a TH1*
00232   TH1D* Spectrum() const {return fSpectrum;}
00233 
00234   // LEGACY FUNCTIONS
00235   // These are functions that exist because of entrenched code using NuMatrixspectrum
00236   // they only exist for NuMatrixSpectrum as people should move to the proper syntax
00237   Double_t PoT(void) const { return GetPOT(); }
00238   void ResetPoT(Double_t pot) { ResetPOT(pot); }
00239   void ScaleToPot(Double_t pot) { ScaleToPOT(pot); }
00240   
00241 protected:
00242   // These crash CINT for me when it builds the dictionary. Who can say why?
00243 #ifndef __CINT__
00244   TArrayD fPoiUp;
00245   TArrayD fPoiDn;
00246 #endif
00247 
00248   bool fDoPoisson;
00249 
00254   void FillOscCache(double dmsq) const;
00255   void FillOscCache_3flavour(Double_t dm232, Double_t t23, Double_t dm221, Double_t t12, Double_t t13,  Double_t delta, int inu) const;
00256   double threeFlavourWeight(double energy, double dm232, double t23, double dm221, double t12, double t13, double delta, int inu);
00257 
00259   static double fgCachedDmsq;
00261   static NuMatrixSpectrum* fgCachedOsc;
00262 
00263   ClassDef(NuMatrixSpectrum, 0);
00264 };
00265 
00266 #endif //NUMATRIXSPECTRUM_H