PWGDQ/dielectron/AliDielectronSignalBase.h

   1 #ifndef ALIDIELECTRONSIGNALBASE_H
   2 #define ALIDIELECTRONSIGNALBASE_H
   3
   4 /* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   5  * See cxx source for full Copyright notice                               */
   6
   7 //#############################################################
   8 //#                                                           #
   9 //#         Class AliDielectronSignalBase                       #
  10 //#         Manage Cuts on the legs of the pair               #
  11 //#                                                           #
  12 //#  Authors:                                                 #
  13 //#   Anton     Andronic, GSI / A.Andronic@gsi.de             #
  14 //#   Ionut C.  Arsene,   GSI / I.C.Arsene@gsi.de             #
  15 //#   Julian    Book,     Uni Ffm / Julian.Book@cern.ch       #
  16 //#   Frederick Kramer,   Uni Ffm, / Frederick.Kramer@cern.ch #
  17 //#   Magnus    Mager,    CERN / Magnus.Mager@cern.ch         #
  18 //#   WooJin J. Park,     GSI / W.J.Park@gsi.de               #
  19 //#   Jens      Wiechula, Uni HD / Jens.Wiechula@cern.ch      #
  20 //#                                                           #
  21 //#############################################################
  22
  23
  24 #include <TNamed.h>
  25 #include <TVectorT.h>
  26 #include <TMath.h>
  27 #include <TH1F.h>
  28 #include <TF1.h>
  29
  30 class TObjArray;
  31 class TPaveText;
  32
  33 class AliDielectronSignalBase : public TNamed {
  34 public:
  35   enum EBackgroundMethod {
  36     kFittedMC = 0,
  37     kFitted,
  38     kLikeSign,
  39     kLikeSignArithm,
  40     kLikeSignRcorr,
  41     kLikeSignArithmRcorr,
  42     kLikeSignFit,
  43     kEventMixing,
  44     kEventMixingFit,
  45     kRotation
  46   };
  47
  48   enum ESignalExtractionMethod {
  49     kBinCounting = 0,
  50     kMCScaledMax,
  51     kMCScaledInt,
  52     kMCFitted,
  53     kCrystalBall,
  54     kGaus
  55   };
  56
  57   AliDielectronSignalBase();
  58   AliDielectronSignalBase(const char*name, const char* title);
  59
  60   virtual ~AliDielectronSignalBase();
  61
  62   void SetMCSignalShape(TH1F* hist) { fgHistSimPM=hist; }
  63   void SetParticleOfInterest(Int_t pdgcode) { fPOIpdg=pdgcode; }
  64   void SetIntegralRange(Double_t min, Double_t max) {fIntMin=min;fIntMax=max;}
  65   void SetFitRange(Double_t min, Double_t max) {fFitMin=min; fFitMax=max;}
  66   void SetRebin(Int_t factor) {fRebin = factor;}
  67   void SetMethod(EBackgroundMethod method) {fMethod = method;}
  68   Int_t GetMethod() const { return (Int_t)fMethod; }
  69   void SetExtractionMethod(ESignalExtractionMethod method) {fPeakMethod = method;}
  70   Int_t GetExtractionMethod() const { return (Int_t)fPeakMethod; }
  71
  72   const TVectorD& GetValues() const {return fValues;}
  73   const TVectorD& GetErrors() const {return fErrors;}
  74
  75   Double_t GetIntegralMin()          const { return fIntMin; }
  76   Double_t GetIntegralMax()          const { return fIntMax; }
  77   Double_t GetSignal()               const { return fValues(0);}
  78   Double_t GetSignalError()          const { return fErrors(0);}
  79   Double_t GetBackground()           const { return fValues(1);}
  80   Double_t GetBackgroundError()      const { return fErrors(1);}
  81   Double_t GetSignificance()         const { return fValues(2);}
  82   Double_t GetSignificanceError()    const { return fErrors(2);}
  83   Double_t GetSB()                   const { return fValues(3);}
  84   Double_t GetSBError()              const { return fErrors(3);}
  85   Double_t GetMass()                 const { return fValues(4);}
  86   Double_t GetMassError()            const { return fErrors(4);}
  87   Double_t GetMassWidth()            const { return fValues(5);}
  88   Double_t GetMassWidthError()       const { return fErrors(5);}
  89   static const char* GetValueName(Int_t i) { return (i>=0&&i<6)?fgkValueNames[i]:""; }
  90
  91   TH1* GetSignalHistogram()      const {return fHistSignal;}
  92   TH1* GetBackgroundHistogram()  const {return fHistBackground;}
  93   TH1* GetUnlikeSignHistogram()  const {return fHistDataPM;}
  94   TH1* GetRfactorHistogram()     const {return fHistRfactor;}
  95   TObject* GetPeakShape()        const {return fgPeakShape;}
  96
  97   void SetScaleRawToBackground(Double_t intMin, Double_t intMax) { fScaleMin=intMin; fScaleMax=intMax; }
  98   void SetScaleRawToBackground(Double_t intMin, Double_t intMax, Double_t intMin2, Double_t intMax2) { fScaleMin=intMin; fScaleMax=intMax; fScaleMin2=intMin2; fScaleMax2=intMax2; }
  99   Double_t GetScaleMin()        const { return fScaleMin;    }
 100   Double_t GetScaleMax()        const { return fScaleMax;    }
 101   Double_t GetScaleMin2()       const { return fScaleMin2;   }
 102   Double_t GetScaleMax2()       const { return fScaleMax2;   }
 103   Double_t GetScaleFactor()     const { return fScaleFactor; }
 104   Int_t GetParticleOfInterest() const { return fPOIpdg;      }
 105
 106   void SetMixingCorrection(Bool_t mixcorr=kTRUE) { fMixingCorr=mixcorr; }
 107
 108   static Double_t ScaleHistograms(TH1* histRaw, TH1* histBackground, Double_t intMin, Double_t intMax);
 109   static Double_t ScaleHistograms(TH1* histRaw, TH1* histBackground, Double_t intMin, Double_t intMax, Double_t intMin2, Double_t intMax2);
 110
 111   virtual void Print(Option_t *option="") const;
 112
 113   /**
 114   This function needs to be implemented by the signal extraction classes.
 115   Here all the work should be done.
 116
 117   The signal extraction is done on the mass spectra.
 118   The TObjArray should contain the Inv. Mass spectra of the 10 possible combinations
 119      for single and mixed events defined in AliDielectron.cxx
 120   */
 121   virtual void Process(TObjArray * const /*arrhist*/) = 0;
 122   TObject* DescribePeakShape(ESignalExtractionMethod method, Bool_t replaceValErr=kFALSE,  TH1F *mcShape=0x0);
 123
 124 protected:
 125
 126   TH1 *fHistSignal;                  // histogram of pure signal
 127   TH1 *fHistBackground;              // histogram of background (fitted=0, like-sign=1, event mixing=2)
 128   TH1 *fHistDataPM;                  // histogram of selected +- pair candidates
 129   TH1 *fHistDataPP;                  // histogram of selected ++ pair candidates
 130   TH1 *fHistDataMM;                  // histogram of selected -- pair candidates
 131   TH1 *fHistDataME;                  // histogram of selected +- pair candidates from mixed event
 132   TH1 *fHistRfactor;                 // histogram of R factors
 133
 134   TVectorD fValues;                  // values
 135   TVectorD fErrors;                  // value errors
 136
 137   Double_t fIntMin;                  // signal extraction range min
 138   Double_t fIntMax;                  // signal extraction range max
 139   Double_t fFitMin;                  // fit range lowest inv. mass
 140   Double_t fFitMax;                  // fit range highest inv. mass
 141
 142   Int_t fRebin;                      // histogram rebin factor
 143   EBackgroundMethod fMethod;         // method for background substraction
 144   Double_t fScaleMin;                // min for scaling of raw and background histogram
 145   Double_t fScaleMax;                // max for scaling of raw and background histogram
 146   Double_t fScaleMin2;                // min for scaling of raw and background histogram
 147   Double_t fScaleMax2;                // max for scaling of raw and background histogram
 148   Double_t fScaleFactor;             // scale factor of raw to background histogram scaling
 149   Bool_t fMixingCorr;                // switch for bin by bin correction with R factor
 150
 151   ESignalExtractionMethod fPeakMethod;  // method for peak description and signal extraction
 152   static TObject *fgPeakShape;       // histogram or function used to describe the extracted signal
 153
 154   Bool_t fProcessed;                 // flag
 155   Int_t  fPOIpdg;                    // pdg code particle of interest
 156   static TH1F* fgHistSimPM;          // simulated peak shape
 157
 158   void SetSignificanceAndSOB();      // calculate the significance and S/B
 159   void SetFWHM();                    // calculate the fwhm
 160   static const char* fgkValueNames[6];  //value names
 161
 162   TPaveText* DrawStats(Double_t x1=0., Double_t y1=0., Double_t x2=0., Double_t y2=0.);
 163
 164   AliDielectronSignalBase(const AliDielectronSignalBase &c);
 165   AliDielectronSignalBase &operator=(const AliDielectronSignalBase &c);
 166
 167   ClassDef(AliDielectronSignalBase,4)         // base and abstract class for signal extraction
 168 };
 169
 170 inline void AliDielectronSignalBase::SetSignificanceAndSOB()
 171 {
 172   //
 173   // Calculate S/B and significance
 174   //
 175   // Signal/Background
 176   fValues(3) = (fValues(1)>0 ? fValues(0)/fValues(1) : 0);
 177   Float_t epsSig = (fValues(0)>0 ? fErrors(0)/fValues(0) : 0);
 178   Float_t epsBknd = (fValues(1)>0 ? fErrors(1)/fValues(1) : 0);
 179   fErrors(3) = fValues(3)*TMath::Sqrt(epsSig*epsSig + epsBknd*epsBknd);
 180   // Significance
 181   fValues(2) = ((fValues(0)+fValues(1))>0 ? fValues(0)/TMath::Sqrt(fValues(0)+fValues(1)) : 0);
 182   Float_t s = (fValues(0)>0?fValues(0):0); Float_t b = fValues(1);
 183   Float_t se = fErrors(0); Float_t be = fErrors(1);
 184   // fErrors(2) = ((s+b)>0 ? TMath::Sqrt((s*(s+2*b)*(s+2*b)+b*s*s)/(4*TMath::Power(s+b,3))) : 0); // old implementation
 185   fErrors(2) = ((s+b)>0 ? fValues(2)*TMath::Sqrt(be*be + TMath::Power(se*(s+2*b)/s, 2)) / 2 / (s+b) : 0);
 186 }
 187
 188 inline void AliDielectronSignalBase::SetFWHM()
 189 {
 190   // calculate the fwhm
 191   if(!fgPeakShape) return;
 192
 193   // case for TF1
 194   if(fgPeakShape->IsA() == TF1::Class()) {
 195     TF1* fit  = (TF1*) fgPeakShape->Clone("fit");
 196     TF1* pfit = (TF1*) fit->Clone("pfit");
 197     TF1* mfit = (TF1*) fit->Clone("mfit");
 198     for (Int_t i=0; i<fit->GetNpar(); i++) {
 199       pfit->SetParameter(i,fit->GetParameter(i) + fit->GetParError(i));
 200       mfit->SetParameter(i,fit->GetParameter(i) - fit->GetParError(i));
 201     }
 202     Double_t maxX   = fit->GetMaximumX();
 203     Double_t maxY   = fit->GetHistogram()->GetMaximum();
 204     Double_t xAxMin = fit->GetXmin();
 205     Double_t xAxMax = fit->GetXmax();
 206     // fwhms
 207     Double_t fwhmMin  = fit->GetX(.5*maxY, xAxMin, maxX);
 208     Double_t fwhmMax  = fit->GetX(.5*maxY, maxX, xAxMax);
 209     Double_t pfwhmMin = pfit->GetX(.5*maxY, xAxMin, maxX);
 210     Double_t pfwhmMax = pfit->GetX(.5*maxY, maxX, xAxMax);
 211     Double_t mfwhmMin = mfit->GetX(.5*maxY, xAxMin, maxX);
 212     Double_t mfwhmMax = mfit->GetX(.5*maxY, maxX, xAxMax);
 213     Double_t pError = TMath::Abs( (fwhmMax-fwhmMin) - (pfwhmMax-pfwhmMin) );
 214     Double_t mError = TMath::Abs( (fwhmMax-fwhmMin) - (mfwhmMax-mfwhmMin) );
 215     fValues(5) = (fwhmMax-fwhmMin);
 216     fErrors(5) = (pError>=mError ? pError : mError);
 217     delete fit;
 218     delete pfit;
 219     delete mfit;
 220   }
 221   else if(fgPeakShape->IsA() == TH1F::Class()) {
 222     // th1 calculation
 223     TH1F *hist = (TH1F*) fgPeakShape->Clone("hist");
 224     Int_t bin1 = hist->FindFirstBinAbove(hist->GetMaximum()/2);
 225     Int_t bin2 = hist->FindLastBinAbove(hist->GetMaximum()/2);
 226     fValues(5) = hist->GetBinCenter(bin2) - hist->GetBinCenter(bin1);
 227     fErrors(5) = 0.0; // not defined
 228     delete hist;
 229   }
 230
 231 }
 232
 233 #endif