]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGDQ/dielectron/AliDielectronSignalBase.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronSignalBase.cxx
... / ...
CommitLineData
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/*
21Base class for signal extraction from a histogram or an array of histograms
22The histogram is assumed to be an inv. mass spectrum,
23the array of histograms is assumed to be an array with inv. mass histograms
24resulting 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
40ClassImp(AliDielectronSignalBase)
41
42TH1F* AliDielectronSignalBase::fgHistSimPM=0x0;
43TObject* AliDielectronSignalBase::fgPeakShape=0x0;
44const char* AliDielectronSignalBase::fgkValueNames[6] = {
45 "Signal","Background","Significance","Signal/Background","Mass","MassWidth"};
46
47AliDielectronSignalBase::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//______________________________________________
80AliDielectronSignalBase::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//______________________________________________
113AliDielectronSignalBase::~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//______________________________________________
128TPaveText* 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//______________________________________________
159void 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//______________________________________________
175Double_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//______________________________________________
201Double_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
232TObject* 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}