major dielectron update (included also the data and plotting macros for paper)
[u/mrichter/AliRoot.git] / PWG3 / dielectron / AliDielectronSignalFunc.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 SignalFunc                                  //
18 //                                                                       //
19 //                                                                       //
20 /*
21 Dielectron signal extraction class using functions as input.
22
23 A function to describe the signal as well as one to describe the background
24 has to be deployed by the user. Alternatively on of the default implementaions
25 can be used.
26
27 */
28 //                                                                       //
29 ///////////////////////////////////////////////////////////////////////////
30
31 #include <TF1.h>
32 #include <TH1.h>
33 #include <TGraph.h>
34 #include <TMath.h>
35 #include <TString.h>
36 #include <TPaveText.h>
37 #include <TList.h>
38 #include <TFitResult.h>
39 //#include <../hist/hist/src/TF1Helper.h>
40
41 #include <AliLog.h>
42
43 #include "AliDielectronSignalFunc.h"
44
45 ClassImp(AliDielectronSignalFunc)
46
47 AliDielectronSignalFunc::AliDielectronSignalFunc() :
48 AliDielectronSignalBase(),
49 fFuncSignal(0x0),
50 fFuncBackground(0x0),
51 fFuncSigBack(0x0),
52 fParMass(1),
53 fParMassWidth(2),
54 fFitOpt("SMNQE"),
55 fUseIntegral(kFALSE)
56 {
57   //
58   // Default Constructor
59   //
60   
61 }
62
63 //______________________________________________
64 AliDielectronSignalFunc::AliDielectronSignalFunc(const char* name, const char* title) :
65 AliDielectronSignalBase(name, title),
66 fFuncSignal(0x0),
67 fFuncBackground(0x0),
68 fFuncSigBack(0x0),
69 fParMass(1),
70 fParMassWidth(2),
71 fFitOpt("SMNQE"),
72 fUseIntegral(kFALSE)
73 {
74   //
75   // Named Constructor
76   //
77 }
78
79 //______________________________________________
80 AliDielectronSignalFunc::~AliDielectronSignalFunc()
81 {
82   //
83   // Default Destructor
84   //
85   if(fFuncSignal) delete fFuncSignal;
86   if(fFuncBackground) delete fFuncBackground;
87   if(fFuncSigBack) delete fFuncSigBack;
88 }
89
90
91 //______________________________________________
92 void AliDielectronSignalFunc::Process(TObjArray * const arrhist)
93 {
94   //
95   // Fit the invariant mass histograms and retrieve the signal and background
96   //
97   switch(fMethod) {
98   case kFitted :
99     ProcessFit(arrhist);
100     break;
101     
102   case kLikeSign :
103     ProcessLS(arrhist);
104     break;
105     
106   case kEventMixing :
107     ProcessEM(arrhist);
108     break;
109     
110   default :
111     AliError("Background substraction method not known!");
112   }
113 }
114
115
116 //______________________________________________
117 void AliDielectronSignalFunc::ProcessFit(TObjArray * const arrhist) {
118   //
119   // Fit the +- invariant mass distribution only
120   // Here we assume that the combined fit function is a sum of the signal and background functions
121   //    and that the signal function is always the first term of this sum
122   //
123   
124   fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
125   fHistDataPM->Sumw2();
126   if(fRebin>1)
127     fHistDataPM->Rebin(fRebin);
128   
129   fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal",
130                          fHistDataPM->GetXaxis()->GetNbins(),
131                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
132   fHistBackground = new TH1F("HistBackground", "Like-sign contribution",
133                              fHistDataPM->GetXaxis()->GetNbins(),
134                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
135   
136   // the starting parameters of the fit function and their limits can be tuned
137   // by the user in its macro
138   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
139   TFitResultPtr pmFitPtr = fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
140   //TFitResult *pmFitResult = pmFitPtr.Get(); // used only with TF1Helper
141   fFuncSignal->SetParameters(fFuncSigBack->GetParameters());
142   fFuncBackground->SetParameters(fFuncSigBack->GetParameters()+fFuncSignal->GetNpar());
143   
144   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
145     Double_t m = fHistDataPM->GetBinCenter(iBin);
146     Double_t pm = fHistDataPM->GetBinContent(iBin);
147     Double_t epm = fHistDataPM->GetBinError(iBin);
148     Double_t bknd = fFuncBackground->Eval(m);
149     Double_t ebknd = 0;
150     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
151 /* TF1Helper problem on alien compilation
152       for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
153         TF1 gradientIpar("gradientIpar",
154                          ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
155         TF1 gradientJpar("gradientJpar",
156                          ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
157         ebknd += pmFitResult->CovMatrix(iPar,jPar)*
158           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
159           (iPar==jPar ? 1.0 : 2.0);
160       }
161 */
162     }
163     Double_t signal = pm-bknd;
164     Double_t error = TMath::Sqrt(epm*epm+ebknd);
165     fHistSignal->SetBinContent(iBin, signal);
166     fHistSignal->SetBinError(iBin, error);
167     fHistBackground->SetBinContent(iBin, bknd);
168     fHistBackground->SetBinError(iBin, TMath::Sqrt(ebknd));
169   }
170   
171   if(fUseIntegral) {
172     // signal
173     fValues(0) = fFuncSignal->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
174     fErrors(0) = 0;
175     for(Int_t iPar=0; iPar<fFuncSignal->GetNpar(); iPar++) {
176 /* TF1Helper problem on alien compilation
177       for(Int_t jPar=iPar; jPar<fFuncSignal->GetNpar(); jPar++) {
178         TF1 gradientIpar("gradientIpar",
179                          ROOT::TF1Helper::TGradientParFunction(iPar,fFuncSignal),0,0,0);
180         TF1 gradientJpar("gradientJpar",
181                          ROOT::TF1Helper::TGradientParFunction(jPar,fFuncSignal),0,0,0);
182         fErrors(0) += pmFitResult->CovMatrix(iPar,jPar)*
183           gradientIpar.Integral(fIntMin,fIntMax)*gradientJpar.Integral(fIntMin,fIntMax)*
184           (iPar==jPar ? 1.0 : 2.0);
185       }
186 */
187     }
188     // background
189     fValues(1) = fFuncBackground->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
190     fErrors(1) = 0;
191     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
192 /* TF1Helper problem on alien compilation
193       for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
194         TF1 gradientIpar("gradientIpar",
195                          ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
196         TF1 gradientJpar("gradientJpar",
197                          ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
198         fErrors(1) += pmFitResult->CovMatrix(iPar,jPar)*
199           gradientIpar.Integral(fIntMin, fIntMax)*gradientJpar.Integral(fIntMin, fIntMax)*
200           (iPar==jPar ? 1.0 : 2.0);
201       }
202 */
203     }
204   }
205   else {
206     // signal
207     fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
208                                                fHistSignal->FindBin(fIntMax), fErrors(0));
209     // background
210     fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
211       fHistBackground->FindBin(fIntMax),
212       fErrors(1));
213   }
214   // S/B and significance
215   SetSignificanceAndSOB();
216   fValues(4) = fFuncSigBack->GetParameter(fParMass);
217   fErrors(4) = fFuncSigBack->GetParError(fParMass);
218   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
219   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
220   
221   fProcessed = kTRUE;
222 }
223
224 //______________________________________________
225 void AliDielectronSignalFunc::ProcessLS(TObjArray * const arrhist) {
226   //
227   // Substract background using the like-sign spectrum
228   //
229   fHistDataPP = (TH1F*)(arrhist->At(0))->Clone("histPP");  // ++    SE
230   fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
231   fHistDataMM = (TH1F*)(arrhist->At(2))->Clone("histMM");  // --    SE
232   if (fRebin>1) {
233     fHistDataPP->Rebin(fRebin);
234     fHistDataPM->Rebin(fRebin);
235     fHistDataMM->Rebin(fRebin);
236   }
237   fHistDataPP->Sumw2();
238   fHistDataPM->Sumw2();
239   fHistDataMM->Sumw2();
240   
241   fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal",
242                          fHistDataPM->GetXaxis()->GetNbins(),
243                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
244   fHistBackground = new TH1F("HistBackground", "Like-sign contribution",
245                              fHistDataPM->GetXaxis()->GetNbins(),
246                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
247   
248   // fit the +- mass distribution
249   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
250   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
251   // declare the variables where the like-sign fit results will be stored
252 //   TFitResult *ppFitResult = 0x0;
253 //   TFitResult *mmFitResult = 0x0;
254   // fit the like sign background
255   TF1 *funcClonePP = (TF1*)fFuncBackground->Clone("funcClonePP");
256   TF1 *funcCloneMM = (TF1*)fFuncBackground->Clone("funcCloneMM");
257   fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
258   fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
259 //   TFitResultPtr ppFitPtr = fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
260 //   ppFitResult = ppFitPtr.Get();
261   fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
262   fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
263 //   TFitResultPtr mmFitPtr = fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
264 //   mmFitResult = mmFitPtr.Get();
265   
266   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
267     Double_t m = fHistDataPM->GetBinCenter(iBin);
268     Double_t pm = fHistDataPM->GetBinContent(iBin);
269     Double_t pp = funcClonePP->Eval(m);
270     Double_t mm = funcCloneMM->Eval(m);
271     Double_t epm = fHistDataPM->GetBinError(iBin);
272     Double_t epp = 0;
273     for(Int_t iPar=0; iPar<funcClonePP->GetNpar(); iPar++) {
274 /* TF1Helper problem on alien compilation
275       for(Int_t jPar=iPar; jPar<funcClonePP->GetNpar(); jPar++) {
276         TF1 gradientIpar("gradientIpar",
277                          ROOT::TF1Helper::TGradientParFunction(iPar,funcClonePP),0,0,0);
278         TF1 gradientJpar("gradientJpar",
279                          ROOT::TF1Helper::TGradientParFunction(jPar,funcClonePP),0,0,0);
280         epp += ppFitResult->CovMatrix(iPar,jPar)*
281           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
282           (iPar==jPar ? 1.0 : 2.0);
283       }
284 */
285     }
286     Double_t emm = 0;
287     for(Int_t iPar=0; iPar<funcCloneMM->GetNpar(); iPar++) {
288 /* TF1Helper problem on alien compilation
289       for(Int_t jPar=iPar; jPar<funcCloneMM->GetNpar(); jPar++) {
290         TF1 gradientIpar("gradientIpar",
291                          ROOT::TF1Helper::TGradientParFunction(iPar,funcCloneMM),0,0,0);
292         TF1 gradientJpar("gradientJpar",
293                          ROOT::TF1Helper::TGradientParFunction(jPar,funcCloneMM),0,0,0);
294         emm += mmFitResult->CovMatrix(iPar,jPar)*
295           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
296           (iPar==jPar ? 1.0 : 2.0);
297       }
298 */
299     }
300     
301     Double_t signal = pm-2.0*TMath::Sqrt(pp*mm);
302     Double_t background = 2.0*TMath::Sqrt(pp*mm);
303     // error propagation on the signal calculation above
304     Double_t esignal = TMath::Sqrt(epm*epm+(mm/pp)*epp+(pp/mm)*emm);
305     Double_t ebackground = TMath::Sqrt((mm/pp)*epp+(pp/mm)*emm);
306     fHistSignal->SetBinContent(iBin, signal);
307     fHistSignal->SetBinError(iBin, esignal);
308     fHistBackground->SetBinContent(iBin, background);
309     fHistBackground->SetBinError(iBin, ebackground);
310   }
311   
312   // signal
313   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
314                                              fHistSignal->FindBin(fIntMax), fErrors(0));
315   // background
316   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
317                                                  fHistBackground->FindBin(fIntMax),
318                                                  fErrors(1));
319   // S/B and significance
320   SetSignificanceAndSOB();
321   fValues(4) = fFuncSigBack->GetParameter(fParMass);
322   fErrors(4) = fFuncSigBack->GetParError(fParMass);
323   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
324   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
325   
326   fProcessed = kTRUE;
327 }
328
329 //______________________________________________
330 void AliDielectronSignalFunc::ProcessEM(TObjArray * const arrhist) {
331   //
332   // Substract background with the event mixing technique
333   //
334   arrhist->GetEntries();   // just to avoid the unused parameter warning
335   AliError("Event mixing for background substraction method not implemented!");
336 }
337
338 //______________________________________________
339 void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig, TF1 * const back,
340                                            Int_t parM, Int_t parMres)
341 {
342   //
343   // Set the signal, background functions and combined fit function
344   // Note: The process method assumes that the first n parameters in the
345   //       combined fit function correspond to the n parameters of the signal function
346   //       and the n+1 to n+m parameters to the m parameters of the background function!!!
347   
348   if (!sig||!back||!combined) {
349     AliError("Both, signal and background function need to be set!");
350     return;
351   }
352   fFuncSignal=sig;
353   fFuncBackground=back;
354   fFuncSigBack=combined;
355   fParMass=parM;
356   fParMassWidth=parMres;
357 }
358
359 //______________________________________________
360 void AliDielectronSignalFunc::SetDefaults(Int_t type)
361 {
362   //
363   // Setup some default functions:
364   // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
365   // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
366   // type = 2: half gaussian, half exponential signal function
367   // type = 3: Crystal-Ball function
368   // type = 4: Crystal-Ball signal + exponential background
369   //
370   
371   if (type==0){
372     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
373     fFuncBackground=new TF1("DieleBackground","pol1",2.5,4);
374     fFuncSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
375     
376     fFuncSigBack->SetParameters(1,3.1,.05,2.5,1);
377     fFuncSigBack->SetParLimits(0,0,10000000);
378     fFuncSigBack->SetParLimits(1,3.05,3.15);
379     fFuncSigBack->SetParLimits(2,.02,.1);
380   }
381   else if (type==1){
382     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
383     fFuncBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
384     fFuncSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
385     
386     fFuncSigBack->SetParameters(1,3.1,.05,1,2.5,1);
387     fFuncSigBack->SetParLimits(0,0,10000000);
388     fFuncSigBack->SetParLimits(1,3.05,3.15);
389     fFuncSigBack->SetParLimits(2,.02,.1);
390   }
391   else if (type==2){
392     // half gaussian, half exponential signal function
393     // exponential background
394     fFuncSignal = new TF1("DieleSignal","(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))",2.5,4);
395     fFuncBackground = new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
396     fFuncSigBack = new TF1("DieleCombined","(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))+[4]*exp(-(x-[5])/[6])+[7]",2.5,4);
397     fFuncSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
398     
399     fFuncSigBack->SetParLimits(0,0,10000000);
400     fFuncSigBack->SetParLimits(1,3.05,3.15);
401     fFuncSigBack->SetParLimits(2,.02,.1);
402     fFuncSigBack->FixParameter(6,2.5);
403     fFuncSigBack->FixParameter(7,0);
404   }
405 }
406
407
408 //______________________________________________
409 void AliDielectronSignalFunc::Draw(const Option_t* option)
410 {
411   //
412   // Draw the fitted function
413   //
414   
415   TString drawOpt(option);
416   drawOpt.ToLower();
417   
418   Bool_t optStat=drawOpt.Contains("stat");
419   
420   fFuncSigBack->SetNpx(200);
421   fFuncSigBack->SetRange(fIntMin,fIntMax);
422   fFuncBackground->SetNpx(200);
423   fFuncBackground->SetRange(fIntMin,fIntMax);
424   
425   TGraph *grSig=new TGraph(fFuncSigBack);
426   grSig->SetFillColor(kGreen);
427   grSig->SetFillStyle(3001);
428   
429   TGraph *grBack=new TGraph(fFuncBackground);
430   grBack->SetFillColor(kRed);
431   grBack->SetFillStyle(3001);
432   
433   grSig->SetPoint(0,grBack->GetX()[0],grBack->GetY()[0]);
434   grSig->SetPoint(grSig->GetN()-1,grBack->GetX()[grBack->GetN()-1],grBack->GetY()[grBack->GetN()-1]);
435   
436   grBack->SetPoint(0,grBack->GetX()[0],0.);
437   grBack->SetPoint(grBack->GetN()-1,grBack->GetX()[grBack->GetN()-1],0.);
438   
439   fFuncSigBack->SetRange(fFitMin,fFitMax);
440   fFuncBackground->SetRange(fFitMin,fFitMax);
441   
442   if (!drawOpt.Contains("same")){
443     if (fHistDataPM){
444       fHistDataPM->Draw();
445       grSig->Draw("f");
446     } else {
447       grSig->Draw("af");
448     }
449   } else {
450     grSig->Draw("f");
451   }
452   if(fMethod==kFitted) grBack->Draw("f");
453   fFuncSigBack->Draw("same");
454   fFuncSigBack->SetLineWidth(2);
455   if(fMethod==kLikeSign) {
456     fHistDataPP->SetLineWidth(2);
457     fHistDataPP->SetLineColor(6);
458     fHistDataPP->Draw("same");
459     fHistDataMM->SetLineWidth(2);
460     fHistDataMM->SetLineColor(8);
461     fHistDataMM->Draw("same");
462   }
463   
464   if(fMethod==kFitted)
465     fFuncBackground->Draw("same");
466   
467   if (optStat) DrawStats();
468   
469 }