Change Mult binning scheme
[u/mrichter/AliRoot.git] / 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 }