CbmRoot
PairAnalysisSpectrum.cxx
Go to the documentation of this file.
1 // PairAnalysisSpectrum //
3 // //
4 // //
5 /*
6 
7  post-processing class to extract signals, efficiency, apply corrections,
8  calculate systematics and describe the spectra by some functions/models.
9 
10  Example:
11 
12  PairAnalysisSpectrum *spectrum = new PairAnalysisSpectrum("legend header","systematics");
13 
14  spectrum->SetParticleOfInterest(pdgcode);
15  spectrum->SetVariable("Pt",PairAnalysisHelper::MakeLinBinning(20,0.,2.) );
16 
17  spectrum->SetSystMethod( PairAnalysisSpectrum::kSystMax );
18 
19  // add input spectra coming from PairAnalysisHistos
20  spectrum->AddInput( histos->DrawSame("pM_Pt","nomc goff"), // raw invariant mass spectrum
21  "like-sign", // unique string
22  histos->DrawSame("pPt","onlymc goff sel","phi") // optional: MC spectra for efficiency
23  );
24  spectrum->AddExtractor( sig ); // signal extraction (see PairAnalysisSignalExt)
25 
26  .... add more input as much as you want
27 
28 
29  // process all inputs
30  spectrum->Process();
31 
32 
33  // draw spectra using TTree::Draw command + some extra arguments (see Draw)
34  // see Extraction for content
35  spectrum->Draw("s/eff:var","","leg logY syst P");
36 
37  // fit the spectrum by some predefined function (see PairAnalysisFunction)
38  spectrum->SetFitRange(0.,2.);
39  spectrum->SetDefault( PairAnalysisFunction::kBoltzmann );
40  spectrum->SetFitOption("RN0");
41  spectrum->Fit("leg L");
42 
43  */
44 // //
46 
47 //#include <TObject.h>
48 #include <TEventList.h>
49 #include <TROOT.h>
50 #include <TTree.h>
51 
52 #include <TCanvas.h>
53 #include <TPad.h>
54 
55 #include <TDatabasePDG.h>
56 #include <TF1.h>
57 #include <TFile.h>
58 #include <TFormula.h>
59 #include <TGraphErrors.h>
60 #include <TH1.h>
61 #include <TList.h>
62 #include <TPaveText.h>
63 #include <TProfile.h>
64 #include <TProfile2D.h>
65 #include <TVectorT.h>
66 
67 #include "PairAnalysisHF.h"
68 #include "PairAnalysisHistos.h"
69 #include "PairAnalysisStyler.h"
70 #include "PairAnalysisVarManager.h"
71 
72 
73 #include "PairAnalysisSignalExt.h"
74 
75 #include "PairAnalysisSpectrum.h"
76 
77 ClassImp(PairAnalysisSpectrum)
78 
79  //______________________________________________
80  PairAnalysisSpectrum::PairAnalysisSpectrum()
81  : PairAnalysisSpectrum("spectrum", "title") {
85 }
86 
87 //______________________________________________
88 PairAnalysisSpectrum::PairAnalysisSpectrum(const char* name, const char* title)
89  : PairAnalysisFunction(name, title)
90  , fRawInput(0)
91  , fMCInput(0)
92  , fMCTruth(0)
93  , fExtractor(0) {
97  fTree = new TTree("PAPa", "PAPa-Spectrum");
98  fExt = new Extraction;
99  for (Int_t i = 0; i < 100; i++)
100  fInputKeys[i] = "";
101 }
102 
103 //______________________________________________
104 PairAnalysisSpectrum::~PairAnalysisSpectrum() {
108  if (fResults) delete fResults;
109  if (fExtractions) delete fExtractions;
110 }
111 
112 //______________________________________________
113 void PairAnalysisSpectrum::AddInput(TObjArray* raw,
114  TString identifier,
115  TObjArray* mc,
116  TObjArray* truth) {
121  fRawInput.Add(raw);
122  if (mc) fMCInput.Add(mc);
123  if (truth) fMCTruth.Add(truth);
124  fInputKeys[fIdx] = identifier;
125  fIdx++;
126 }
127 
128 //______________________________________________
129 void PairAnalysisSpectrum::Init() {
133  fTree->Branch("Extraction", &fExt);
134  // fTree ->Branch("hSignal","TH1",&fHistSignal,32000,0);
135 }
136 
137 //______________________________________________
138 void PairAnalysisSpectrum::Process() {
142  Init();
143 
144  Int_t i = -1;
145 
146  TIter nextRaw(&fRawInput); // raw content
147  TObject* obj = NULL;
148  PairAnalysisHF* hf = NULL;
149  PairAnalysisHistos* h = NULL;
150 
151  TObject* objMC = NULL;
152  PairAnalysisHF* hfMC = NULL;
153  PairAnalysisHistos* hMC = NULL;
154 
155  TObject* objMCtruth = NULL;
156 
158  while ((obj = nextRaw())) {
159  i++;
160  TObject::Info("Process", "Check Extraction for %s", fInputKeys[i].Data());
161  TObject::Info("Process",
162  "------------------------------------------------------------"
163  "-----------");
164  TObject::Info("Process", "Input type: %s \n", obj->ClassName());
165 
166  if (fVarBinning) {
167  TObject::Info("Process",
168  "Binning provided for %s from %.3f to %.3f",
169  fVar.Data(),
170  fVarBinning->Min(),
171  fVarBinning->Max());
172  // fVarBinning->Print();
173  }
174 
175  TObjArray* histArr = dynamic_cast<TObjArray*>(obj);
176  if (!(histArr)) {
177  if (!(h = dynamic_cast<PairAnalysisHistos*>(obj))) {
178  if (!(hf = dynamic_cast<PairAnalysisHF*>(obj))) {
179  TObject::Error("Process", "No input format found");
180  continue;
181  }
182  }
183  }
184 
185  // get extractor
186  PairAnalysisSignalExt* sig =
187  dynamic_cast<PairAnalysisSignalExt*>(fExtractor.At(i));
188  if (!sig) continue;
189 
190 
192  objMC = fMCInput.At(i);
193  TObjArray* histArrMC = NULL;
194  if (objMC) {
195  TObject::Info("Process", "Input MC type: %s \n", objMC->ClassName());
196  histArrMC = dynamic_cast<TObjArray*>(objMC);
197  if (!(histArrMC)) {
198  if (!(hMC = dynamic_cast<PairAnalysisHistos*>(objMC))) {
199  if (!(hfMC = dynamic_cast<PairAnalysisHF*>(objMC))) {
200  TObject::Error("Process", "No MC input format found");
201  //continue;
202  }
203  }
204  }
205  }
206 
208  objMCtruth = fMCTruth.At(i);
209  TObjArray* histArrMCt = NULL;
210  if (objMCtruth) {
211  TObject::Info(
212  "Process", "Input MC truth type: %s \n", objMCtruth->ClassName());
213  histArrMCt = dynamic_cast<TObjArray*>(objMCtruth);
214  }
215 
216 
218  if (!fVarBinning) {
219 
220  // get raw input histograms
221  if (!histArr && h)
222  histArr = h->DrawSame("pM-wghtWeight",
223  "nomc goff"); //NOTE w/o "can" it crashes
224  // if(!histArr && h) histArr = h->DrawSame("pM","can nomc goff"); //NOTE w/o "can" it crashes
225  // histArr->Print();
226  // get mc input histograms TODO: think about integration
227  if (!histArrMC && hMC)
228  histArrMC =
229  hMC->DrawSame("pM", "onlymc eff goff"); //NOTE w/o "can" it crashes
230 
231 
232  // process raw signal extraction
233  sig->Process(histArr);
234  if (gErrorIgnoreLevel < kWarning) sig->Print("");
235 
236  // validate signal extraction
237  if (sig->GetSignal() < 0.) continue;
238 
239  // fill the tree
240  fExt->setup = fInputKeys[i];
241  fExt->setupId = i;
242  fExt->poi = sig->GetParticleOfInterest();
243  fExt->var = 1.; //integrated
244  fExt->varE = 0.;
245  fExt->s = sig->GetSignal();
246  fExt->sE = sig->GetSignalError();
247  fExt->b = sig->GetBackground();
248  fExt->bE = sig->GetBackgroundError();
249  fExt->sb = sig->GetSB();
250  fExt->sbE = sig->GetSBError();
251  fExt->sgn = sig->GetSignificance();
252  fExt->sgnE = sig->GetSignificanceError();
253  fExt->HistSignal =
254  NULL; //dynamic_cast<TH1F*>(sig->GetSignalHistogram());
255  fExt->eff = 1.; //TODO: calculate this
256  fExt->effE = 0.; //TODO: calculate this
257  fExt->signal = sig; //NULL;
258  fExt->sref = 1.; //TODO: calculate this
259  fExt->srefE = 0.; //TODO: calculate this
260 
261  fTree->Fill();
262  } else {
263 
265  if (h && !h->SetCutClass(fInputKeys[i].Data())) continue;
266 
267  if (!histArr && h)
268  histArr = h->DrawSame(Form("pM_%s", fVar.Data()), "nomc goff");
269  TH2* histPM = (TH2*) sig->FindObject(histArr, PairAnalysis::kSEPM);
270  if (!histPM) return;
271 
272  TObjArray tmpArr;
273  tmpArr.SetOwner(kFALSE);
274 
276  TH1* histMC = NULL;
277  if (!histArrMC && hMC) {
278  histArrMC = hMC->DrawSame(Form("p%s", fVar.Data()),
279  "goff sel",
280  "phi"); // TODO: add search using fPOIpdg
281  }
282  if (histArrMC) {
283  TH1* tmpMCnom = (TH1*) histArrMC->At(0);
284  if (histArrMC->GetEntriesFast() < 2) return;
285  TH1* tmpMCden = (TH1*) histArrMC->At(1);
287  histMC = tmpMCnom->Rebin(
288  fVarBinning->GetNrows() - 1, "effMC", fVarBinning->GetMatrixArray());
289  TH1* histMCden = tmpMCden->Rebin(fVarBinning->GetNrows() - 1,
290  "effMCden",
291  fVarBinning->GetMatrixArray());
292  histMC->Divide(histMCden);
293  delete histMCden;
294  }
296  if (histMC) TObject::Info("Process", "MC histogram found and rebinned");
297 
299  TH1* histMCtruth = NULL;
300  if (histArrMCt) {
301  TH1* tmpMCtrue = (TH1*) histArrMCt->At(0);
302  if (!tmpMCtrue) return;
304  histMCtruth = tmpMCtrue->Rebin(fVarBinning->GetNrows() - 1,
305  "sMCtrue",
306  fVarBinning->GetMatrixArray());
307  }
309  if (histMCtruth)
310  TObject::Info("Process",
311  "MCtruth reference histogram found and rebinned");
312 
313  // loop over all bins
314  for (Int_t bin = 0; bin < fVarBinning->GetNrows() - 1; bin++) {
315  tmpArr.Clear();
316 
317  // var bin limits
318  Double_t xLo = fVarBinning->GetMatrixArray()[bin];
319  Double_t xHi = fVarBinning->GetMatrixArray()[bin + 1] - 0.000001;
320  // axis limits
321  Int_t binLo = histPM->GetYaxis()->FindBin(xLo);
322  Int_t binHi = histPM->GetYaxis()->FindBin(xHi);
323  // found bin limits
324  Double_t fndLo = histPM->GetYaxis()->GetBinLowEdge(binLo);
325  Double_t fndHi = histPM->GetYaxis()->GetBinLowEdge(binHi + 1);
326  //printf("binning requested, found of %s: %.3f-%.3f, %.3f-%.3f \n", fVar.Data(), xLo,xHi, fndLo,fndHi );
327  TObject::Info("Process",
328  "Bin %d: %.3f < %s < %.3f",
329  bin,
330  fndLo,
331  fVar.Data(),
332  fndHi);
333 
334  // fill array for signal extraction
335  for (Int_t ih = 0; ih < histArr->GetEntriesFast(); ih++) {
336  TH2* hist = (TH2*) histArr->UncheckedAt(ih);
337  tmpArr.Add(hist->ProjectionX(hist->GetTitle(), binLo, binHi, "e"));
338  }
339  // tmpArr.Print();
340 
341  // process
342  sig->Process(&tmpArr);
343  if (gErrorIgnoreLevel < kWarning) sig->Print("");
344 
345  // validate signal extraction
346  if (sig->GetSignal() < 0.) continue;
347  if (TMath::IsNaN(sig->GetSignalError())) continue;
348 
349  // fill the tree
350  fExt->setup = fInputKeys[i];
351  fExt->setupId = i;
352  fExt->poi = sig->GetParticleOfInterest();
353  fExt->var = (fndHi - fndLo) / 2 + fndLo; // center of the found bin
354  fExt->varE = (fndHi - fndLo) / 2; // found bin width
355  fExt->s = sig->GetSignal();
356  fExt->sE = sig->GetSignalError();
357  fExt->b = sig->GetBackground();
358  fExt->bE = sig->GetBackgroundError();
359  fExt->sb = sig->GetSB();
360  fExt->sbE = sig->GetSBError();
361  fExt->sgn = sig->GetSignificance();
362  fExt->sgnE = sig->GetSignificanceError();
363  fExt->HistSignal =
364  NULL; //dynamic_cast<TH1F*>(sig->GetSignalHistogram());
365  fExt->eff = (histMC ? histMC->GetBinContent(bin + 1) : 1.);
366  fExt->effE = (histMC ? histMC->GetBinError(bin + 1) : 0.);
367  fExt->signal = sig;
368  fExt->sref = (histMCtruth ? histMCtruth->GetBinContent(bin + 1) : 0.);
369  fExt->srefE = (histMCtruth ? histMCtruth->GetBinError(bin + 1) : 0.);
370 
371  fTree->Fill();
372 
373  // set variable
374  // hf->AddCutVariable( (EValueTypes)
375  // PairAnalysisVarManager::GetValueType(fVar.Data()),
376  // fVarBinning->GetMatrixArray()[bin],
377  // fVarBinning->GetMatrixArray()[bin+1]
378  // );
379 
380  // get histogram array
381  //...
382 
383  } //end binning
384 
386  if (histMC) delete histMC;
387  if (histMCtruth) delete histMCtruth;
388 
389  } // end 2D
390 
391  } // end raw content
392 
393  // fTree->Print();
394 }
395 
396 //______________________________________________
397 void PairAnalysisSpectrum::DrawSpectrum(const char* varexp,
398  const char* selection,
399  Option_t* option) {
400  //
401  // TTree draw alias
402  //
415 
416  TString optString(option);
417  optString.ToLower();
418  printf("Plot spectrum: '%s' \t selection: '%s' \t options: '%s' \n",
419  varexp,
420  selection,
421  optString.Data());
422  Bool_t optLegFull = optString.Contains("legf");
423  optString.ReplaceAll("legf", "");
424  Bool_t optLeg = optString.Contains("leg");
425  optString.ReplaceAll("leg", "");
426  Bool_t optSyst = optString.Contains("syst");
427  optString.ReplaceAll("syst", "");
428  Bool_t optPrint = optString.Contains("print");
429  optString.ReplaceAll("print", "");
430  Bool_t optSamePad = optString.Contains("samepad");
431  optString.ReplaceAll("samepad", "");
432 
434  Long64_t n = 1;
435 
438 
440  TString ckey(varexp);
441  Int_t ndim = ckey.CountChar(':') + 1;
442 
444  TObjArray* carr = ckey.Tokenize(":");
445  carr->SetOwner();
446  // delete carr;
447 
449  TObjArray* oarr = ckey.Tokenize("+-*/:");
450  oarr->SetOwner();
451  TString fkey = ((TObjString*) oarr->At(0))->GetString();
452  delete oarr;
453 
455  ckey.ReplaceAll("/", "#");
456  TCanvas* c = NULL;
457  if (optSamePad)
458  c = gPad->GetCanvas();
459  else
460  c = (TCanvas*) gROOT->FindObject(ckey.Data());
461  if (!c) {
462  TObject::Info("Draw", "create new canvas: '%s'", ckey.Data());
463  c = new TCanvas(ckey.Data(), ckey.Data());
464  }
465  c->cd();
466 
468  TObject* obj;
469  Int_t nobj = 0;
470  TList* prim = gPad->GetListOfPrimitives();
472  for (Int_t io = 0; io < prim->GetSize(); io++) {
473  obj = prim->At(io);
474  if (obj->InheritsFrom(TGraph::Class()) && obj != prim->At(io + 1)) nobj++;
475  }
476 
477  TLegend* leg = 0;
478  if ((optLeg && !nobj)) {
479  // if ( (optLeg && optTask && !nobj) || (optLeg && !optTask && !optDet) ) {
480  leg = new TLegend(0. + gPad->GetLeftMargin() + gStyle->GetTickLength("Y"),
481  0. + gPad->GetBottomMargin() + gStyle->GetTickLength("X"),
482  1. - gPad->GetRightMargin() + gStyle->GetTickLength("Y"),
483  1. - gPad->GetTopMargin() + gStyle->GetTickLength("X"),
484  GetName(),
485  "nbNDC");
486  } else if (optLeg && nobj) {
487  leg = (TLegend*) prim->FindObject("TPave");
488  }
489 
490  Info("Draw", "Basics: nobj: %d \t leg: %p", nobj, leg);
491 
492  // logaritmic style
493  if (optString.Contains("logx")) gPad->SetLogx();
494  if (optString.Contains("logy")) gPad->SetLogy();
495  if (optString.Contains("logz")) gPad->SetLogz();
496  optString.ReplaceAll("logx", "");
497  optString.ReplaceAll("logy", "");
498  optString.ReplaceAll("logz", "");
499 
500  if (ndim < 3
501  && !ckey.Contains("signal->")) { // build own histogram with labels
502 
503  // printf("try to get errors for %d dimensional input \n",ndim);
504  TString varkey(varexp);
505 
506  // get event list
507  n = fTree->Draw(">>elist", selection);
508  if (!n) {
509  delete carr;
510  return;
511  }
512  TEventList* elist = (TEventList*) gDirectory->Get("elist");
513  if (elist) {
514  elist->SetReapplyCut(kTRUE); // important!
515  // elist->Print("all");
516  fTree->SetEventList(elist);
517  }
518 
519  TH1D* hist = 0x0;
520  // set up a proper histogram with labels
521  if (varkey.Contains("setup")) {
522 
523  hist = new TH1D(); // activate buffer
524  hist->SetNameTitle(varexp, selection);
525 
526  //read strings for all entries
527  for (Long64_t i = 0; i < n; i++) {
528  fTree->GetEntry(i);
529  hist->Fill((fExt->setup).Data(), 1.);
530  // printf(" read from tree: %s \n",(fExt->setup).Data());
531  }
532 
533  hist->Draw("AXIS");
534  hist->GetXaxis()->SetRange(1, n);
535  }
536 
537  // get errors from tree
538  TString errkey = ((TObjString*) carr->At(0))->GetString();
539  errkey.ReplaceAll(fkey, fkey + "E");
540  // if(ndim>1) varkey.ReplaceAll(fkey+":",fkey+":"+fkey+"E:");
541  // else varkey.ReplaceAll(fkey,fkey+":"+fkey+"E");
542  // if(ndim>1) varkey.ReplaceAll(fkey+":",fkey+":"+fkey+"E:");
543  // else
544  if (!fkey.Contains("E")) {
545  varkey.Append(errkey.Prepend(":"));
546  Info("Draw",
547  " Appended '%s' by '%s' for error caluclation",
548  varkey.Data(),
549  errkey.Data());
550  }
551 
553  //printf("execute collect/draw command for %s \n",varkey.Data());
554  fTree->Draw(varkey, selection, "goff");
555 
557  Double_t* xval = new Double_t[fTree->GetSelectedRows()];
558  for (Int_t ix = 0; ix < fTree->GetSelectedRows(); ix++)
559  xval[ix] = 1.;
560 
562  TGraphErrors* gr = 0x0;
563  if (ndim > 1)
564  gr = new TGraphErrors(fTree->GetSelectedRows(),
565  fTree->GetV2(),
566  fTree->GetV1(),
567  0,
568  fTree->GetV3());
569  else
570  gr = new TGraphErrors(
571  fTree->GetSelectedRows(), fTree->GetV1(), xval, fTree->GetV2(), 0);
572  delete[] xval;
573 
574  if (!gr) {
575  delete carr;
576  return;
577  }
578 
581  TString sel(selection);
582  sel.ReplaceAll("setup.Contains", "");
583  sel.ReplaceAll("(\"", "*");
584  sel.ReplaceAll("\")", "*");
585  sel.ReplaceAll("setup==", "");
586  sel.ReplaceAll("\"", "");
587  gr->SetName(Form("%s", sel.Data()));
588  if (sel.Contains("setupId==") && sel.Length() < 11) {
589  sel.ReplaceAll("setupId==", "");
590  Int_t iId = sel.Atoi();
591  gr->SetName(Form("%s", fInputKeys[iId].Data()));
592  }
593  if (optSyst) gr->SetName(GetTitle());
594  if (fkey.Contains("ref")) gr->SetName("MC-truth");
595 
596  // sort x-values
597  gr->Sort();
598  TGraphErrors* grE = NULL; // statistical
599  TGraphErrors* grS = NULL; // systematic
600  TGraphErrors* grC = NULL; // stat + syst
601  if (optSyst
602  && fVarBinning) { //TODO: implement systematic calculation w/o binning
603  grE = new TGraphErrors(); // statistical graph
604  grS = new TGraphErrors(); // systematics graph
605  grC = new TGraphErrors(); // systematics graph
606  Double_t* gx = gr->GetX();
607  Double_t* gy = gr->GetY();
608  Double_t* gye = gr->GetEY();
609 
610  // loop over all variable bins
611  Int_t first = 0; //TMath::BinarySearch(gr->GetN(),gx,xLo); //first bin
612  Int_t ibin = 0;
613  while (first < gr->GetN()) {
614  // for(ibin=0; ibin<fVarBinning->GetNrows()-1; ibin++) {
615 
616  Int_t nsys = 0; // counter #systematics
617  Double_t ysys = 0.; // y-vlaue of systematic = mean value
618  Double_t esys = 0.; // uncertainty depends on method
619  // calculate mean
620  Double_t xvar = gx[first];
621  for (Int_t i = first; i < gr->GetN(); i++) {
622  // printf("gx[i]:%f xvar:%f xLo:%f\n",gx[i],xvar,xLo);
623  // if(TMath::Abs(gx[i]-xvar)>1.e-8) break;
624  if ((gx[i] - xvar) > 1.e-8) break;
625  // printf("graph entry %d, found index first: %d with value %.3f \t y: %.3f+-%.3f \n",
626  // i,first,xvar,gy[i],gye[i]);
627  nsys++;
628  ysys += gy[i];
629  }
630  ysys /= (nsys ? nsys : 1); // protect for zero division
631  // printf("bin %d, found index first: %d, %.3f \t y:%3.f \t nsys:%d\n",ibin,first,xvar,ysys,nsys);
632 
633  // y-syst. uncertainty
634  // printf("syst %f , <y> %f\n",esys,ysys);
635  for (Int_t i = 0; i < nsys; i++) {
636  Int_t j = first + i;
637  // if(gx[j]!=xLo) break; // check should not be needed
638  Double_t uce = 0.;
639  switch (fSystMthd) {
640  case kBarlow:
641  // I. calc uncorr. stat. error from sub/superset w.r.t. first measurement
642  uce = TMath::Sqrt(
643  TMath::Abs(gye[j] * gye[j] - gye[first] * gye[first]));
644  // II. calc max. deviation w.r.t. mean y-value incl. 1sigma* 0.9 of I.
645  // NOTE: 0.9 can be change to a max value of 1->1sigma, 0.9 is more consevative
646  esys = TMath::Max(esys, TMath::Abs(ysys - gy[j]) - 0.9 * uce);
647  break;
648  case kSystMax:
649  esys = TMath::Max(esys, TMath::Abs(ysys - gy[j]));
650  break;
651  case kSystRMS: esys += gy[j] * gy[j]; break;
652  }
653  // printf("bin error %f \t syst %f from abs %f \n",gye[j],esys, TMath::Abs( gy[j] ));
654  }
655 
656  // normalisation
657  switch (fSystMthd) {
658  case kBarlow: /* nothing to be done */ break;
659  case kSystMax: /* nothing to be done */ break;
660  case kSystRMS:
661  esys =
662  TMath::Sqrt(TMath::Abs(esys / (nsys ? nsys : 1) - ysys * ysys));
663  break;
664  }
665 
666  // fill statistical and systematic graph values and errors
667  grE->SetPoint(ibin, xvar, ysys); // mean
668  grE->SetPointError(ibin, 0.0, gye[first]); // stat.uncert. of first set
669 
670  Double_t boxW = (fVarBinning->Max() - fVarBinning->Min())
671  / (fVarBinning->GetNrows() - 1);
672  grS->SetPoint(ibin, xvar, ysys); // mean
673  grS->SetPointError(ibin, boxW * 0.35, esys); // systematic value
674 
675  // calculate err = sqrt(stat.**2 + syst**2)
676  grC->SetPoint(ibin, xvar, ysys); // mean
677  grC->SetPointError(ibin,
678  boxW * 0.35,
679  TMath::Sqrt(esys * esys + gye[first] * gye[first]));
680 
681  // increase index counter
682  first += nsys;
683  ibin++;
684  } //next bin
685  // grS->Print();
686  } else {
687  grC = new TGraphErrors(*gr);
688  }
689 
690  if (optPrint) grC->Print();
691 
692  Info("Draw", " Draw object with options: '%s'", optString.Data());
694  gr->Draw((optString + "A").Data());
695  else {
696  gr->Draw((optString + "same").Data());
697 
698  // set axis maximum
699  Double_t* valE = grC->GetEY();
700  Double_t* val = grC->GetY();
701  Int_t npnts = grC->GetN();
702  Int_t idx[1000];
703  TMath::Sort(npnts, val, idx, kTRUE); // kFALSE=increasing numbers
704 
705  Double_t errmin =
706  (TMath::IsNaN(valE[idx[npnts - 1]]) ? 0. : valE[idx[npnts - 1]]);
707  Double_t min = (val[idx[npnts - 1]] - errmin) * 0.9;
708  Double_t tmpmin = PairAnalysisStyler::GetFirstHistogram()->GetMinimum();
710  (tmpmin < min ? tmpmin : min));
711  Double_t errmax = (TMath::IsNaN(valE[idx[0]]) ? 0. : valE[idx[0]]);
712  Double_t max = (val[idx[0]] + errmax) * 1.1;
713  Double_t tmpmax = PairAnalysisStyler::GetFirstHistogram()->GetMaximum();
715  (tmpmax > max ? tmpmax : max));
716  }
717 
718  // draw systemtaic graph ontop
719  if (grS) {
720  PairAnalysisStyler::Style(grE, nobj);
721  grE->Draw((optString + "A").Data());
722  PairAnalysisStyler::Style(grS, nobj);
723  grS->SetFillColor(grS->GetLineColor());
724  grS->SetFillStyle(kFEmpty);
725  grS->Draw("2same");
726  }
727 
728  // legend
729  TString legOpt = optString + "L";
730  legOpt.ReplaceAll("hist", "");
731  legOpt.ReplaceAll("scat", "");
732  if (legOpt.Contains("col")) legOpt = "";
733  legOpt.ReplaceAll("z", "");
734  legOpt.ReplaceAll("e", "");
735 
736  TString legkey = gr->GetName();
737  if (leg) leg->AddEntry(gr, gr->GetName(), legOpt.Data());
738 
739  fSignal = grC;
740  PairAnalysisStyler::Style(fSignal, nobj);
741 
742  } else {
743  // execute tree draw command
744  fTree->Draw(varexp, selection, optString.Data());
745  fprintf(stderr, "use plain TTree::Draw command \n");
746  return;
747  }
748 
749 
751  // printf("modify axis titles \n");
752  if (!optSamePad) {
753  UInt_t varx = PairAnalysisVarManager::GetValueType(fVar.Data());
754  TString var(varexp);
755  TObjArray* arr = var.Tokenize(":");
756  arr->SetOwner();
757  TString xt = "";
758  TString yt = "Entries";
759  xt = ((TObjString*) arr->At(0))->GetString();
760  if (xt.EqualTo("sb"))
762  else if (xt.EqualTo("s"))
764  else if (xt.EqualTo("b"))
766  else if (xt.EqualTo("sgn"))
768  else if (xt.EqualTo("var"))
769  xt = Form("%s %s",
772  else if (xt.Contains("var")) {
773  xt.ReplaceAll(
774  "varE", Form("#Delta%s", PairAnalysisVarManager::GetValueLabel(varx)));
775  xt.ReplaceAll("var", PairAnalysisVarManager::GetValueLabel(varx));
776  }
777 
778  if (arr->GetEntriesFast() < 2) {
779  PairAnalysisStyler::GetFirstHistogram()->SetXTitle(xt.Data());
780  PairAnalysisStyler::GetFirstHistogram()->SetYTitle(yt.Data());
781  } else {
782  PairAnalysisStyler::GetFirstHistogram()->SetYTitle(xt.Data());
783  xt = ((TObjString*) arr->At(1))->GetString();
784  if (xt.EqualTo("sb"))
786  else if (xt.EqualTo("s"))
788  else if (xt.EqualTo("b"))
790  else if (xt.EqualTo("sgn"))
792  else if (xt.EqualTo("var"))
793  xt = Form("%s %s",
796  else if (xt.Contains("var")) {
797  xt.ReplaceAll(
798  "varE",
799  Form("#Delta%s", PairAnalysisVarManager::GetValueLabel(varx)));
800  xt.ReplaceAll("var", PairAnalysisVarManager::GetValueLabel(varx));
801  }
802  PairAnalysisStyler::GetFirstHistogram()->SetXTitle(xt.Data());
803  }
804 
806  if (arr) delete arr;
807  }
808 
810  if (carr) delete carr;
811 
813  if (fVarBinning)
815  fVarBinning->Min(), fVarBinning->Max(), "X");
816  // PairAnalysisStyler::GetFirstHistogram()->GetXaxis()->SetNdivisions(, 0, 0, kFALSE);
817 
818 
819  // legend
820  if (leg) {
822  leg, optLegFull); // coordinates, margins, fillstyle, fontsize
823  if (!nobj) leg->Draw(); // only draw the legend once
825  }
826 
828  fTree->SetEventList(0);
829 }
830 
831 //______________________________________________
832 void PairAnalysisSpectrum::Write() {
836  TFile* fout = new TFile("test.root", "RECREATE");
837  fout->cd();
838  // fTree->Print();
839  fTree->Write();
840  fout->Close();
841 }
842 
843 //______________________________________________
844 void PairAnalysisSpectrum::Fit(TString drawoption) {
852 
853  drawoption.ToLower();
854  Bool_t optLeg = drawoption.Contains("leg");
855  drawoption.ReplaceAll("leg", "");
856 
857  // fSignal->Print();
858 
859  Info("Fit", "Spectrum fit method: %s", fFuncSigBack->GetName());
860  Int_t fitResult = fSignal->Fit(fFuncSigBack, (fFitOpt + "EX0").Data());
861 
862  // warning in case of fit issues
863  if (fitResult != 0) {
864  Error("Fit", "fit has error/issue (%d)", fitResult);
865  return;
866  }
867 
869  fFuncSigBack->SetLineColor(fSignal->GetLineColor());
870  // fFuncSigBack->SetLineStyle(kDashed);
871 
872  // PairAnalysisStyler::Style(fit, PairAnalysisStyler::kFit);
873  TF1* fit = fFuncSigBack->DrawCopy((drawoption + "same").Data());
874 
876  fDof = fFuncSigBack->GetNDF();
877  if (fDof) fChi2Dof = fFuncSigBack->GetChisquare() / fFuncSigBack->GetNDF();
878 
880  if (optLeg) {
881  TList* prim = gPad->GetListOfPrimitives();
882  TLegend* leg = (TLegend*) prim->FindObject("TPave");
883  TString legkey = fFuncSigBack->GetName();
885  if (leg) {
886  leg->AddEntry(fit, fFuncSigBack->GetName(), drawoption.Data());
889  }
890  }
891 }
PairAnalysisSignalExt::GetSignificance
Double_t GetSignificance() const
Definition: PairAnalysisSignalExt.h:130
PairAnalysisSpectrum.h
PairAnalysisFunction
Definition: PairAnalysisFunction.h:22
PairAnalysisStyler::SetLegendAttributes
void SetLegendAttributes(TLegend *leg, Bool_t fill=kFALSE)
Definition: PairAnalysisStyler.cxx:466
PairAnalysisSignalExt::Print
void Print(Option_t *option="") const
Definition: PairAnalysisSignalExt.cxx:214
PairAnalysisSignalExt::GetValueName
static const char * GetValueName(Int_t i)
Definition: PairAnalysisSignalExt.h:140
PairAnalysisStyler::kFit
@ kFit
Definition: PairAnalysisStyler.h:25
PairAnalysisSignalExt
Definition: PairAnalysisSignalExt.h:25
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
gr
Char_t * gr
Definition: CbmEvDisTracks.cxx:289
PairAnalysisVarManager::GetValueType
static UInt_t GetValueType(const char *valname)
Definition: PairAnalysisVarManager.cxx:355
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition: L1/vectors/P4_F32vec4.h:33
PairAnalysisSignalExt.h
PairAnalysisVarManager.h
PairAnalysisSignalExt::GetSBError
Double_t GetSBError() const
Definition: PairAnalysisSignalExt.h:133
h
Data class with information on a STS local track.
PairAnalysisFunction::GetParticleOfInterest
Int_t GetParticleOfInterest() const
Definition: PairAnalysisFunction.h:85
PairAnalysisVarManager::GetValueLabel
static const char * GetValueLabel(Int_t i)
Definition: PairAnalysisVarManager.h:380
PairAnalysisSignalExt::GetBackground
Double_t GetBackground() const
Definition: PairAnalysisSignalExt.h:128
PairAnalysisSignalExt::Process
void Process(TObjArray *const arrhist)
Definition: PairAnalysisSignalExt.cxx:596
PairAnalysisHF
Definition: PairAnalysisHF.h:24
PairAnalysis::kSEPM
@ kSEPM
Definition: PairAnalysis.h:28
Extraction
Definition: PairAnalysisSpectrum.h:32
first
bool first
Definition: LKFMinuit.cxx:143
PairAnalysisSignalExt::FindObject
TObject * FindObject(TObjArray *arrhist, PairAnalysis::EPairType type) const
Definition: PairAnalysisSignalExt.h:283
PairAnalysisSignalExt::GetBackgroundError
Double_t GetBackgroundError() const
Definition: PairAnalysisSignalExt.h:129
PairAnalysisStyler::LoadStyle
void LoadStyle()
Definition: PairAnalysisStyler.cxx:49
PairAnalysisStyler.h
PairAnalysisHF.h
PairAnalysisStyler::GetFirstHistogram
TH1 * GetFirstHistogram()
Definition: PairAnalysisStyler.cxx:605
PairAnalysisSignalExt::GetSB
Double_t GetSB() const
Definition: PairAnalysisSignalExt.h:132
PairAnalysisSignalExt::GetSignificanceError
Double_t GetSignificanceError() const
Definition: PairAnalysisSignalExt.h:131
ClassImp
ClassImp(PairAnalysisSpectrum) PairAnalysisSpectrum
Definition: PairAnalysisSpectrum.cxx:77
PairAnalysisSignalExt::GetSignal
Double_t GetSignal() const
Definition: PairAnalysisSignalExt.h:126
PairAnalysisStyler::Style
void Style(TObject *obj, Int_t idx=0)
Definition: PairAnalysisStyler.cxx:262
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: L1/vectors/P4_F32vec4.h:36
PairAnalysisHistos.h
PairAnalysisVarManager::GetValueUnit
static const char * GetValueUnit(Int_t i)
Definition: PairAnalysisVarManager.h:383
PairAnalysisSignalExt::GetSignalError
Double_t GetSignalError() const
Definition: PairAnalysisSignalExt.h:127