small fix to remove TF1Helper which sneaked in again, inadvertently...
[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> //not supposed to be used!
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(); //not used when TF1Helper out
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     /* to be revised ... TF1Helper
151     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
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     */ // TF1Helper
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     /* to be revised ... TF1Helper
176     for(Int_t iPar=0; iPar<fFuncSignal->GetNpar(); iPar++) {
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     */ //TF1Helper
188     // background
189     fValues(1) = fFuncBackground->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
190     fErrors(1) = 0;
191     /* to be revised... TF1Helper
192     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
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     */ // TF1Helper
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   TFitResultPtr ppFitPtr = fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
259   ppFitResult = ppFitPtr.Get();
260   fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
261   TFitResultPtr mmFitPtr = fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
262   mmFitResult = mmFitPtr.Get();
263   
264   /* to be revised ... TF1Helper
265   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
266     Double_t m = fHistDataPM->GetBinCenter(iBin);
267     Double_t pm = fHistDataPM->GetBinContent(iBin);
268     Double_t pp = funcClonePP->Eval(m);
269     Double_t mm = funcCloneMM->Eval(m);
270     Double_t epm = fHistDataPM->GetBinError(iBin);
271     Double_t epp = 0;
272     for(Int_t iPar=0; iPar<funcClonePP->GetNpar(); iPar++) {
273       for(Int_t jPar=iPar; jPar<funcClonePP->GetNpar(); jPar++) {
274         TF1 gradientIpar("gradientIpar",
275                          ROOT::TF1Helper::TGradientParFunction(iPar,funcClonePP),0,0,0);
276         TF1 gradientJpar("gradientJpar",
277                          ROOT::TF1Helper::TGradientParFunction(jPar,funcClonePP),0,0,0);
278         epp += ppFitResult->CovMatrix(iPar,jPar)*
279           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
280           (iPar==jPar ? 1.0 : 2.0);
281       }
282     }
283     
284     Double_t emm = 0;
285     for(Int_t iPar=0; iPar<funcCloneMM->GetNpar(); iPar++) {
286       for(Int_t jPar=iPar; jPar<funcCloneMM->GetNpar(); jPar++) {
287         TF1 gradientIpar("gradientIpar",
288                          ROOT::TF1Helper::TGradientParFunction(iPar,funcCloneMM),0,0,0);
289         TF1 gradientJpar("gradientJpar",
290                          ROOT::TF1Helper::TGradientParFunction(jPar,funcCloneMM),0,0,0);
291         emm += mmFitResult->CovMatrix(iPar,jPar)*
292           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
293           (iPar==jPar ? 1.0 : 2.0);
294       }
295     }
296     
297     Double_t signal = pm-2.0*TMath::Sqrt(pp*mm);
298     Double_t background = 2.0*TMath::Sqrt(pp*mm);
299     // error propagation on the signal calculation above
300     Double_t esignal = TMath::Sqrt(epm*epm+(mm/pp)*epp+(pp/mm)*emm);
301     Double_t ebackground = TMath::Sqrt((mm/pp)*epp+(pp/mm)*emm);
302     fHistSignal->SetBinContent(iBin, signal);
303     fHistSignal->SetBinError(iBin, esignal);
304     fHistBackground->SetBinContent(iBin, background);
305     fHistBackground->SetBinError(iBin, ebackground);
306   }
307   */  // TF1Helper
308
309   // signal
310   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
311                                              fHistSignal->FindBin(fIntMax), fErrors(0));
312   // background
313   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
314                                                  fHistBackground->FindBin(fIntMax),
315                                                  fErrors(1));
316   // S/B and significance
317   SetSignificanceAndSOB();
318   fValues(4) = fFuncSigBack->GetParameter(fParMass);
319   fErrors(4) = fFuncSigBack->GetParError(fParMass);
320   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
321   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
322   
323   fProcessed = kTRUE;
324 }
325
326 //______________________________________________
327 void AliDielectronSignalFunc::ProcessEM(TObjArray * const arrhist) {
328   //
329   // Substract background with the event mixing technique
330   //
331   arrhist->GetEntries();   // just to avoid the unused parameter warning
332   AliError("Event mixing for background substraction method not implemented!");
333 }
334
335 //______________________________________________
336 void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig, TF1 * const back,
337                                            Int_t parM, Int_t parMres)
338 {
339   //
340   // Set the signal, background functions and combined fit function
341   // Note: The process method assumes that the first n parameters in the
342   //       combined fit function correspond to the n parameters of the signal function
343   //       and the n+1 to n+m parameters to the m parameters of the background function!!!
344   
345   if (!sig||!back||!combined) {
346     AliError("Both, signal and background function need to be set!");
347     return;
348   }
349   fFuncSignal=sig;
350   fFuncBackground=back;
351   fFuncSigBack=combined;
352   fParMass=parM;
353   fParMassWidth=parMres;
354 }
355
356 //______________________________________________
357 void AliDielectronSignalFunc::SetDefaults(Int_t type)
358 {
359   //
360   // Setup some default functions:
361   // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
362   // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
363   // type = 2: half gaussian, half exponential signal function
364   // type = 3: Crystal-Ball function
365   // type = 4: Crystal-Ball signal + exponential background
366   //
367   
368   if (type==0){
369     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
370     fFuncBackground=new TF1("DieleBackground","pol1",2.5,4);
371     fFuncSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
372     
373     fFuncSigBack->SetParameters(1,3.1,.05,2.5,1);
374     fFuncSigBack->SetParLimits(0,0,10000000);
375     fFuncSigBack->SetParLimits(1,3.05,3.15);
376     fFuncSigBack->SetParLimits(2,.02,.1);
377   }
378   else if (type==1){
379     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
380     fFuncBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
381     fFuncSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
382     
383     fFuncSigBack->SetParameters(1,3.1,.05,1,2.5,1);
384     fFuncSigBack->SetParLimits(0,0,10000000);
385     fFuncSigBack->SetParLimits(1,3.05,3.15);
386     fFuncSigBack->SetParLimits(2,.02,.1);
387   }
388   else if (type==2){
389     // half gaussian, half exponential signal function
390     // exponential background
391     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);
392     fFuncBackground = new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
393     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);
394     fFuncSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
395     
396     fFuncSigBack->SetParLimits(0,0,10000000);
397     fFuncSigBack->SetParLimits(1,3.05,3.15);
398     fFuncSigBack->SetParLimits(2,.02,.1);
399     fFuncSigBack->FixParameter(6,2.5);
400     fFuncSigBack->FixParameter(7,0);
401   }
402 }
403
404
405 //______________________________________________
406 void AliDielectronSignalFunc::Draw(const Option_t* option)
407 {
408   //
409   // Draw the fitted function
410   //
411   
412   TString drawOpt(option);
413   drawOpt.ToLower();
414   
415   Bool_t optStat=drawOpt.Contains("stat");
416   
417   fFuncSigBack->SetNpx(200);
418   fFuncSigBack->SetRange(fIntMin,fIntMax);
419   fFuncBackground->SetNpx(200);
420   fFuncBackground->SetRange(fIntMin,fIntMax);
421   
422   TGraph *grSig=new TGraph(fFuncSigBack);
423   grSig->SetFillColor(kGreen);
424   grSig->SetFillStyle(3001);
425   
426   TGraph *grBack=new TGraph(fFuncBackground);
427   grBack->SetFillColor(kRed);
428   grBack->SetFillStyle(3001);
429   
430   grSig->SetPoint(0,grBack->GetX()[0],grBack->GetY()[0]);
431   grSig->SetPoint(grSig->GetN()-1,grBack->GetX()[grBack->GetN()-1],grBack->GetY()[grBack->GetN()-1]);
432   
433   grBack->SetPoint(0,grBack->GetX()[0],0.);
434   grBack->SetPoint(grBack->GetN()-1,grBack->GetX()[grBack->GetN()-1],0.);
435   
436   fFuncSigBack->SetRange(fFitMin,fFitMax);
437   fFuncBackground->SetRange(fFitMin,fFitMax);
438   
439   if (!drawOpt.Contains("same")){
440     if (fHistDataPM){
441       fHistDataPM->Draw();
442       grSig->Draw("f");
443     } else {
444       grSig->Draw("af");
445     }
446   } else {
447     grSig->Draw("f");
448   }
449   if(fMethod==kFitted) grBack->Draw("f");
450   fFuncSigBack->Draw("same");
451   fFuncSigBack->SetLineWidth(2);
452   if(fMethod==kLikeSign) {
453     fHistDataPP->SetLineWidth(2);
454     fHistDataPP->SetLineColor(6);
455     fHistDataPP->Draw("same");
456     fHistDataMM->SetLineWidth(2);
457     fHistDataMM->SetLineColor(8);
458     fHistDataMM->Draw("same");
459   }
460   
461   if(fMethod==kFitted)
462     fFuncBackground->Draw("same");
463   
464   if (optStat) DrawStats();
465   
466 }