PWGDQ/dielectron/AliDielectronSignalBase.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 ///////////////////////////////////////////////////////////////////////////
  17 //                Dielectron SignalBase                                  //
  18 //                                                                       //
  19 //                                                                       //
  20 /*
  21 Base class for signal extraction from a histogram or an array of histograms
  22 The histogram is assumed to be an inv. mass spectrum,
  23 the array of histograms is assumed to be an array with inv. mass histograms
  24 resulting from single and mixed events, as defined in AliDielectron.cxx
  25
  26 */
  27 //                                                                       //
  28 ///////////////////////////////////////////////////////////////////////////
  29
  30
  31 #include <TVectorT.h>
  32 #include <TPaveText.h>
  33 #include <TF1.h>
  34 #include <TH1.h>
  35 #include <TDatabasePDG.h>
  36
  37 #include "AliDielectronSignalFunc.h"
  38 #include "AliDielectronSignalBase.h"
  39
  40 ClassImp(AliDielectronSignalBase)
  41
  42 TH1F* AliDielectronSignalBase::fgHistSimPM=0x0;
  43 TObject* AliDielectronSignalBase::fgPeakShape=0x0;
  44 const char* AliDielectronSignalBase::fgkValueNames[6] = {
  45   "Signal","Background","Significance","Signal/Background","Mass","MassWidth"};
  46
  47 AliDielectronSignalBase::AliDielectronSignalBase() :
  48   TNamed(),
  49   fHistSignal(0),
  50   fHistBackground(0),
  51   fHistDataPM(0),
  52   fHistDataPP(0),
  53   fHistDataMM(0),
  54   fHistDataME(0),
  55   fHistRfactor(0),
  56   fValues(6),
  57   fErrors(6),
  58   fIntMin(0),
  59   fIntMax(0),
  60   fFitMin(0),
  61   fFitMax(0),
  62   fRebin(1),
  63   fMethod(kLikeSign),
  64   fScaleMin(0.),
  65   fScaleMax(0.),
  66   fScaleMin2(0.),
  67   fScaleMax2(0.),
  68   fScaleFactor(1.),
  69   fMixingCorr(kFALSE),
  70   fPeakMethod(kBinCounting),
  71   fProcessed(kFALSE),
  72   fPOIpdg(443)
  73 {
  74   //
  75   // Default Constructor
  76   //
  77 }
  78
  79 //______________________________________________
  80 AliDielectronSignalBase::AliDielectronSignalBase(const char* name, const char* title) :
  81   TNamed(name, title),
  82   fHistSignal(0),
  83   fHistBackground(0),
  84   fHistDataPM(0),
  85   fHistDataPP(0),
  86   fHistDataMM(0),
  87   fHistDataME(0),
  88   fHistRfactor(0),
  89   fValues(6),
  90   fErrors(6),
  91   fIntMin(0),
  92   fIntMax(0),
  93   fFitMin(0),
  94   fFitMax(0),
  95   fRebin(1),
  96   fMethod(kLikeSign),
  97   fScaleMin(0.),
  98   fScaleMax(0.),
  99   fScaleMin2(0.),
 100   fScaleMax2(0.),
 101   fScaleFactor(1.),
 102   fMixingCorr(kFALSE),
 103   fPeakMethod(kBinCounting),
 104   fProcessed(kFALSE),
 105   fPOIpdg(443)
 106 {
 107   //
 108   // Named Constructor
 109   //
 110 }
 111
 112 //______________________________________________
 113 AliDielectronSignalBase::~AliDielectronSignalBase()
 114 {
 115   //
 116   // Default Destructor
 117   //
 118   if(fHistSignal) delete fHistSignal;
 119   if(fHistBackground) delete fHistBackground;
 120   if (fHistDataPP) delete fHistDataPP;
 121   if (fHistDataPM) delete fHistDataPM;
 122   if (fHistDataMM) delete fHistDataMM;
 123   if (fHistDataME) delete fHistDataME;
 124   if (fHistRfactor)delete fHistRfactor;
 125 }
 126
 127 //______________________________________________
 128 TPaveText* AliDielectronSignalBase::DrawStats(Double_t x1/*=0.*/, Double_t y1/*=0.*/, Double_t x2/*=0.*/, Double_t y2/*=0.*/)
 129 {
 130   //
 131   // Draw extracted values in a TPaveText
 132   // with the corners x1,y2,x2,y2
 133   //
 134   if (TMath::Abs(x1)<1e-20&&TMath::Abs(x2)<1e-20){
 135     x1=.6;
 136     x2=.9;
 137     y1=.7;
 138     y2=.9;
 139   }
 140   TPaveText *t=new TPaveText(x1,y1,x2,y2,"brNDC");
 141   t->SetFillColor(kWhite);
 142   t->SetBorderSize(1);
 143   t->SetTextAlign(12);
 144   t->AddText(Form("Range  : %.2f - %.2f GeV/c^{2}", fIntMin, fIntMax));
 145   t->AddText(Form("Signal : %.1f #pm %.1f", fValues(0), fErrors(0)));
 146   t->AddText(Form("Backgnd: %.1f #pm %.1f", fValues(1), fErrors(1)));
 147   t->AddText(Form("Signif.: %.2f #pm %.2f", fValues(2), fErrors(2)));
 148   t->AddText(Form("S/B    : %.2f #pm %.2f", fValues(3), fErrors(3)));
 149   if(fValues(4)>0)
 150     t->AddText(Form("Mass: %.2f #pm %.2f GeV/c^{2}", fValues(4), fErrors(4)));
 151   if(fValues(5)>0)
 152     t->AddText(Form("Mass res.: %.1f #pm %.1f MeV/c^{2}", 1000*fValues(5), 1000*fErrors(5)));
 153   t->Draw();
 154
 155   return t;
 156 }
 157
 158 //______________________________________________
 159 void AliDielectronSignalBase::Print(Option_t */*option*/) const
 160 {
 161   //
 162   // Print the statistics
 163   //
 164   printf("Signal : %.5g #pm %.5g\n",fValues(0), fErrors(0));
 165   printf("Backgnd: %.5g #pm %.5g\n",fValues(1), fErrors(1));
 166   printf("Signif.: %.5g #pm %.5g\n",fValues(2), fErrors(2));
 167   printf("SoB    : %.5g #pm %.5g\n",fValues(3), fErrors(3));
 168   if(fValues(4)>0)
 169     printf("Mass: %.5g #pm %.5g\n", fValues(4), fErrors(4));
 170   if(fValues(5)>0)
 171     printf("Mass res.: %.5g #pm %.5g\n", fValues(5), fErrors(5));
 172 }
 173
 174 //______________________________________________
 175 Double_t AliDielectronSignalBase::ScaleHistograms(TH1* histRaw, TH1* histBackground, Double_t intMin, Double_t intMax)
 176 {
 177   //
 178   // scale histBackground to match the integral of histRaw in the interval intMin, intMax
 179   //
 180
 181   //protect using over and underflow bins in normalisation calculation
 182   if (intMin<histRaw->GetXaxis()->GetXmin()) intMin=histRaw->GetXaxis()->GetXmin();
 183   if (intMin<histBackground->GetXaxis()->GetXmin()) intMin=histBackground->GetXaxis()->GetXmin();
 184
 185   if (intMax>histRaw->GetXaxis()->GetXmax())
 186     intMax=histRaw->GetXaxis()->GetXmax()-histRaw->GetBinWidth(histRaw->GetNbinsX())/2.;
 187   if (intMax>histBackground->GetXaxis()->GetXmax())
 188     intMax=histBackground->GetXaxis()->GetXmax()-histBackground->GetBinWidth(histBackground->GetNbinsX())/2.;
 189
 190   Double_t intRaw  = histRaw->Integral(histRaw->FindBin(intMin),histRaw->FindBin(intMax));
 191   Double_t intBack = histBackground->Integral(histBackground->FindBin(intMin),histBackground->FindBin(intMax));
 192   Double_t scaleFactor=intBack>0?intRaw/intBack:0.;
 193   if (intBack>0){
 194     histBackground->Sumw2();
 195     histBackground->Scale(scaleFactor);
 196   }
 197
 198   return scaleFactor;
 199 }
 200 //______________________________________________
 201 Double_t AliDielectronSignalBase::ScaleHistograms(TH1* histRaw, TH1* histBackground, Double_t intMin, Double_t intMax, Double_t intMin2, Double_t intMax2)
 202 {
 203   //
 204   // scale histBackground to match the integral of histRaw in the interval intMin, intMax and intMin2, intMax2
 205   //
 206
 207   if(intMin2==intMax2) return (ScaleHistograms(histRaw, histBackground, intMin, intMax));
 208
 209   //protect using over and underflow bins in normalisation calculation
 210   if (intMin<histRaw->GetXaxis()->GetXmin()) intMin=histRaw->GetXaxis()->GetXmin();
 211   if (intMin<histBackground->GetXaxis()->GetXmin()) intMin=histBackground->GetXaxis()->GetXmin();
 212
 213   if (intMax2>histRaw->GetXaxis()->GetXmax())
 214     intMax2=histRaw->GetXaxis()->GetXmax()-histRaw->GetBinWidth(histRaw->GetNbinsX())/2.;
 215   if (intMax2>histBackground->GetXaxis()->GetXmax())
 216     intMax2=histBackground->GetXaxis()->GetXmax()-histBackground->GetBinWidth(histBackground->GetNbinsX())/2.;
 217
 218   Double_t intRaw  = histRaw->Integral(histRaw->FindBin(intMin),histRaw->FindBin(intMax));
 219   Double_t intBack = histBackground->Integral(histBackground->FindBin(intMin),histBackground->FindBin(intMax));
 220   intRaw  += histRaw->Integral(histRaw->FindBin(intMin2),histRaw->FindBin(intMax2));
 221   intBack += histBackground->Integral(histBackground->FindBin(intMin2),histBackground->FindBin(intMax2));
 222
 223   Double_t scaleFactor=intBack>0?intRaw/intBack:0.;
 224   if (intBack>0){
 225     histBackground->Sumw2();
 226     histBackground->Scale(scaleFactor);
 227   }
 228
 229   return scaleFactor;
 230 }
 231
 232 TObject* AliDielectronSignalBase::DescribePeakShape(ESignalExtractionMethod method, Bool_t replaceValErr,  TH1F *mcShape) {
 233   //
 234   // Describe the extracted peak by the selected method and overwrite signal etc if needed
 235   //
 236
 237   fPeakMethod=method;
 238   Double_t data=0.;
 239   Double_t mc=0.;
 240   Double_t massPOI=TDatabasePDG::Instance()->GetParticle(fPOIpdg)->Mass();
 241   Double_t nPOI     = fHistSignal->GetBinContent(fHistSignal->FindBin(massPOI));
 242   Double_t binWidth = fHistSignal->GetBinWidth(  fHistSignal->FindBin(massPOI));
 243   TF1 *fit=0x0;
 244   Int_t parMass =-1;
 245   Int_t parSigma=-1;
 246
 247   // do the scaling/fitting
 248   switch(fPeakMethod) {
 249   case kBinCounting: /*nothing needs to be done*/
 250     break;
 251
 252   case kMCScaledMax:
 253     if(!mcShape) { printf(" ERROR: No MC histogram passed. Returning. \n"); return 0x0; }
 254     data = fHistSignal->GetBinContent(fHistSignal->FindBin(massPOI));
 255     mc   = mcShape->GetBinContent(fHistSignal->FindBin(massPOI));
 256     mcShape->Scale(data / mc );
 257     break;
 258
 259   case kMCScaledInt:
 260     if(!mcShape) { printf(" ERROR: No MC histogram passed. Returning. \n"); return 0x0; }
 261     if(mcShape->GetBinWidth(1)!=fHistSignal->GetBinWidth(1))
 262       printf(" WARNING: MC and signal histogram have different bin widths. \n");
 263     data = fHistSignal->Integral(fHistSignal->FindBin(fIntMin),fHistSignal->FindBin(fIntMax));
 264     mc   = mcShape->Integral(mcShape->FindBin(fIntMin),mcShape->FindBin(fIntMax));
 265     mcShape->Scale(data / mc );
 266     break;
 267
 268   case kMCFitted:
 269     if(!mcShape && !fgHistSimPM) { printf(" ERROR: No MC histogram passed or set. Returning. \n"); return 0x0; }
 270     if(!fgHistSimPM) fgHistSimPM=mcShape;
 271     fit = new TF1("fitMC",AliDielectronSignalFunc::PeakFunMC,fFitMin,fFitMax,1);
 272     fit->SetParNames("N");
 273     fHistSignal->Fit(fit,"RNI0");
 274     break;
 275
 276   case kCrystalBall:
 277     fit = new TF1("fitCB",AliDielectronSignalFunc::PeakFunCB,fFitMin,fFitMax,5);
 278     fit->SetParNames("alpha","n","meanx","sigma","N");
 279     //  fit->SetParameters(-.2,5.,gMjpsi,.06,20);
 280     //  fit->SetParameters(1.,3.6,gMjpsi,.08,700);
 281     fit->SetParameters(0.4, 4.0, massPOI, 0.025, 1.3*nPOI);
 282     fit->SetParLimits(0, 0.0,           1.           );
 283     fit->SetParLimits(1, 0.01,          10.          );
 284     fit->SetParLimits(2, massPOI-0.02,  massPOI+0.02 );
 285     fit->SetParLimits(3, 0.001,          0.2         );
 286     fit->SetParLimits(4, 0.2*nPOI,      2.0*nPOI     );
 287     parMass=2;
 288     parSigma=3;
 289     fHistSignal->Fit(fit,"RNI0");
 290     break;
 291
 292   case kGaus:
 293     fit = new TF1("fitGaus",AliDielectronSignalFunc::PeakFunGaus,fFitMin,fFitMax,3);
 294     //fit = new TF1("fitGaus","gaus",fFitMin,fFitMax);
 295     fit->SetParNames("N","meanx","sigma");
 296     fit->SetParameters(1.3*nPOI, massPOI, 0.025);
 297     fit->SetParLimits(0, 0.2*nPOI,      2.0*nPOI     );
 298     fit->SetParLimits(1, massPOI-0.02, massPOI+0.02);
 299     fit->SetParLimits(2, 0.001,         1.           );
 300     parMass=1;
 301     parSigma=2;
 302     fHistSignal->Fit(fit,"RNI0");
 303     break;
 304
 305   }
 306
 307   // overwrite values and errors if requested
 308   if(replaceValErr) {
 309     switch(fPeakMethod) {
 310     case kBinCounting:
 311       fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin), fHistSignal->FindBin(fIntMax), fErrors(0));
 312       SetSignificanceAndSOB();
 313       break;
 314     case kMCScaledMax:
 315     case kMCScaledInt:
 316       fValues(0) = mcShape->IntegralAndError(mcShape->FindBin(fIntMin), mcShape->FindBin(fIntMax), fErrors(0));
 317       SetSignificanceAndSOB();
 318       break;
 319
 320     case kMCFitted:
 321     case kCrystalBall:
 322     case kGaus:
 323       fValues(0) = fit->Integral(fIntMin, fIntMax)/binWidth;
 324       fErrors(0) = fit->IntegralError(fIntMin, fIntMax)/binWidth;
 325       SetSignificanceAndSOB();
 326       break;
 327     }
 328     // set mass position and width
 329     if(parMass>=0) {
 330       fValues(4) = fit->GetParameter(parMass);
 331       fErrors(4) = fit->GetParError(parMass);
 332     }
 333     if(parSigma>=0) {
 334       fValues(5) = fit->GetParameter(parSigma);
 335       fErrors(5) = fit->GetParError(parSigma);
 336     }
 337     else {
 338       // calculate FWHM
 339       SetFWHM();
 340     }
 341   }
 342
 343   // set the peak method obj
 344   switch(fPeakMethod) {
 345   case kBinCounting:
 346     fgPeakShape=(TH1F*)fHistSignal->Clone("BinCount");
 347     break;
 348   case kMCScaledMax:
 349   case kMCScaledInt:
 350     fgPeakShape=mcShape;
 351     break;
 352   case kMCFitted:
 353   case kCrystalBall:
 354   case kGaus:
 355     fgPeakShape=fit;
 356     break;
 357   }
 358
 359   return fgPeakShape;
 360
 361 }