Add fast merging option (Diego)
[u/mrichter/AliRoot.git] / PWG3 / dielectron / AliDielectronSignalExt.cxx
CommitLineData
572b0139 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// //
18// Dielectron SignalExt //
19// //
20/*
21Dielectron signal extraction class using functions as input.
22
23A function to describe the signal as well as one to describe the background
24has to be deployed by the user. Alternatively on of the default implementaions
25can be used.
26
27*/
28// //
29///////////////////////////////////////////////////////////////////////////
572b0139 30#include <TF1.h>
31#include <TH1.h>
bc75eeb5 32#include <TH2F.h>
33#include <TLatex.h>
34#include <TLegend.h>
572b0139 35#include <TCanvas.h>
36#include <TMath.h>
37#include <TString.h>
bc75eeb5 38#include <TLine.h>
572b0139 39
40#include <AliLog.h>
41
42#include "AliDielectronSignalExt.h"
43
44ClassImp(AliDielectronSignalExt)
45
46AliDielectronSignalExt::AliDielectronSignalExt() :
bc75eeb5 47 AliDielectronSignalBase()
572b0139 48{
49 //
50 // Default Constructor
51 //
572b0139 52}
53
54//______________________________________________
55AliDielectronSignalExt::AliDielectronSignalExt(const char* name, const char* title) :
bc75eeb5 56 AliDielectronSignalBase(name, title)
572b0139 57{
58 //
59 // Named Constructor
60 //
61}
62
63//______________________________________________
64AliDielectronSignalExt::~AliDielectronSignalExt()
65{
66 //
67 // Default Destructor
68 //
572b0139 69}
70
71//______________________________________________
72void AliDielectronSignalExt::Process(TObjArray* const arrhist)
73{
74 //
75 // signal subtraction. support like-sign subtraction and event mixing method
76 //
77 switch ( fMethod ){
bc75eeb5 78 case kLikeSign :
45b2b1b8 79 case kLikeSignArithm :
572b0139 80 ProcessLS(arrhist); // process like-sign subtraction method
81 break;
82
bc75eeb5 83 case kEventMixing :
572b0139 84 ProcessEM(arrhist); // process event mixing method
85 break;
86
2a14a7b1 87 case kRotation:
88 ProcessRotation(arrhist);
89 break;
90
572b0139 91 default :
92 AliWarning("Subtraction method not supported. Please check SetMethod() function.");
93 }
94}
95
96//______________________________________________
97void AliDielectronSignalExt::ProcessLS(TObjArray* const arrhist)
98{
99 //
100 // signal subtraction
101 //
2a14a7b1 102 fHistDataPP = (TH1*)(arrhist->At(0))->Clone("histPP"); // ++ SE
103 fHistDataPM = (TH1*)(arrhist->At(1))->Clone("histPM"); // +- SE
104 fHistDataMM = (TH1*)(arrhist->At(2))->Clone("histMM"); // -- SE
bc75eeb5 105 fHistDataPP->Sumw2();
106 fHistDataPM->Sumw2();
107 fHistDataMM->Sumw2();
554e40f8 108 fHistDataPP->SetDirectory(0);
109 fHistDataPM->SetDirectory(0);
110 fHistDataMM->SetDirectory(0);
572b0139 111
bc75eeb5 112 // rebin the histograms
113 if (fRebin>1) {
114 fHistDataPP->Rebin(fRebin);
115 fHistDataPM->Rebin(fRebin);
116 fHistDataMM->Rebin(fRebin);
117 }
118
2a14a7b1 119 fHistSignal = new TH1D("HistSignal", "Like-Sign substracted signal",
bc75eeb5 120 fHistDataPM->GetXaxis()->GetNbins(),
121 fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
554e40f8 122 fHistSignal->SetDirectory(0);
2a14a7b1 123 fHistBackground = new TH1D("HistBackground", "Like-sign contribution",
bc75eeb5 124 fHistDataPM->GetXaxis()->GetNbins(),
125 fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
554e40f8 126 fHistBackground->SetDirectory(0);
127
572b0139 128 // fill out background and subtracted histogram
bc75eeb5 129 for(Int_t ibin=1; ibin<=fHistDataPM->GetXaxis()->GetNbins(); ibin++) {
130 Float_t pm = fHistDataPM->GetBinContent(ibin);
131 Float_t pp = fHistDataPP->GetBinContent(ibin);
132 Float_t mm = fHistDataMM->GetBinContent(ibin);
133 Float_t epm = fHistDataPM->GetBinError(ibin);
134
135 Float_t background = 2*TMath::Sqrt(pp*mm);
136 Float_t ebackground = TMath::Sqrt(mm+pp);
45b2b1b8 137 if (fMethod==kLikeSignArithm){
138 //Arithmetic mean instead of geometric
139 background=(pp+mm);
140 ebackground=TMath::Sqrt(pp+mm);
141 if (TMath::Abs(ebackground)<1e-30) ebackground=1;
142 }
554e40f8 143// Float_t signal = pm - background;
144// Float_t error = TMath::Sqrt(epm*epm+mm+pp);
bc75eeb5 145
554e40f8 146 fHistSignal->SetBinContent(ibin, pm);
147 fHistSignal->SetBinError(ibin, epm);
bc75eeb5 148 fHistBackground->SetBinContent(ibin, background);
149 fHistBackground->SetBinError(ibin, ebackground);
572b0139 150 }
fb7d2d99 151 //scale histograms to match integral between fScaleMin and fScaleMax
554e40f8 152 if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
153
154 //subract background
155 fHistSignal->Add(fHistBackground,-1);
156
bc75eeb5 157 // signal
158 fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
159 fHistSignal->FindBin(fIntMax), fErrors(0));
160 // background
161 fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
48609e3d 162 fHistBackground->FindBin(fIntMax),
163 fErrors(1));
bc75eeb5 164 // S/B and significance
572b0139 165 SetSignificanceAndSOB();
166
bc75eeb5 167 fProcessed = kTRUE;
572b0139 168}
169
170//______________________________________________
171void AliDielectronSignalExt::ProcessEM(TObjArray* const arrhist)
172{
173 //
174 // event mixing. not implemented yet.
175 //
176 printf("event mixing method is not implemented yet. Like-sign method will be used.\n");
177 ProcessLS(arrhist);
178}
179
180//______________________________________________
2a14a7b1 181void AliDielectronSignalExt::ProcessRotation(TObjArray* const arrhist)
182{
183 //
184 // signal subtraction
185 //
186 fHistDataPM = (TH1*)(arrhist->At(1))->Clone("histPM"); // +- SE
187 if (!fHistDataPM){
188 AliError("Unlike sign histogram not available. Cannot extract the signal.");
189 return;
190 }
191 fHistDataPM->Sumw2();
192
193 fHistBackground=(TH1*)arrhist->At(10)->Clone("histRotation");
194 if (!fHistBackground){
195 AliError("Histgram from rotation not available. Cannot extract the signal.");
196 delete fHistDataPM;
197 fHistDataPM=0x0;
198 return;
199 }
ba15fdfb 200 fHistBackground->Sumw2();
201
202 // rebin the histograms
203 if (fRebin>1) {
204 fHistDataPM->Rebin(fRebin);
205 fHistBackground->Rebin(fRebin);
206 }
2a14a7b1 207
554e40f8 208 //scale histograms to match integral between fScaleMin and fScaleMax
209 if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax);
210
2a14a7b1 211 fHistSignal=(TH1*)fHistDataPM->Clone("histSignal");
212 fHistSignal->Add(fHistBackground,-1.);
213
214 // signal
215 fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
216 fHistSignal->FindBin(fIntMax), fErrors(0));
217 // background
218 fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
219 fHistBackground->FindBin(fIntMax),
220 fErrors(1));
221 // S/B and significance
222 SetSignificanceAndSOB();
223
224 fProcessed = kTRUE;
225
226}
227
228//______________________________________________
572b0139 229void AliDielectronSignalExt::Draw(const Option_t* option)
230{
231 //
232 // Draw the fitted function
233 //
234 TString drawOpt(option);
235 drawOpt.ToLower();
236
bc75eeb5 237 Float_t minY = 0.001;
238 Float_t maxY = 1.2*fHistDataPM->GetMaximum();
239 Float_t minX = 1.001*fHistDataPM->GetXaxis()->GetXmin();
240 Float_t maxX = 0.999*fHistDataPM->GetXaxis()->GetXmax();
241 Int_t binSize = Int_t(1000*fHistDataPM->GetBinWidth(1)); // in MeV
242 Float_t minMinY = fHistSignal->GetMinimum();
243
48609e3d 244 TCanvas *cSub = new TCanvas(Form("%s", fName.Data()),Form("%s", fTitle.Data()),1400,1000);
bc75eeb5 245 cSub->SetLeftMargin(0.15);
246 cSub->SetRightMargin(0.0);
247 cSub->SetTopMargin(0.002);
248 cSub->SetBottomMargin(0.0);
249 cSub->Divide(2,2,0.,0.);
572b0139 250 cSub->Draw();
251
bc75eeb5 252 TVirtualPad* pad = cSub->cd(1);
253 pad->SetLeftMargin(0.15);
254 pad->SetRightMargin(0.0);
255 pad->SetTopMargin(0.005);
256 pad->SetBottomMargin(0.0);
257 TH2F *range1=new TH2F("range1","",10,minX,maxX,10,minY,maxY);
258 range1->SetStats(kFALSE);
259 range1->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
260 range1->GetYaxis()->CenterTitle();
261 range1->GetYaxis()->SetLabelSize(0.05);
262 range1->GetYaxis()->SetTitleSize(0.06);
263 range1->GetYaxis()->SetTitleOffset(0.8);
264 range1->Draw();
265 fHistDataPM->SetLineColor(1);
266 fHistDataPM->SetLineWidth(2);
267 // fHistDataPM->SetMarkerStyle(21);
268 fHistDataPM->Draw("Psame");
269 TLatex *latex = new TLatex();
270 latex->SetNDC();
271 latex->SetTextSize(0.05);
272 latex->DrawLatex(0.2, 0.95, "Background un-substracted");
273 TLine line;
274 line.SetLineWidth(1);
275 line.SetLineStyle(2);
276 line.DrawLine(fIntMin, minY, fIntMin, maxY);
277 line.DrawLine(fIntMax, minY, fIntMax, maxY);
278
279 pad = cSub->cd(2);
280 pad->SetLeftMargin(0.);
281 pad->SetRightMargin(0.005);
282 pad->SetTopMargin(0.005);
283 pad->SetBottomMargin(0.0);
284 TH2F *range2=new TH2F("range2","",10,minX,maxX,10,minY,maxY);
285 range2->SetStats(kFALSE);
286 range2->Draw();
287 fHistBackground->SetLineColor(4);
288 fHistBackground->SetLineWidth(2);
289 // fHistBackground->SetMarkerColor(4);
290 // fHistBackground->SetMarkerStyle(6);
291 fHistBackground->Draw("Psame");
292 latex->DrawLatex(0.05, 0.95, "Like-sign background");
293 line.DrawLine(fIntMin, minY, fIntMin, maxY);
294 line.DrawLine(fIntMax, minY, fIntMax, maxY);
295 TLegend *legend = new TLegend(0.65, 0.70, 0.98, 0.98);
296 legend->SetFillColor(0);
297 legend->SetMargin(0.15);
298 legend->AddEntry(fHistDataPM, "N_{+-}", "l");
299 legend->AddEntry(fHistDataPP, "N_{++}", "l");
300 legend->AddEntry(fHistDataMM, "N_{--}", "l");
301 legend->AddEntry(fHistSignal, "N_{+-} - 2 #sqrt{N_{++} #times N_{--}}", "l");
302 legend->AddEntry(fHistBackground, "2 #sqrt{N_{++} #times N_{--}}", "l");
303 legend->Draw();
304
572b0139 305
bc75eeb5 306 pad = cSub->cd(3);
307 pad->SetLeftMargin(0.15);
308 pad->SetRightMargin(0.0);
309 pad->SetTopMargin(0.0);
310 pad->SetBottomMargin(0.15);
311 TH2F *range3=new TH2F("range3","",10,minX,maxX,10,minMinY,maxY);
312 range3->SetStats(kFALSE);
313 range3->GetYaxis()->SetTitle(Form("entries [counts per %d MeV bin]", binSize));
314 range3->GetYaxis()->CenterTitle();
315 range3->GetYaxis()->SetLabelSize(0.05);
316 range3->GetYaxis()->SetTitleSize(0.06);
317 range3->GetYaxis()->SetTitleOffset(0.8);
318 range3->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
319 range3->GetXaxis()->CenterTitle();
320 range3->GetXaxis()->SetLabelSize(0.05);
321 range3->GetXaxis()->SetTitleSize(0.06);
322 range3->GetXaxis()->SetTitleOffset(1.0);
323 range3->Draw();
324 fHistDataPM->Draw("Psame");
325 fHistDataPP->SetLineWidth(2);
326 fHistDataPP->SetLineColor(6);
327 fHistDataMM->SetLineWidth(2);
328 fHistDataMM->SetLineColor(8);
329 fHistDataPP->Draw("Psame");
330 fHistDataMM->Draw("Psame");
331 line.DrawLine(minX, 0.,maxX, 0.);
332 line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
333 line.DrawLine(fIntMax, minMinY, fIntMax, maxY);
334
335 pad = cSub->cd(4);
336 pad->SetLeftMargin(0.0);
337 pad->SetRightMargin(0.005);
338 pad->SetTopMargin(0.0);
339 pad->SetBottomMargin(0.15);
340 TH2F *range4=new TH2F("range4","",10,minX,maxX,10,minMinY,maxY);
341 range4->SetStats(kFALSE);
342 range4->GetXaxis()->SetTitle("inv. mass [GeV/c^{2}]");
343 range4->GetXaxis()->CenterTitle();
344 range4->GetXaxis()->SetLabelSize(0.05);
345 range4->GetXaxis()->SetTitleSize(0.06);
346 range4->GetXaxis()->SetTitleOffset(1.0);
347 range4->Draw();
348 fHistSignal->SetLineWidth(2);
349 fHistSignal->SetLineColor(2);
350 fHistSignal->Draw("Psame");
351 latex->DrawLatex(0.05, 0.95, "Like-sign background substracted");
352 if(fProcessed) DrawStats(0.05, 0.6, 0.5, 0.9);
353 line.DrawLine(minX, 0.,maxX, 0.);
354 line.DrawLine(fIntMin, minMinY, fIntMin, maxY);
355 line.DrawLine(fIntMax, minMinY, fIntMax, maxY);
48609e3d 356
357 cSub->SaveAs(Form("%s_summary.png", fName.Data()));
572b0139 358}
359