CbmRoot
Loading...
Searching...
No Matches
PairAnalysisSignalExt.cxx
Go to the documentation of this file.
1/* Copyright (C) 2015-2020 Justus-Liebig-Universitaet Giessen, Giessen
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Julian Book [committer] */
4
6// Dielectron SignalExt //
7// //
8// //
9/*
10Ext class for signal extraction from a histogram or an array of histograms
11The histogram is assumed to be an inv. mass spectrum,
12the array of histograms is assumed to be an array with inv. mass histograms
13resulting from single and mixed events, as defined in PairAnalysis.cxx
14
15*/
16// //
18
19
20#include <TCollection.h>
21#include <TDatabasePDG.h>
22#include <TF1.h>
23#include <TH1.h>
24#include <TList.h>
25#include <TPaveText.h>
26#include <TProfile.h>
27#include <TProfile2D.h>
28#include <TROOT.h>
29#include <TVectorT.h>
30//#include <TH1.h>
31#include <TCanvas.h>
32#include <TGaxis.h>
33#include <TH2F.h>
34#include <TLatex.h>
35#include <TLegend.h>
36#include <TLine.h>
37#include <TMath.h>
38#include <TString.h>
39
40#include "PairAnalysis.h"
41#include "PairAnalysisFunction.h"
42#include "PairAnalysisHelper.h"
43#include "PairAnalysisSignalExt.h"
44#include "PairAnalysisSignalMC.h"
45#include "PairAnalysisStyler.h"
46
47ClassImp(PairAnalysisSignalExt)
48
49 TH1F* PairAnalysisSignalExt::fgHistSimPM = 0x0;
50TObject* PairAnalysisSignalExt::fgPeakShape = 0x0;
51const char* PairAnalysisSignalExt::fgkValueNames[7] = {"S", "B", "S/#sqrt{S+B}", "S/B",
52 "Mass", "MassWidth", "ChiSqNDFmatch"};
53const char* PairAnalysisSignalExt::fgkBackgroundMethodNames[11] = {"FittedMC",
54 "Fitted",
55 "like-sign",
56 "like-sign (arithm.)",
57 "like-sign #times R(#Deltam)",
58 "like-sign (arithm.) #times R(#Deltam)",
59 "fitted like-sign",
60 "mixed event",
61 "fitted mixed event",
62 "track rotation",
63 "cocktail"};
64
65PairAnalysisSignalExt::PairAnalysisSignalExt() : PairAnalysisSignalExt("PairAnalysisSignalExt", "title")
66{
67 //
68 // Default Constructor
69 //
70}
71
72//______________________________________________
73PairAnalysisSignalExt::PairAnalysisSignalExt(const char* name, const char* title)
74 : PairAnalysisFunction(name, title)
75 , fValues(7)
76 , fErrors(7)
77{
78 //
79 // Named Constructor
80 //
81}
82
83//______________________________________________
84PairAnalysisSignalExt::PairAnalysisSignalExt(const PairAnalysisSignalExt& c)
85 : PairAnalysisFunction(c.GetName(), c.GetTitle())
86 , fArrHists(c.fArrHists)
87 , fArrCocktail(c.fArrCocktail)
88 , fHistSignal(c.GetSignalHistogram())
89 , fHistSB(c.GetSoverBHistogram())
90 , fHistSgn(c.GetSignificanceHistogram())
91 , fHistBackground(c.GetBackgroundHistogram())
92 , fHistCocktail(c.GetCocktailHistogram())
93 , fHistDataPM(c.GetUnlikeSignHistogram())
94 , fHistRfactor(c.GetRfactorHistogram())
95 , fHistSignalMC(c.GetMCSignalShape())
96 , fValues(c.GetValues())
97 , fErrors(c.GetErrors())
98 , fIntMin(c.GetIntegralMin())
99 , fIntMax(c.GetIntegralMax())
100 , fRebin(c.GetRebin())
101 ,
102 // fRebinStat(1.), //TODO
103 // fBinLimits(0x0), //TODO
104 fMethod(c.GetMethod())
105 , fScaleMin(c.GetScaleMin())
106 , fScaleMax(c.GetScaleMax())
107 , fScaleMin2(c.GetScaleMin2())
108 , fScaleMax2(c.GetScaleMax2())
109 , fScaleFactor(c.GetScaleFactor())
110 , fPeakMethod(c.GetExtractionMethod())
111// fExtrFunc(0x0) //TODO: needed
112{
113 //
114 // Copy Constructor
115 //
116}
117
118//______________________________________________
119PairAnalysisSignalExt::~PairAnalysisSignalExt()
120{
121 //
122 // Default Destructor
123 //
124 // TODO: add new datamembers
125 if (fArrHists) delete fArrHists;
126 if (fArrCocktail) delete fArrCocktail;
127 if (fHistSignal) delete fHistSignal;
128 if (fHistSB) delete fHistSB;
129 if (fHistSgn) delete fHistSgn;
130 if (fHistBackground) delete fHistBackground;
131 if (fHistCocktail) delete fHistCocktail;
132 if (fHistDataPP) delete fHistDataPP;
133 if (fHistDataPM) delete fHistDataPM;
134 if (fHistDataMM) delete fHistDataMM;
135 if (fHistDataME) delete fHistDataME;
136 if (fHistRfactor) delete fHistRfactor;
137 if (fHistSignalMC) delete fHistSignalMC;
138 if (fExtrFunc) delete fExtrFunc;
139 if (fBinLimits) delete fBinLimits;
140}
141
142//______________________________________________
143TPaveText* PairAnalysisSignalExt::DrawStats(Double_t x1 /*=0.*/, Double_t y1 /*=0.*/, Double_t x2 /*=0.*/,
144 Double_t y2 /*=0.*/, TString opt /*="pRnbsSmrc"*/)
145{
146 //
147 // Draw extracted values in a TPaveText
148 //
149 // with the corners x1,y2,x2,y2
150 // use option string to select information displayed:
151 // p := Particle name using fPOIpdg
152 // R := integration Range
153 // n := Number of signal counts
154 // b := Background
155 // s := signal-to-bgrd
156 // S := Significance
157 // m := mass position
158 // r := mass resolution (sigma or fwhm)
159 // c := matching Chi2/ndf
160 //
161 if (TMath::Abs(x1) < 1e-20 && TMath::Abs(x2) < 1e-20) {
162 x1 = .6;
163 x2 = .9;
164 y1 = .7;
165 y2 = .9;
166 }
167 if (y1 < 0.5) y2 = y1 + 0.025 * opt.Length();
168 else
169 y1 = y2 - 0.025 * opt.Length();
170
171 // particle of interest
172 // TParticlePDG *fPOI = TDatabasePDG::Instance()->GetParticle(fPOIpdg);
173
174 TPaveText* t = new TPaveText(x1, y1, x2, y2, "brNDC");
175 t->SetFillColor(0);
176 t->SetFillStyle(kFEmpty);
177 t->SetBorderSize(0);
178 t->SetTextAlign(12);
179 if (opt.Contains("p")) t->AddText(Form("#bf{%s}", (fPOI ? fPOI->GetName() : "particle not found")));
180 if (opt.Contains("R")) t->AddText(Form("%.2f < %s < %.2f GeV/c^{2}", fIntMin, "m_{inv}", fIntMax));
181 if (opt.Contains("n")) t->AddText(Form("%s: %.4g #pm %.4g", fgkValueNames[0], fValues(0), fErrors(0)));
182 if (opt.Contains("b")) t->AddText(Form("%s: %.4g #pm %.4g", fgkValueNames[1], fValues(1), fErrors(1)));
183 Int_t smallS2B = (fValues(3) < 0.1 ? 1 : 100);
184 if (opt.Contains("s"))
185 t->AddText(Form("%s: %.2f #pm %.2f%s", fgkValueNames[3], fValues(3) * 100 / smallS2B, fErrors(3) * 100 / smallS2B,
186 smallS2B > 1 ? "" : "%"));
187 if (opt.Contains("S")) t->AddText(Form("%s: %.1f #pm %.1f", fgkValueNames[2], fValues(2), fErrors(2)));
188 if (opt.Contains("m") && fValues(4) > 0) t->AddText(Form("Mass: %.2f #pm %.2f GeV/c^{2}", fValues(4), fErrors(4)));
189 if (opt.Contains("r") && fValues(5) > 0)
190 t->AddText(Form("Mass res.: %.1f #pm %.1f MeV/c^{2}", 1000 * fValues(5), 1000 * fErrors(5)));
191 if (opt.Contains("c") && fValues(6) > 0) t->AddText(Form("Match chi2/ndf.: %.1f", fValues(6)));
192 t->Draw();
193
194 return t;
195}
196
197//______________________________________________
198void PairAnalysisSignalExt::Print(Option_t* /*option*/) const
199{
200 //
201 // Print the statistics
202 //
203 printf("Signal : %.5g #pm %.5g (+-%.3g%%)\n", fValues(0), fErrors(0), 100 * fErrors(0) / fValues(0));
204 printf("Backgnd: %.5g #pm %.5g (+-%.3g%%)\n", fValues(1), fErrors(1), 100 * fErrors(1) / fValues(1));
205 printf("Signif : %.5g #pm %.5g (+-%.3g%%)\n", fValues(2), fErrors(2), 100 * fErrors(2) / fValues(2));
206 printf("SoB : %.5g #pm %.5g (+-%.3g%%)\n", fValues(3), fErrors(3), 100 * fErrors(3) / fValues(3));
207 if (fValues(4) > 0) printf("Mass: %.5g #pm %.5g\n", fValues(4), fErrors(4));
208 if (fValues(5) > 0) printf("Mass res.: %.5g #pm %.5g\n", fValues(5), fErrors(5));
209 if (fValues(6) > 0) printf("Match chi2/ndf: %.5g \n", fValues(6));
210 printf("Scale : %.5g \n", fScaleFactor);
211 printf("Mass int.: %.5g - %.5g \n", fIntMin, fIntMax);
212}
213
214//______________________________________________
215Double_t PairAnalysisSignalExt::ScaleHistograms(TH1* histRaw, TH1* histBackground, Double_t intMin, Double_t intMax)
216{
217 //
218 // scale histBackground to match the integral of histRaw in the interval intMin, intMax
219 //
220
221 //protect using over and underflow bins in normalisation calculation
222 if (intMin < histRaw->GetXaxis()->GetXmin()) intMin = histRaw->GetXaxis()->GetXmin();
223 if (intMin < histBackground->GetXaxis()->GetXmin()) intMin = histBackground->GetXaxis()->GetXmin();
224
225 if (intMax > histRaw->GetXaxis()->GetXmax())
226 intMax = histRaw->GetXaxis()->GetXmax() - histRaw->GetBinWidth(histRaw->GetNbinsX()) / 2.;
227 if (intMax > histBackground->GetXaxis()->GetXmax())
228 intMax = histBackground->GetXaxis()->GetXmax() - histBackground->GetBinWidth(histBackground->GetNbinsX()) / 2.;
229
230 Double_t intRaw = histRaw->Integral(histRaw->FindBin(intMin), histRaw->FindBin(intMax));
231 Double_t intBack = histBackground->Integral(histBackground->FindBin(intMin), histBackground->FindBin(intMax));
232 fScaleFactor = intBack > 0 ? intRaw / intBack : 1.;
233 // scale
234 if (histBackground->GetDefaultSumw2()) histBackground->Sumw2();
235 histBackground->Scale(fScaleFactor);
236
237 return fScaleFactor;
238}
239
240//______________________________________________
241Double_t PairAnalysisSignalExt::ScaleHistograms(TH1* histRaw, TH1* histBackground, Double_t intMin, Double_t intMax,
242 Double_t intMin2, Double_t intMax2)
243{
244 //
245 // scale histBackground to match the integral of "histRaw" in the interval "intMin", "intMax" and "intMin2", "intMax2"
246 //
247
248 if (TMath::Abs(intMin2 - intMax2) < 0.00001) return (ScaleHistograms(histRaw, histBackground, intMin, intMax));
249
250 //protect using over and underflow bins in normalisation calculation
251 if (intMin < histRaw->GetXaxis()->GetXmin()) intMin = histRaw->GetXaxis()->GetXmin();
252 if (intMin < histBackground->GetXaxis()->GetXmin()) intMin = histBackground->GetXaxis()->GetXmin();
253
254 if (intMax2 > histRaw->GetXaxis()->GetXmax())
255 intMax2 = histRaw->GetXaxis()->GetXmax() - histRaw->GetBinWidth(histRaw->GetNbinsX()) / 2.;
256 if (intMax2 > histBackground->GetXaxis()->GetXmax())
257 intMax2 = histBackground->GetXaxis()->GetXmax() - histBackground->GetBinWidth(histBackground->GetNbinsX()) / 2.;
258
259 Double_t intRaw = histRaw->Integral(histRaw->FindBin(intMin), histRaw->FindBin(intMax));
260 Double_t intBack = histBackground->Integral(histBackground->FindBin(intMin), histBackground->FindBin(intMax));
261 intRaw += histRaw->Integral(histRaw->FindBin(intMin2), histRaw->FindBin(intMax2));
262 intBack += histBackground->Integral(histBackground->FindBin(intMin2), histBackground->FindBin(intMax2));
263
264 fScaleFactor = intBack > 0 ? intRaw / intBack : 1.;
265 // scale
266 if (histBackground->GetDefaultSumw2()) histBackground->Sumw2();
267 histBackground->Scale(fScaleFactor);
268
269 return fScaleFactor;
270}
271
272//______________________________________________
273TH1* PairAnalysisSignalExt::MergeObjects(TH1* obj1, TH1* obj2, Double_t val)
274{
275 //
276 // function to merge all TH1 inherited objects
277 // (needed because TProfile::Add with negative 'val' does not what is needed)
278 // what about TProfile2D (inherits from TH2D) ?
279 //
280 if (!obj1) return obj2;
281 if (!obj2) return obj1;
282
283 TString key = Form("%s_%s", obj1->GetName(), "Signal");
284 if (obj1->IsA() == TProfile2D::Class() && obj2->IsA() == TProfile2D::Class()) {
285 if (val >= 0.) {
286 //summation
287 ((TProfile2D*) obj1)->Add(((TProfile2D*) obj2), val);
288 return (TH1*) obj1->Clone(key.Data()); //TOCHECK: clone needed??
289 }
290 else
291 return 0x0;
292 }
293 else if (obj1->IsA() == TProfile::Class() && obj2->IsA() == TProfile::Class()) {
294 if (val >= 0.) {
295 //summation
296 ((TProfile*) obj1)->Add(((TProfile*) obj2), val);
297 return (TH1*) obj1->Clone(key.Data()); //TOCHECK: clone needed??
298 }
299 else {
300 // ONLY 1D subtraction
301 TH1D* histSign = new TH1D(key.Data(), "", obj1->GetXaxis()->GetNbins(), obj1->GetXaxis()->GetXmin(),
302 obj1->GetXaxis()->GetXmax());
303 if (histSign->GetDefaultSumw2()) histSign->Sumw2();
304 histSign->SetDirectory(0);
305
306 for (Int_t i = 0; i <= obj1->GetNbinsX(); i++) {
307 histSign->SetBinContent(i, obj1->GetBinContent(i) - obj2->GetBinContent(i));
308 histSign->SetBinError(
309 i, TMath::Sqrt(obj1->GetBinError(i) * obj1->GetBinError(i) - obj2->GetBinError(i) * obj2->GetBinError(i)));
310 }
311 return histSign;
312 }
313
314 /*
315 TList listH;
316 TString listHargs;
317 listHargs.Form("((TCollection*)0x%lx)", (ULong_t)&listH);
318 Int_t error = 0;
319 listH.Add(obj2);
320 obj1->Execute("Merge", listHargs.Data(), &error);
321 */
322 }
323 else { // Histograms
324 if (val < 0.) { // subtraction
325 TH1* histSign = (TH1*) obj1->Clone(key.Data());
326 histSign->Add(obj2, val);
327 return histSign;
328 }
329 else { // merge
330 obj1->Add(obj2, val);
331 return obj1;
332 }
333 }
334}
335
336//______________________________________________
337TObject* PairAnalysisSignalExt::DescribePeakShape(ESignalExtractionMethod method, Bool_t replaceValErr, TH1F* mcShape)
338{
339 //
340 // Describe the extracted peak by the selected "method" and overwrite signal etc if activated
341 //
342
343 if (replaceValErr) fPeakMethod = method;
344 Double_t data = 0.;
345 Double_t mc = 0.;
346 Double_t massPOI = TDatabasePDG::Instance()->GetParticle(fPOIpdg)->Mass();
347 Double_t sigPOI = massPOI * 0.02;
348 Double_t nPOI = fHistSignal->GetBinContent(fHistSignal->FindBin(massPOI));
349 Double_t binWidth = fHistSignal->GetBinWidth(fHistSignal->FindBin(massPOI));
350 TF1* fit = 0x0;
351 Int_t fitResult = 0;
352 Int_t parMass = -1;
353 Int_t parSigma = -1;
354 if (!fExtrFunc) fExtrFunc = new PairAnalysisFunction();
355 //PairAnalysisSignalFit *fExtrFunc = new PairAnalysisSignalFit();
356 // PairAnalysisSignalFunc fct;// = 0;//new PairAnalysisSignalFunc();
357
358 Info("DescribePeakShape", "Signal extraction method: %d", (Int_t) fPeakMethod);
359
360 // do the scaling/fitting
361 switch (method) {
362 case kBinCounting: /*nothing needs to be done*/ break;
363
364 case kMCScaledMax:
365 if (!mcShape) {
366 Error("DescribePeakShape", "No MC histogram passed. Returning.");
367 return 0x0;
368 }
369 data = fHistSignal->GetBinContent(fHistSignal->FindBin(massPOI));
370 mc = mcShape->GetBinContent(fHistSignal->FindBin(massPOI));
371 mcShape->Scale(data / mc);
372 break;
373
374 case kMCScaledInt:
375 if (!mcShape) {
376 Error("DescribePeakShape", "No MC histogram passed. Returning.");
377 return 0x0;
378 }
379 if (mcShape->GetBinWidth(1) != fHistSignal->GetBinWidth(1))
380 printf(" WARNING: MC and signal histogram have different bin widths. \n");
381 data = fHistSignal->Integral(fHistSignal->FindBin(fIntMin), fHistSignal->FindBin(fIntMax));
382 mc = mcShape->Integral(mcShape->FindBin(fIntMin), mcShape->FindBin(fIntMax));
383 mcShape->Scale(data / mc);
384 break;
385
386 case kMCFitted:
387 if (!mcShape && !fgHistSimPM) {
388 printf(" ERROR: No MC histogram passed or set. Returning. \n");
389 return 0x0;
390 }
391 if (!fgHistSimPM) fgHistSimPM = mcShape;
392 if (fgHistSimPM->GetBinWidth(1) != fHistSignal->GetBinWidth(1)) fgHistSimPM->Rebin(fRebin);
393 // fit = new TF1("fitMC",PairAnalysisSignalFunc::PeakFunMC,fFitMin,fFitMax,1);
394 fit = new TF1("MC", fExtrFunc, &PairAnalysisFunction::PeakFunMC, fFitMin, fFitMax, 1);
395 fit->SetParNames("N");
396 fitResult = fHistSignal->Fit(fit, "RNI0Q");
397 break;
398
399 case kCrystalBall:
400 fit = new TF1("Crystal Ball", fExtrFunc, &PairAnalysisFunction::PeakFunCB, fFitMin, fFitMax, 5);
401 fit->SetParNames("alpha", "n", "meanx", "sigma", "N");
402 // fit->SetParameters(-.2,5.,gMjpsi,.06,20);
403 // fit->SetParameters(1.,3.6,gMjpsi,.08,700);
404 fit->SetParameters(0.4, 4.0, massPOI, sigPOI, 1.3 * nPOI);
405 fit->SetParLimits(0, 0.0, 1.);
406 fit->SetParLimits(1, 0.01, 10.);
407 fit->SetParLimits(2, massPOI - sigPOI, massPOI + sigPOI);
408 fit->SetParLimits(3, sigPOI / 5, 5 * sigPOI);
409 fit->SetParLimits(4, 0.2 * nPOI, 2.0 * nPOI);
410 parMass = 2;
411 parSigma = 3;
412 fitResult = fHistSignal->Fit(fit, "RNI0Q");
413 break;
414
415 case kGaus:
416 fit = new TF1("Gaussisan", fExtrFunc, &PairAnalysisFunction::PeakFunGauss, fFitMin, fFitMax, 3);
417 //fit = new TF1("fitGaus","gaus",fFitMin,fFitMax);
418 fit->SetParNames("N", "meanx", "sigma");
419 fit->SetParameters(1.3 * nPOI, massPOI, sigPOI);
420 fit->SetParLimits(0, 0.2 * nPOI, 2.0 * nPOI);
421 fit->SetParLimits(1, massPOI - sigPOI, massPOI + sigPOI);
422 fit->SetParLimits(2, sigPOI / 5, 5 * sigPOI);
423 parMass = 1;
424 parSigma = 2;
425 // fit->Print("V");
426 fitResult = fHistSignal->Fit(fit, "RNI0");
427 break;
428
429 case kUserFunc:
430 if (!fExtrFunc || !fExtrFunc->GetCombinedFunction()) {
431 Fatal("DescribePeakShape", "Function class not added!");
432 return 0x0;
433 }
434 fit = fExtrFunc->GetCombinedFunction();
435 fitResult = fHistSignal->Fit(fit, "RNI0Q");
436 break;
437 }
438
439 // warning in case of fit issues
440 if (fitResult != 0) Warning("DescripePeakShape", "fit has error/issue (%d)", fitResult);
441
442 // store chi2/ndf of the fit
443 if (fit) fValues(6) = (fit->GetNDF() ? fit->GetChisquare() / fit->GetNDF() : -1.);
444
445 // overwrite values and errors if requested
446 if (replaceValErr) {
447 switch (method) {
448 case kBinCounting:
449 fValues(0) =
450 fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin), fHistSignal->FindBin(fIntMax), fErrors(0));
451 SetSignificanceAndSOB();
452 break;
453 case kMCScaledMax:
454 case kMCScaledInt:
455 fValues(0) = mcShape->IntegralAndError(mcShape->FindBin(fIntMin), mcShape->FindBin(fIntMax), fErrors(0));
456 SetSignificanceAndSOB();
457 break;
458
459 case kMCFitted:
460 case kCrystalBall:
461 case kGaus:
462 case kUserFunc:
463 fValues(0) = fit->Integral(fIntMin, fIntMax) / binWidth;
464 fErrors(0) = fit->IntegralError(fIntMin, fIntMax) / binWidth;
465 SetSignificanceAndSOB();
466 break;
467 }
468
469 // set mass position and width
470 if (parMass >= 0) {
471 fValues(4) = fit->GetParameter(parMass);
472 fErrors(4) = fit->GetParError(parMass);
473 }
474 if (parSigma >= 0) {
475 fValues(5) = fit->GetParameter(parSigma);
476 fErrors(5) = fit->GetParError(parSigma);
477 }
478 else {
479 // calculate FWHM
480 SetFWHM();
481 }
482 }
483
484 // set the peak method obj
485 switch (method) {
486 case kBinCounting:
487 if (replaceValErr) fgPeakShape = (TH1F*) fHistSignal->Clone("BinCount");
488 // delete fExtrFunc;
489 return (TH1F*) fHistSignal->Clone("BinCountReturn");
490 break;
491 case kMCScaledMax:
492 case kMCScaledInt:
493 if (replaceValErr) fgPeakShape = mcShape;
494 // delete fExtrFunc;
495 return mcShape;
496 break;
497 case kMCFitted:
498 case kCrystalBall:
499 case kGaus:
500 if (fgHistSimPM) fit->SetName(Form("mcShapeFunc-%s", fgHistSimPM->GetName()));
501 if (replaceValErr) fgPeakShape = fit;
502 // delete fExtrFunc;
503 fPeakIsTF1 = kTRUE;
504 return fit;
505 case kUserFunc:
506 fit->SetName(Form("UserFunc"));
507 if (replaceValErr) fgPeakShape = fit;
508 fit->Print();
509 fPeakIsTF1 = kTRUE;
510 break;
511 }
512
513 // printf("true integration range: %f %f \n",
514 // fHistSignal->GetBinLowEdge(fHistSignal->FindBin(fIntMin)),
515 // fHistSignal->GetBinLowEdge(fHistSignal->FindBin(fIntMax))+fHistSignal->GetBinWidth(fHistSignal->FindBin(fIntMax)));
516 //delete fExtrFunc;
517 // if(fgPeakShape->IsA()==TF1::Class()) fPeakIsTF1=kTRUE;
518 if (replaceValErr && fgPeakShape->IsA() == TF1::Class()) fPeakIsTF1 = kTRUE;
519 return fgPeakShape;
520}
521
522//______________________________________________
523void PairAnalysisSignalExt::Process(TObjArray* const arrhist)
524{
525 //
526 // process signal subtraction
527 //
528 fArrHists = arrhist;
529 fArrHists->SetOwner(kFALSE);
530
531 // calculate optimal binning if configured
532 if (fRebinStat < 1. && fBinLimits == 0x0) {
533 fBinLimits =
534 PairAnalysisHelper::MakeStatBinLimits((TH1*) FindObject(arrhist, PairAnalysis::EPairType::kSEPM), fRebinStat);
535 }
536
537
538 // clean up spectra
539 if (fProcessed && 0) { //TODO: not needed??
540 if (fHistDataPP) {
541 delete fHistDataPP;
542 fHistDataPP = 0x0;
543 }
544 if (fHistDataPM) {
545 delete fHistDataPM;
546 fHistDataPM = 0x0;
547 }
548 if (fHistDataMM) {
549 delete fHistDataMM;
550 fHistDataMM = 0x0;
551 }
552
553 if (fHistDataME) {
554 delete fHistDataME;
555 fHistDataME = 0x0;
556 }
557 if (fHistRfactor) {
558 delete fHistRfactor;
559 fHistRfactor = 0x0;
560 }
561
562 if (fHistMixPP) {
563 delete fHistMixPP;
564 fHistMixPP = 0x0;
565 }
566 if (fHistMixPM) {
567 delete fHistMixPM;
568 fHistMixPM = 0x0;
569 }
570 if (fHistMixMM) {
571 delete fHistMixMM;
572 fHistMixMM = 0x0;
573 }
574 if (fHistMixMP) {
575 delete fHistMixMP;
576 fHistMixMP = 0x0;
577 }
578
579 if (fHistDataTR) {
580 delete fHistDataTR;
581 fHistDataTR = 0x0;
582 }
583
584 if (fHistSignal) {
585 delete fHistSignal;
586 fHistSignal = 0x0;
587 }
588 if (fHistBackground) {
589 delete fHistBackground;
590 fHistBackground = 0x0;
591 }
592 if (fHistSB) {
593 delete fHistSB;
594 fHistSB = 0x0;
595 }
596 if (fHistSgn) {
597 delete fHistSgn;
598 fHistSgn = 0x0;
599 }
600 if (fHistCocktail) {
601 delete fHistCocktail;
602 fHistCocktail = 0x0;
603 }
604 }
605
607 // SE ++
608 fHistDataPP = (TH1*) FindObject(arrhist, PairAnalysis::EPairType::kSEPP);
609 if (fHistDataPP) {
610 if (fBinLimits) {
611 fHistDataPP = fHistDataPP->Rebin(fBinLimits->GetSize() - 1, "histPP", fBinLimits->GetArray());
612 fHistDataPP->Scale(1., "width");
613 }
614 else
615 fHistDataPP->Clone("histPP");
616 if (fHistDataPP->GetDefaultSumw2()) fHistDataPP->Sumw2();
617 fHistDataPP->SetDirectory(0);
618 if (fRebin > 1) fHistDataPP->Rebin(fRebin);
619 }
620 // SE +-
621 fHistDataPM = (TH1*) FindObject(arrhist, PairAnalysis::EPairType::kSEPM);
622 if (fHistDataPM) {
623 if (fBinLimits) {
624 fHistDataPM = fHistDataPM->Rebin(fBinLimits->GetSize() - 1, "histPM", fBinLimits->GetArray());
625 fHistDataPM->Scale(1., "width");
626 }
627 else
628 fHistDataPM->Clone("histPM");
629 if (fHistDataPM->GetDefaultSumw2()) fHistDataPM->Sumw2();
630 fHistDataPM->SetDirectory(0);
631 if (fRebin > 1) fHistDataPM->Rebin(fRebin);
632 fHistDataPM->SetYTitle((fBinLimits ? "dN/dm" : "Counts"));
633 }
634 // SE --
635 fHistDataMM = (TH1*) FindObject(arrhist, PairAnalysis::EPairType::kSEMM);
636 if (fHistDataMM) {
637 if (fBinLimits) {
638 fHistDataMM = fHistDataMM->Rebin(fBinLimits->GetSize() - 1, "histMM", fBinLimits->GetArray());
639 fHistDataMM->Scale(1., "width");
640 }
641 else
642 fHistDataMM->Clone("histMM");
643 if (fHistDataMM->GetDefaultSumw2()) fHistDataMM->Sumw2();
644 fHistDataMM->SetDirectory(0);
645 if (fRebin > 1) fHistDataMM->Rebin(fRebin);
646 }
647 // ME ++
648 fHistMixPP = (TH1*) FindObject(arrhist, PairAnalysis::EPairType::kMEPP);
649 if (fHistMixPP) {
650 if (fBinLimits) {
651 fHistMixPP = fHistMixPP->Rebin(fBinLimits->GetSize() - 1, "mixPP", fBinLimits->GetArray());
652 fHistMixPP->Scale(1., "width");
653 }
654 else
655 fHistMixPP->Clone("mixPP");
656 if (fHistMixPP->GetDefaultSumw2()) fHistMixPP->Sumw2();
657 fHistMixPP->SetDirectory(0);
658 if (fRebin > 1) fHistMixPP->Rebin(fRebin);
659 }
660 // ME +-
661 fHistMixPM = (TH1*) FindObject(arrhist, PairAnalysis::EPairType::kMEPM);
662 if (fHistMixPM) {
663 if (fBinLimits) {
664 fHistMixPM = fHistMixPM->Rebin(fBinLimits->GetSize() - 1, "mixPM", fBinLimits->GetArray());
665 fHistMixPM->Scale(1., "width");
666 }
667 else
668 fHistMixPM->Clone("mixPM");
669 if (fHistMixPM->GetDefaultSumw2()) fHistMixPM->Sumw2();
670 fHistMixPM->SetDirectory(0);
671 if (fRebin > 1) fHistMixPM->Rebin(fRebin);
672 }
673 // ME -+
674 fHistMixMP = (TH1*) FindObject(arrhist, PairAnalysis::EPairType::kMEMP);
675 if (fHistMixMP) {
676 if (fBinLimits) {
677 fHistMixMP = fHistMixMP->Rebin(fBinLimits->GetSize() - 1, "mixMP", fBinLimits->GetArray());
678 fHistMixMP->Scale(1., "width");
679 }
680 else
681 fHistMixMP->Clone("mixMP");
682 if (fHistMixMP->GetDefaultSumw2()) fHistMixMP->Sumw2();
683 fHistMixMP->SetDirectory(0);
684 if (fRebin > 1) fHistMixMP->Rebin(fRebin);
685
686 if (fHistMixPM) fHistMixMP->Add(fHistMixPM); // merge ME +- and -+
687 }
688 // ME --
689 fHistMixMM = (TH1*) FindObject(arrhist, PairAnalysis::EPairType::kMEMM);
690 if (fHistMixMM) {
691 if (fBinLimits) {
692 fHistMixMM = fHistMixMM->Rebin(fBinLimits->GetSize() - 1, "mixMM", fBinLimits->GetArray());
693 fHistMixMM->Scale(1., "width");
694 }
695 else
696 fHistMixMM->Clone("mixMM");
697 if (fHistMixMM->GetDefaultSumw2()) fHistMixMM->Sumw2();
698 fHistMixMM->SetDirectory(0);
699 if (fRebin > 1) fHistMixMM->Rebin(fRebin);
700 }
701 // TR +-
702 fHistDataTR = (TH1*) FindObject(arrhist, PairAnalysis::EPairType::kSEPMRot);
703 if (fHistDataTR) {
704 if (fBinLimits) {
705 fHistDataTR = fHistDataTR->Rebin(fBinLimits->GetSize() - 1, "histTR", fBinLimits->GetArray());
706 fHistDataTR->Scale(1., "width");
707 }
708 else
709 fHistDataTR->Clone("histTR");
710 if (fHistDataTR->GetDefaultSumw2()) fHistDataTR->Sumw2();
711 fHistDataTR->SetDirectory(0);
712 if (fRebin > 1) fHistDataTR->Rebin(fRebin);
713 }
714
715 // init histograms for R-factor, subtracted signal, background
716 fHistSignal = new TH1D("HistSignal", "Substracted signal", fHistDataPM->GetXaxis()->GetNbins(),
717 fHistDataPM->GetXaxis()->GetXbins()->GetArray());
718 fHistSignal->SetXTitle(fHistDataPM->GetXaxis()->GetTitle());
719 fHistSignal->SetYTitle(fHistDataPM->GetYaxis()->GetTitle());
720 if (fHistSignal->GetDefaultSumw2()) fHistSignal->Sumw2();
721 fHistSignal->SetDirectory(0);
722 fHistBackground = new TH1D("HistBackground", "Background contribution", fHistDataPM->GetXaxis()->GetNbins(),
723 fHistDataPM->GetXaxis()->GetXbins()->GetArray());
724 if (fHistBackground->GetDefaultSumw2()) fHistBackground->Sumw2();
725 fHistBackground->SetDirectory(0);
726 fHistRfactor = new TH1D("HistRfactor", "Rfactor;;N^{mix}_{+-}/2#sqrt{N^{mix}_{++} N^{mix}_{--}}",
727 fHistDataPM->GetXaxis()->GetNbins(), fHistDataPM->GetXaxis()->GetXbins()->GetArray());
728 if (fHistRfactor->GetDefaultSumw2()) fHistRfactor->Sumw2();
729 fHistRfactor->SetDirectory(0);
730
731 // cocktail
732 if (fArrCocktail && fArrCocktail->GetEntriesFast()) {
733 // printf("rebin %d cocktail histograms\n",fArrCocktail->GetEntriesFast());
734 fHistCocktail = new TH1D("HistCocktail", "Cocktail contribution", fHistDataPM->GetXaxis()->GetNbins(),
735 fHistDataPM->GetXaxis()->GetXbins()->GetArray());
736 fHistCocktail->SetXTitle(fHistDataPM->GetXaxis()->GetTitle());
737 fHistCocktail->SetYTitle(fHistDataPM->GetYaxis()->GetTitle());
738 if (fHistCocktail->GetDefaultSumw2()) fHistCocktail->Sumw2();
739 fHistCocktail->SetDirectory(0);
740
741 // loop over all ingredients
742 for (Int_t i = 0; i < fArrCocktail->GetEntriesFast(); i++) {
743 TH1* htmp = static_cast<TH1*>(fArrCocktail->UncheckedAt(i));
744 if (fBinLimits) {
745 htmp = htmp->Rebin(fBinLimits->GetSize() - 1, htmp->GetTitle(), fBinLimits->GetArray());
746 if (htmp->GetDefaultSumw2()) htmp->Sumw2();
747 htmp->SetDirectory(0);
748 htmp->Scale(1., "width");
749 fArrCocktail->AddAt(htmp, i);
750 }
751 // else htmp->Clone("histTR");
752 // if(fHistDataTR->GetDefaultSumw2()) fHistDataTR->Sumw2();
753 // fHistDataTR->SetDirectory(0);
754 if (fRebin > 1) htmp->Rebin(fRebin);
755 fHistCocktail->Add(htmp, +1.);
756 }
757 }
758
760 /*
761 for(Int_t i=0; i<arrhist->GetEntriesFast(); i++) {
762 TH1* htmp = static_cast<TH1*>(arrhist->UncheckedAt(i));
763 if( htmp->GetNbinsX()!=fHistSignal->GetNbinsX() ) {
764 TArrayD *limits = PairAnalysisHelper::MakeStatBinLimits( htmp , fRebinStat);
765 htmp = htmp->Rebin(limits->GetSize()-1,htmp->GetTitle(),limits->GetArray());
766 if(htmp->GetDefaultSumw2()) htmp->Sumw2();
767 htmp->SetDirectory(0);
768 htmp->Scale(1.,"width");
769 arrhist->AddAt(htmp,i);
770 }
771 }
772 */
773
774 // process method
775 switch (fMethod) {
776 case kLikeSign:
777 case kLikeSignArithm:
778 case kLikeSignRcorr:
779 case kLikeSignArithmRcorr:
780 ProcessLS(); // process like-sign subtraction method
781 break;
782
783 case kEventMixing:
784 ProcessEM(); // process event mixing method
785 break;
786
787 case kRotation: ProcessTR(); break;
788
789 case kCocktail:
790 fCocktailSubtr = kTRUE;
791 ProcessCocktail();
792 break;
793
794 default: Warning("Process", "Subtraction method not supported. Please check SetMethod() function.");
795 }
796
797
798 // set S/B and Significance histos
799 if (!fPeakIsTF1) fHistSB = (TH1*) fgPeakShape->Clone("histSB");
800 else {
801 fHistSB = (TH1*) fHistSignal->Clone("histSB");
802 fHistSB->Reset("CEIS");
803 fHistSB->Eval((TF1*) fgPeakShape);
804 }
805 fHistSB->Divide(fHistBackground);
806 fHistSB->SetYTitle("S/B");
807
808 //significance
809 fHistSgn = (TH1*) fHistSignal->Clone("histSB");
810 fHistSgn->Reset("CEIS");
811 Double_t s = 0.;
812 Double_t b = 0.;
813 for (Int_t i = 1; i <= fHistSgn->GetNbinsX(); i++) {
814
815 if (!fPeakIsTF1) s = static_cast<TH1*>(fgPeakShape)->GetBinContent(i);
816 else
817 s = static_cast<TF1*>(fgPeakShape)->Eval(fHistSgn->GetBinCenter(i));
818 b = fHistBackground->GetBinContent(i);
819
820 if (s + b < 1.e-6) continue;
821 fHistSgn->SetBinContent(i, s / TMath::Sqrt(s + b));
822 }
823 fHistSgn->SetYTitle("S/#sqrt{S+B}");
824 // fErrors(2) = ((s+b)>0 ? fValues(2)*TMath::Sqrt(be*be + TMath::Power(se*(s+2*b)/s, 2)) / 2 / (s+b) : 0);
825}
826
827//______________________________________________
828void PairAnalysisSignalExt::ProcessLS()
829{
830 //
831 // signal subtraction suing the like-sign method
832 //
833
837
838 // protections
839 if (!fHistDataPP || !fHistDataMM) return;
840
841 // fill background histogram
842 for (Int_t ibin = 1; ibin <= fHistDataPM->GetXaxis()->GetNbins(); ibin++) {
843 Float_t pp = fHistDataPP->GetBinContent(ibin);
844 Float_t ppe = fHistDataPP->GetBinError(ibin);
845 Float_t mm = fHistDataMM->GetBinContent(ibin);
846 Float_t mme = fHistDataMM->GetBinError(ibin);
847
848 Float_t background = 2 * TMath::Sqrt(pp * mm);
849 // Float_t ebackground = TMath::Sqrt(mm+pp); // do not use it since LS could be weighted err!=sqrt(entries)
850 Float_t ebackground = TMath::Sqrt(mm / pp * ppe * ppe + pp / mm * mme * mme);
851 // Arithmetic mean instead of geometric
852 if (fMethod == kLikeSignArithm || fMethod == kLikeSignArithmRcorr) {
853 background = (pp + mm);
854 // ebackground=TMath::Sqrt(pp+mm);
855 ebackground = TMath::Sqrt(ppe * ppe + mme * mme);
856 if (TMath::Abs(ebackground) < 1e-30) ebackground = 1;
857 }
858 // set contents
859 fHistBackground->SetBinContent(ibin, background);
860 fHistBackground->SetBinError(ibin, ebackground);
861 }
862
863 //correct LS spectrum bin-by-bin with R factor obtained in mixed events
864 if (fMixingCorr || fMethod == kLikeSignRcorr || fMethod == kLikeSignArithmRcorr) {
865
866 // protections
867 if (!fHistMixPM || fHistMixPP || !fHistMixMM) return;
868
869 // fill R-factor histogram
870 for (Int_t ibin = 1; ibin <= fHistMixPM->GetXaxis()->GetNbins(); ibin++) {
871 Float_t pp = fHistMixPP->GetBinContent(ibin);
872 Float_t ppe = fHistMixPP->GetBinError(ibin);
873 Float_t mm = fHistMixMM->GetBinContent(ibin);
874 Float_t mme = fHistMixMM->GetBinError(ibin);
875 Float_t pm = fHistMixPM->GetBinContent(ibin);
876 Float_t pme = fHistMixPM->GetBinError(ibin);
877
878 Float_t background = 2 * TMath::Sqrt(pp * mm);
879 // do not use it since ME could be weighted err!=sqrt(entries)
880 // Float_t ebackground = TMath::Sqrt(mm+pp);
881 Float_t ebackground = TMath::Sqrt(mm / pp * ppe * ppe + pp / mm * mme * mme);
882 //Arithmetic mean instead of geometric
883 if (fMethod == kLikeSignArithm) {
884 background = (pp + mm);
885 ebackground = TMath::Sqrt(ppe * ppe + mme * mme);
886 if (TMath::Abs(ebackground) < 1e-30) ebackground = 1;
887 }
888
889 Float_t rcon = 1.0;
890 Float_t rerr = 0.0;
891 if (background > 0.0) {
892 rcon = pm / background;
893 rerr = TMath::Sqrt((1. / background) * (1. / background) * pme * pme
894 + (pm / (background * background)) * (pm / (background * background)) * ebackground
895 * ebackground);
896 }
897 fHistRfactor->SetBinContent(ibin, rcon);
898 fHistRfactor->SetBinError(ibin, rerr);
899 }
900 // correction
901 fHistBackground->Multiply(fHistRfactor);
902 }
903
904 //scale histograms to match integral between fScaleMin and fScaleMax
905 // or if fScaleMax < fScaleMin use fScaleMin as scale factor
906 ScaleBackground();
907
908 //subtract background
909 fHistSignal->Add(fHistDataPM);
910 fHistSignal->Add(fHistBackground, -1);
911
912 //subtract cocktail (if added)
913 if (fCocktailSubtr) fHistSignal->Add(fHistCocktail, -1);
914
915 // background
916 fValues(1) =
917 fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin), fHistBackground->FindBin(fIntMax), fErrors(1));
918
919 // signal depending on peak description method
920 DescribePeakShape(fPeakMethod, kTRUE, fgHistSimPM);
921 //printf("%f %f\n",fValues(0),fValues(1));
922 // SetSignificanceAndSOB();
923
924 fProcessed = kTRUE;
925}
926
927//______________________________________________
928void PairAnalysisSignalExt::ProcessEM()
929{
930 //
931 // signal subtraction using event mixing (+-,-+) method
932 //
933
934 if (!fHistMixMP) {
935 Error("ProcessEM", "Mixed event histogram missing");
936 return;
937 }
938
939
940 // fill background histogram
941 for (Int_t ibin = 1; ibin <= fHistDataPM->GetXaxis()->GetNbins(); ibin++) {
942 Float_t pm = fHistMixMP->GetBinContent(ibin);
943 Float_t pme = fHistMixMP->GetBinError(ibin);
944
945 Float_t background = pm;
946 Float_t ebackground = TMath::Sqrt(pme * pme);
947 //TODO: Float_t ebackground = TMath::Power(pm, 3./4);
948
949 // set contents
950 fHistBackground->SetBinContent(ibin, background);
951 fHistBackground->SetBinError(ibin, ebackground);
952 }
953
954 //scale histograms to match integral between fScaleMin and fScaleMax
955 // or if fScaleMax < fScaleMin use fScaleMin as scale factor
956 ScaleBackground();
957
958 //subtract background
959 fHistSignal->Add(fHistDataPM);
960 fHistSignal->Add(fHistBackground, -1);
961 // fHistSignal = MergeObjects(fHistDataPM,fHistBackground,-1.);
962
963 //subtract cocktail (if added)
964 if (fCocktailSubtr) fHistSignal->Add(fHistCocktail, -1);
965
966 // background
967 fValues(1) =
968 fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin), fHistBackground->FindBin(fIntMax), fErrors(1));
969
970 // signal depending on peak description method
971 DescribePeakShape(fPeakMethod, kTRUE, fgHistSimPM);
972
973 fProcessed = kTRUE;
974}
975
976//______________________________________________
977void PairAnalysisSignalExt::ProcessTR()
978{
979 //
980 // signal subtraction using the track-rotation method
981 //
982
983 if (!fHistDataTR) {
984 Error("ProcessTR", "Track rotation histogram missing");
985 return;
986 }
987
988 // fill background histogram
989 for (Int_t ibin = 1; ibin <= fHistDataPM->GetXaxis()->GetNbins(); ibin++) {
990 Float_t pm = fHistDataTR->GetBinContent(ibin);
991 Float_t pme = fHistDataTR->GetBinError(ibin);
992
993 Float_t background = pm;
994 Float_t ebackground = TMath::Sqrt(pme * pme);
995
996 // set contents
997 fHistBackground->SetBinContent(ibin, background);
998 fHistBackground->SetBinError(ibin, ebackground);
999 }
1000
1001 // scale according to number of interations used in track rotation
1002 fHistBackground->Scale(1. / fNiterTR);
1003
1004 //scale histograms to match integral between fScaleMin and fScaleMax
1005 // or if fScaleMax < fScaleMin use fScaleMin as scale factor
1006 ScaleBackground();
1007
1008 //subtract background
1009 fHistSignal->Add(fHistDataPM);
1010 fHistSignal->Add(fHistBackground, -1);
1011
1012 //subtract cocktail (if added)
1013 if (fCocktailSubtr) fHistSignal->Add(fHistCocktail, -1);
1014
1015 // background
1016 fValues(1) =
1017 fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin), fHistBackground->FindBin(fIntMax), fErrors(1));
1018
1019 // signal depending on peak description method
1020 DescribePeakShape(fPeakMethod, kTRUE, fgHistSimPM);
1021
1022 fProcessed = kTRUE;
1023}
1024
1025//______________________________________________
1026void PairAnalysisSignalExt::ProcessCocktail()
1027{
1028 //
1029 // signal subtraction using the cocktail method
1030 //
1031
1032 if (!fArrCocktail) {
1033 Error("ProcessCocktail", "Cocktail histograms missing");
1034 return;
1035 }
1036
1037 // fill background histogram
1038 if (fArrCocktail) fHistBackground = fHistCocktail;
1039
1040 //scale histograms to match integral between fScaleMin and fScaleMax
1041 // or if fScaleMax < fScaleMin use fScaleMin as scale factor
1042 ScaleBackground();
1043
1044 //subtract background
1045 fHistSignal->Add(fHistDataPM);
1046 fHistSignal->Add(fHistBackground, -1);
1047
1048 // background
1049 fValues(1) =
1050 fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin), fHistBackground->FindBin(fIntMax), fErrors(1));
1051
1052 // signal depending on peak description method
1053 DescribePeakShape(fPeakMethod, kTRUE, fgHistSimPM);
1054
1055 fProcessed = kTRUE;
1056}
1057
1058//______________________________________________
1059void PairAnalysisSignalExt::Draw(const Option_t* option)
1060{
1083
1084 Info("Draw", "Signal extraction results for '%s'", fgkBackgroundMethodNames[fMethod]);
1085 TString optString(option);
1086 optString.ToLower();
1087 printf("Plot extraction: bgrd. estimator: '%s' \t options: '%s' \n", fgkBackgroundMethodNames[fMethod],
1088 optString.Data());
1089 Bool_t optTask = optString.Contains("same");
1090 optString.ReplaceAll("same", "");
1091 Bool_t optLegFull = optString.Contains("legf");
1092 optString.ReplaceAll("legf", "");
1093 Bool_t optLeg = optString.Contains("leg");
1094 optString.ReplaceAll("leg", "");
1095 Bool_t optCan = optString.Contains("can");
1096 optString.ReplaceAll("can", "");
1097 Bool_t optLine = optString.Contains("line");
1098 optString.ReplaceAll("line", "");
1099 Bool_t optStat = optString.Contains("stat");
1100 optString.ReplaceAll("stat", "");
1101 Bool_t optSSB = optString.Contains("ssb");
1102 optString.ReplaceAll("ssb", "");
1103 Bool_t optSB = optString.Contains("sb");
1104 optString.ReplaceAll("sb", "");
1105 Bool_t optSgn = optString.Contains("sgn");
1106 optString.ReplaceAll("sgn", "");
1107 Bool_t optOnlyRaw = optString.Contains("onlyraw");
1108 optString.ReplaceAll("onlyraw", "");
1109 Bool_t optOnlySig = optString.Contains("onlysig");
1110 optString.ReplaceAll("onlysig", "");
1111 Bool_t optNoSig = optString.Contains("nosig");
1112 optString.ReplaceAll("nosig", "");
1113 Bool_t optNoBgrd = optString.Contains("nobgrd");
1114 optString.ReplaceAll("nobgrd", "");
1115 Bool_t optOnlyMC = optString.Contains("onlymc");
1116 optString.ReplaceAll("onlymc", "");
1117 Bool_t optNoMC = optString.Contains("nomc");
1118 optString.ReplaceAll("nomc", "");
1119 Bool_t optCocktail = optString.Contains("cocktail");
1120 optString.ReplaceAll("cocktail", "");
1121 optString.ReplaceAll(" ", "");
1122
1124 PairAnalysisStyler::LoadStyle();
1125
1127 if (optLegFull) optLeg = kTRUE;
1128
1130 TCanvas* c = 0;
1131 if (optCan) {
1132 TString key = Form("cSignalExtraction%s", (optSB ? "SB" : (optSgn ? "SGN" : (optSSB ? "SSB" : ""))));
1133 c = (TCanvas*) gROOT->FindObject(key.Data());
1134 if (!c) c = new TCanvas(key.Data(), key.Data());
1135 // c->Clear();
1136 c->cd();
1137 }
1138
1139 Int_t nobj = 0;
1140 TObject* obj;
1142 TList* prim = gPad->GetListOfPrimitives();
1143 for (Int_t io = 0; io < prim->GetSize(); io++) {
1144 obj = prim->At(io);
1145 if (obj->InheritsFrom(TH1::Class()) && obj != prim->At(io + 1)) nobj++;
1146 }
1147
1148 // add or get legend
1149 TLegend* leg = 0;
1150 if ((optLeg && optTask && !nobj) || (optLeg && !optTask)) {
1151 leg = new TLegend(.75, .3, 1. - gPad->GetTopMargin() - gStyle->GetTickLength("X"),
1152 1. - gPad->GetRightMargin() - gStyle->GetTickLength("Y"), fArrHists->GetName(), "nbNDC");
1153 if (optTask) leg->SetHeader("");
1154 }
1155 else if (optLeg && nobj) {
1156 leg = (TLegend*) prim->FindObject("TPave");
1157 }
1158
1159 // logaritmic style
1160 gPad->SetLogx(optString.Contains("logx"));
1161 gPad->SetLogy(optString.Contains("logy"));
1162 gPad->SetLogz(optString.Contains("logz"));
1163 optString.ReplaceAll("logx", "");
1164 optString.ReplaceAll("logy", "");
1165 optString.ReplaceAll("logz", "");
1166
1168 if (optString.Contains("e2") || optString.Contains("e3") || optString.Contains("e4")) {
1169 // printf("set x-error style \n");
1170 gStyle->SetErrorX(0.5);
1171 // PairAnalysisStyler::SetForceFillStyle(1);
1172 }
1173
1174 Int_t i = nobj; // TODO: obsolete?
1175 Info("Draw", "this is object number: %d", nobj);
1176
1177
1178 TString ytitle = fHistDataPM->GetYaxis()->GetTitle();
1179 // add bin width to y-axis title
1180 if (!(fHistDataPM->GetXaxis()->IsVariableBinSize())) {
1181 ytitle += Form("/%.0f MeV/#it{c}^{2}", fHistDataPM->GetBinWidth(1) * 1e3);
1182 }
1183 fHistDataPM->GetYaxis()->SetTitle(ytitle.Data());
1184 fHistSignal->GetYaxis()->SetTitle(ytitle.Data());
1185
1186 // styling
1187 fHistDataPM->SetNameTitle(Form("unlike-sign%s", GetName()), "unlike-sign");
1188 fHistDataPM->UseCurrentStyle();
1189 PairAnalysisStyler::Style(fHistDataPM, static_cast<Int_t>(PairAnalysisStyler::Eidx::kRaw));
1190 if (fPlotMin != fPlotMax) fHistDataPM->SetAxisRange(fPlotMin, fPlotMax, "X");
1191
1192 fHistBackground->SetNameTitle(fgkBackgroundMethodNames[fMethod], fgkBackgroundMethodNames[fMethod]);
1193 fHistBackground->UseCurrentStyle();
1194 PairAnalysisStyler::Style(fHistBackground, static_cast<Int_t>(PairAnalysisStyler::Eidx::kBgrd));
1195
1196 fHistSignal->SetNameTitle(Form("signal%s", GetName()), "signal");
1197 fHistSignal->UseCurrentStyle();
1198 PairAnalysisStyler::Style(fHistSignal, static_cast<Int_t>(PairAnalysisStyler::Eidx::kSig));
1199 if (fPlotMin != fPlotMax) fHistSignal->SetAxisRange(fPlotMin, fPlotMax, "X");
1200
1201 if (fHistCocktail) {
1202 fHistCocktail->SetNameTitle(Form("cocktail%s", GetName()), "cocktail");
1203 fHistCocktail->UseCurrentStyle();
1204 PairAnalysisStyler::Style(fHistCocktail, static_cast<Int_t>(PairAnalysisStyler::Eidx::kCocktail));
1205 if (fPlotMin != fPlotMax) fHistCocktail->SetAxisRange(fPlotMin, fPlotMax, "X");
1206 }
1207
1208 if (optSB) {
1209 fHistSB->SetNameTitle(Form("s2b%s", GetName()), "signal");
1210 fHistSB->UseCurrentStyle();
1211 PairAnalysisStyler::Style(fHistSB, static_cast<Int_t>(PairAnalysisStyler::Eidx::kSig));
1212 if (fPlotMin != fPlotMax) fHistSB->SetAxisRange(fPlotMin, fPlotMax, "X");
1213 }
1214 else if (optSgn) {
1215 fHistSgn->SetNameTitle(Form("sgn%s", GetName()), "signal");
1216 fHistSgn->UseCurrentStyle();
1217 PairAnalysisStyler::Style(fHistSgn, static_cast<Int_t>(PairAnalysisStyler::Eidx::kSig));
1218 if (fPlotMin != fPlotMax) fHistSgn->SetAxisRange(fPlotMin, fPlotMax, "X");
1219 }
1220
1221 // fHistRfactor->UseCurrentStyle();
1222 // fHistRfactor->SetTitle("");
1223 PairAnalysisStyler::Style(fHistRfactor, 3);
1224
1225 fgPeakShape->UseCurrentStyle();
1226 // fgPeakShape->SetTitle("");
1227 PairAnalysisStyler::Style(fgPeakShape, static_cast<Int_t>(PairAnalysisStyler::Eidx::kFit));
1228
1229 // draw stuff
1230 if (c) c->cd(1);
1231 Info("Draw", "Start plotting stuff with option -%s-", optString.Data());
1232 TString drawOpt = (optString.IsNull() ? "EP" : optString);
1233 // draw spectra
1234 if (!optOnlyMC) {
1235
1236 if (optSB && !optOnlyRaw && !optNoSig) {
1237 fHistSB->Draw(i > 0 ? (drawOpt + "same").Data() : drawOpt.Data());
1238 i++;
1239 }
1240 else if (optSgn && !optOnlyRaw && !optNoSig) {
1241 fHistSgn->Draw(i > 0 ? (drawOpt + "same").Data() : drawOpt.Data());
1242 i++;
1243 }
1244 else if (optOnlySig) {
1245 drawOpt = (optString.IsNull() ? "EP0" : optString);
1246 fHistSignal->DrawCopy(i > 0 ? (drawOpt + "same").Data() : drawOpt.Data());
1247 i++;
1248 drawOpt = (optString.IsNull() ? "L same" : optString + "same");
1249 if (fPeakIsTF1) {
1250 static_cast<TF1*>(fgPeakShape)->DrawCopy(drawOpt.Data());
1251 i++;
1252 }
1253 }
1254 else if (!optSB && !optSgn) {
1255 fHistDataPM->DrawCopy(i > 0 ? (drawOpt + "same").Data() : drawOpt.Data());
1256 i++;
1257 if (fMethod == kRotation || fMethod == kEventMixing || fMethod == kCocktail)
1258 drawOpt = (optString.IsNull() ? "HIST" : optString);
1259 if (!optNoBgrd) {
1260 fHistBackground->DrawCopy(i > 0 ? (drawOpt + "same").Data() : drawOpt.Data());
1261 i++;
1262 }
1263 if (!optOnlyRaw && !optNoSig) {
1264 drawOpt = (optString.IsNull() ? "EP0" : optString);
1265 fHistSignal->DrawCopy(i > 0 ? (drawOpt + "same").Data() : drawOpt.Data());
1266 i++;
1267 drawOpt = (optString.IsNull() ? "L same" : optString + "same");
1268 if (fPeakIsTF1) {
1269 static_cast<TF1*>(fgPeakShape)->DrawCopy(drawOpt.Data());
1270 i++;
1271 }
1272 }
1273 }
1274 // add cocktail
1275 if (optCocktail && fHistCocktail) {
1276 drawOpt = (optString.IsNull() ? "HIST" : optString);
1277 if (optSSB) {
1278 fHistCocktail->Divide(fHistDataPM);
1279 fHistCocktail->SetYTitle(Form("S/(S+B)"));
1280 }
1281 if (optSB) {
1282 fHistCocktail->Divide(fHistBackground);
1283 fHistCocktail->SetYTitle(Form("%s", GetValueName(3)));
1284 }
1285 if (optSgn) { FillSignificance(fHistCocktail, fHistCocktail, fHistBackground); }
1286 fHistCocktail->DrawCopy(i > 0 ? (drawOpt + "same").Data() : drawOpt.Data());
1287 i++; //TODO: needed?
1288 // fHistCocktail->DrawCopy((drawOpt+"same").Data()); i++;
1289 }
1290 }
1291
1292 // draw MC signals
1293 if (!optNoMC) {
1294 if (!optSgn) {
1295 TIter nextObj(fArrHists);
1296 TH1* hmc;
1297 Int_t isty = 0;
1298 while ((hmc = dynamic_cast<TH1*>(nextObj()))) {
1299 TString key = hmc->GetTitle(); //hmc->GetName();
1300 TString tit = hmc->GetTitle();
1301 if (key.CountChar('_') != 1) continue; // only reconstr. MC signals
1302 // remove cocktail subtracted signals
1303 if (optOnlySig && fCocktailSubtr && FindObjectByTitle(fArrCocktail, key)) continue;
1304 else if (optOnlySig && !fCocktailSubtr && fArrCocktail && !FindObjectByTitle(fArrCocktail, key))
1305 continue;
1306 else if (optCocktail && FindObjectByTitle(fArrCocktail, key))
1307 continue;
1308 PairAnalysisStyler::Style(hmc, isty++);
1309 // check if rebinning is necessary
1310 if (1 && fHistSignal->GetNbinsX() != hmc->GetNbinsX()) {
1311 if (fBinLimits) {
1312 hmc = hmc->Rebin(fBinLimits->GetSize() - 1, key.Data(), fBinLimits->GetArray());
1313 hmc->SetTitle(tit.Data());
1314 hmc->Scale(1., "width");
1315 }
1316 if (fRebin > 1) hmc->Rebin(fRebin);
1317 }
1318 if (optSSB) {
1319 hmc->Divide(fHistDataPM);
1320 hmc->SetYTitle(Form("S/(S+B)"));
1321 }
1322 if (optSB) {
1323 hmc->Divide(fHistBackground);
1324 hmc->SetYTitle(Form("%s", GetValueName(3)));
1325 }
1326 if (fPlotMin != fPlotMax) hmc->SetAxisRange(fPlotMin, fPlotMax, "X");
1327 hmc->Draw(i > 0 ? "HISTsame" : "HIST");
1328 i++;
1329 }
1330 }
1331 }
1332
1334 Double_t max = -1e10;
1335 Double_t min = +1e10;
1336 //TListIter nextObj(gPad->GetListOfPrimitives(),kIterForward);
1337 TListIter nextObj(gPad->GetListOfPrimitives(), kIterBackward);
1338 // TObject *obj;
1339 while ((obj = nextObj())) {
1340 if (obj->InheritsFrom(TH1::Class())) {
1341 TH1* hobj = static_cast<TH1*>(obj);
1342 max = TMath::Max(max,
1343 PairAnalysisHelper::GetContentMaximum(hobj)); //hobj->GetMaximum();
1344 Double_t tmpmax = max * (gPad->GetLogy() ? 5. : 1.1);
1345 hobj->SetMaximum(tmpmax);
1346
1347 Double_t objmin = PairAnalysisHelper::GetContentMinimum(hobj, kTRUE);
1348 if (gPad->GetLogy() && objmin < 0.) objmin = 0.5;
1349 min = TMath::Min(min, objmin);
1350 Double_t tmpmin = min * (min < 0. ? 1.1 : 0.9);
1351 hobj->SetMinimum(tmpmin);
1352
1353 // Printf("after %s max%f \t min%f \t for logy %.3e >? %.3e",
1354 // hobj->GetTitle(),tmpmax,tmpmin, tmpmax/(tmpmin>0.?tmpmin:1.),TMath::Power(10.,TGaxis::GetMaxDigits()));
1355 // automatically set log option
1356 if (gPad->GetLogy()
1357 && (tmpmax / (tmpmin > 0. ? tmpmin : 1.) > TMath::Power(10., TGaxis::GetMaxDigits())
1358 || tmpmin < TMath::Power(10., -TGaxis::GetMaxDigits()))) {
1359 hobj->GetYaxis()->SetMoreLogLabels(kFALSE);
1360 hobj->GetYaxis()->SetNoExponent(kFALSE);
1361 }
1362 // if(gPad->GetLogy()) break;
1363 }
1364 }
1365
1366 // add legend entries
1367 TListIter nextObjFwd(gPad->GetListOfPrimitives(), kIterForward);
1368 if (optLeg && leg) {
1369 nextObjFwd.Reset();
1370 Int_t ileg = 0;
1371 while ((obj = nextObjFwd())) {
1372 if (obj->InheritsFrom(TH1::Class()) || obj->InheritsFrom(TF1::Class())) {
1373 // TH1 *hleg = static_cast<TH1*>(obj);
1374 if (nobj && ileg < nobj) {
1375 ileg++;
1376 continue;
1377 }
1378
1379 TString histClass = obj->GetTitle(); //hleg->GetName();
1380 if (histClass.IsNull()) histClass = obj->GetName();
1381 histClass.ReplaceAll("Pair.", "");
1382 histClass.ReplaceAll("Pair_", "");
1383 // change default signal names to titles
1384 for (Int_t isig = 0; isig < static_cast<Int_t>(PairAnalysisSignalMC::EDefinedSignal::kNSignals); isig++) {
1385 histClass.ReplaceAll(PairAnalysisSignalMC::fgkSignals[isig][0], PairAnalysisSignalMC::fgkSignals[isig][1]);
1386 }
1387 // remove trailing and leading spaces
1388 histClass.Remove(TString::kBoth, ' ');
1389 // histClass.Remove(TString::kBoth,'.');
1390 // histClass.Remove(TString::kBoth,'_');
1391
1392 // modify legend option
1393 TString legOpt = obj->GetDrawOption();
1394 legOpt.ToLower();
1395 legOpt += "l";
1396 legOpt.ReplaceAll("hist", "");
1397 legOpt.ReplaceAll("same", "");
1398 legOpt.ReplaceAll("scat", "");
1399 legOpt.ReplaceAll("col", "");
1400 legOpt.ReplaceAll("z", "");
1401 legOpt.ReplaceAll("text", "");
1402 legOpt.ReplaceAll("e", "");
1403 // Printf("hist %s legopt %s",histClass.Data(),legOpt.Data());
1404 if (ileg == nobj && optTask) leg->AddEntry((TObject*) 0x0, fArrHists->GetName(), "");
1405 leg->AddEntry(obj, histClass.Data(), legOpt.Data());
1406 ileg++;
1407 }
1408 }
1409 }
1410
1413 if (leg) {
1414 PairAnalysisStyler::SetLegendAttributes(leg, optLegFull);
1415 if (!nobj) leg->Draw(); // was w/o !nobj
1416 }
1417
1418 // baseline
1419 TLine* line = new TLine();
1420 line->SetLineColor(kBlack);
1421 line->SetLineStyle(kDashed);
1422 if (optLine) line->DrawLine(fPlotMin, 0., fPlotMax, 0.);
1423
1424 // draw statistics box
1425 if (optStat)
1426 DrawStats(0.7, gPad->GetBottomMargin() + gStyle->GetTickLength("Y"),
1427 1. - gPad->GetRightMargin() - gStyle->GetTickLength("X"),
1428 gPad->GetBottomMargin() + gStyle->GetTickLength("Y") + 5 * 0.025);
1429
1430
1431 // styling
1432 // gPad->RedrawAxis();
1433}
1434
1435void PairAnalysisSignalExt::ScaleBackground()
1436{
1437 //
1438 // scale histograms according to method and integral limits
1439 //
1440
1441 TH1* scalingRef;
1442 // scalingRef->Sumw2();
1443 switch (fSclMethod) {
1444 case kSclToRaw: scalingRef = fHistDataPM; break;
1445 case kSclToLikeSign:
1446 scalingRef = (TH1*) fHistDataPP->Clone();
1447 scalingRef->Add(fHistDataMM);
1448 break; //arithmetic mean
1449 default: scalingRef = nullptr; break;
1450 }
1451 if (fScaleMax > fScaleMin && fScaleMax2 > fScaleMin2)
1452 fScaleFactor = ScaleHistograms(scalingRef, fHistBackground, fScaleMin, fScaleMax, fScaleMin2, fScaleMax2);
1453 else if (fScaleMax > fScaleMin)
1454 fScaleFactor = ScaleHistograms(scalingRef, fHistBackground, fScaleMin, fScaleMax);
1455 else if (fScaleMin > 0.) {
1456 fScaleFactor = fScaleMin;
1457 fHistBackground->Scale(fScaleFactor);
1458 }
1459
1460 // clean up
1461 if (fSclMethod == kSclToLikeSign && scalingRef) delete scalingRef;
1462}
ClassImp
ClassImp(CbmConverterManager)
PairAnalysisFunction::PeakFunGauss
Double_t PeakFunGauss(const Double_t *x, const Double_t *par)
Definition PairAnalysisFunction.cxx:206
PairAnalysisFunction::GetCombinedFunction
TF1 * GetCombinedFunction() const
Definition PairAnalysisFunction.h:92
PairAnalysisFunction::PeakFunMC
Double_t PeakFunMC(const Double_t *x, const Double_t *par)
Definition PairAnalysisFunction.cxx:105
PairAnalysisFunction::PeakFunCB
Double_t PeakFunCB(const Double_t *x, const Double_t *par)
Definition PairAnalysisFunction.cxx:126
PairAnalysisSignalExt::fPeakMethod
ESignalExtractionMethod fPeakMethod
Definition PairAnalysisSignalExt.h:254
PairAnalysisSignalExt::Print
void Print(Option_t *option="") const
Definition PairAnalysisSignalExt.cxx:198
PairAnalysisSignalExt::GetValueName
static const char * GetValueName(Int_t i)
Definition PairAnalysisSignalExt.h:150
PairAnalysisSignalExt::Process
void Process(TObjArray *const arrhist)
Definition PairAnalysisSignalExt.cxx:523
PairAnalysisSignalExt::Draw
void Draw(const Option_t *option="")
Definition PairAnalysisSignalExt.cxx:1059
PairAnalysisSignalExt::ScaleHistograms
Double_t ScaleHistograms(TH1 *histRaw, TH1 *histBackground, Double_t intMin, Double_t intMax)
Definition PairAnalysisSignalExt.cxx:215
PairAnalysisSignalExt::DrawStats
TPaveText * DrawStats(Double_t x1=0., Double_t y1=0., Double_t x2=0., Double_t y2=0., TString opt="pRnbsSmrc")
Definition PairAnalysisSignalExt.cxx:143
PairAnalysisSignalExt::fExtrFunc
PairAnalysisFunction * fExtrFunc
Definition PairAnalysisSignalExt.h:253
PairAnalysisSignalExt::fgkValueNames
static const char * fgkValueNames[7]
Definition PairAnalysisSignalExt.h:51
PairAnalysisSignalExt::FindObject
TObject * FindObject(TObjArray *arrhist, PairAnalysis::EPairType type) const
Definition PairAnalysisSignalExt.h:274
PairAnalysisSignalExt::FillSignificance
void FillSignificance(TH1 *hfill, TObject *signal, TH1 *hbgrd)
Definition PairAnalysisSignalExt.h:371
PairAnalysisSignalExt::MergeObjects
static TH1 * MergeObjects(TH1 *obj1, TH1 *obj2, Double_t val=+1.)
Definition PairAnalysisSignalExt.cxx:273
PairAnalysisSignalExt::FindObjectByTitle
TObject * FindObjectByTitle(TObjArray *arrhist, TString ref)
Definition PairAnalysisSignalExt.h:288
PairAnalysisSignalExt::fgkBackgroundMethodNames
static const char * fgkBackgroundMethodNames[11]
Definition PairAnalysisSignalExt.h:53
PairAnalysisSignalExt::DescribePeakShape
TObject * DescribePeakShape(ESignalExtractionMethod method=kMCFitted, Bool_t replaceValErr=kFALSE, TH1F *mcShape=0x0)
Definition PairAnalysisSignalExt.cxx:337
PairAnalysisSignalMC::fgkSignals
static const char * fgkSignals[static_cast< int >(EDefinedSignal::kNSignals)][2]
Definition PairAnalysisSignalMC.h:204
PairAnalysisHelper::MakeStatBinLimits
TArrayD * MakeStatBinLimits(TH1 *h, Double_t stat)
Definition PairAnalysisHelper.cxx:198
PairAnalysisHelper::GetContentMinimum
Double_t GetContentMinimum(TH1 *h, Bool_t inclErr=kTRUE)
Definition PairAnalysisHelper.cxx:479
PairAnalysisHelper::GetContentMaximum
Double_t GetContentMaximum(TH1 *h, Bool_t inclErr=kTRUE)
Definition PairAnalysisHelper.cxx:510
PairAnalysisStyler::Style
void Style(TObject *obj, Int_t idx=0)
Definition PairAnalysisStyler.cxx:261
PairAnalysisStyler::SetLegendAttributes
void SetLegendAttributes(TLegend *leg, Bool_t fill=kFALSE)
Definition PairAnalysisStyler.cxx:466