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