48 #include <TEventList.h>
55 #include <TDatabasePDG.h>
59 #include <TGraphErrors.h>
62 #include <TPaveText.h>
64 #include <TProfile2D.h>
80 PairAnalysisSpectrum::PairAnalysisSpectrum()
81 : PairAnalysisSpectrum(
"spectrum",
"title") {
88 PairAnalysisSpectrum::PairAnalysisSpectrum(
const char* name,
const char* title)
97 fTree =
new TTree(
"PAPa",
"PAPa-Spectrum");
99 for (Int_t
i = 0;
i < 100;
i++)
104 PairAnalysisSpectrum::~PairAnalysisSpectrum() {
108 if (fResults)
delete fResults;
109 if (fExtractions)
delete fExtractions;
113 void PairAnalysisSpectrum::AddInput(TObjArray* raw,
122 if (mc) fMCInput.Add(mc);
123 if (truth) fMCTruth.Add(truth);
124 fInputKeys[fIdx] = identifier;
129 void PairAnalysisSpectrum::Init() {
133 fTree->Branch(
"Extraction", &fExt);
138 void PairAnalysisSpectrum::Process() {
146 TIter nextRaw(&fRawInput);
149 PairAnalysisHistos*
h = NULL;
151 TObject* objMC = NULL;
153 PairAnalysisHistos* hMC = NULL;
155 TObject* objMCtruth = NULL;
158 while ((obj = nextRaw())) {
160 TObject::Info(
"Process",
"Check Extraction for %s", fInputKeys[
i].Data());
161 TObject::Info(
"Process",
162 "------------------------------------------------------------"
164 TObject::Info(
"Process",
"Input type: %s \n", obj->ClassName());
167 TObject::Info(
"Process",
168 "Binning provided for %s from %.3f to %.3f",
175 TObjArray* histArr =
dynamic_cast<TObjArray*
>(obj);
177 if (!(
h =
dynamic_cast<PairAnalysisHistos*
>(obj))) {
179 TObject::Error(
"Process",
"No input format found");
192 objMC = fMCInput.At(
i);
193 TObjArray* histArrMC = NULL;
195 TObject::Info(
"Process",
"Input MC type: %s \n", objMC->ClassName());
196 histArrMC =
dynamic_cast<TObjArray*
>(objMC);
198 if (!(hMC =
dynamic_cast<PairAnalysisHistos*
>(objMC))) {
200 TObject::Error(
"Process",
"No MC input format found");
208 objMCtruth = fMCTruth.At(
i);
209 TObjArray* histArrMCt = NULL;
212 "Process",
"Input MC truth type: %s \n", objMCtruth->ClassName());
213 histArrMCt =
dynamic_cast<TObjArray*
>(objMCtruth);
222 histArr =
h->DrawSame(
"pM-wghtWeight",
227 if (!histArrMC && hMC)
229 hMC->DrawSame(
"pM",
"onlymc eff goff");
234 if (gErrorIgnoreLevel < kWarning) sig->
Print(
"");
240 fExt->setup = fInputKeys[
i];
249 fExt->sb = sig->
GetSB();
265 if (
h && !
h->SetCutClass(fInputKeys[
i].Data()))
continue;
268 histArr =
h->DrawSame(Form(
"pM_%s", fVar.Data()),
"nomc goff");
273 tmpArr.SetOwner(kFALSE);
277 if (!histArrMC && hMC) {
278 histArrMC = hMC->DrawSame(Form(
"p%s", fVar.Data()),
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,
291 fVarBinning->GetMatrixArray());
292 histMC->Divide(histMCden);
296 if (histMC) TObject::Info(
"Process",
"MC histogram found and rebinned");
299 TH1* histMCtruth = NULL;
301 TH1* tmpMCtrue = (TH1*) histArrMCt->At(0);
302 if (!tmpMCtrue)
return;
304 histMCtruth = tmpMCtrue->Rebin(fVarBinning->GetNrows() - 1,
306 fVarBinning->GetMatrixArray());
310 TObject::Info(
"Process",
311 "MCtruth reference histogram found and rebinned");
314 for (Int_t bin = 0; bin < fVarBinning->GetNrows() - 1; bin++) {
318 Double_t xLo = fVarBinning->GetMatrixArray()[bin];
319 Double_t xHi = fVarBinning->GetMatrixArray()[bin + 1] - 0.000001;
321 Int_t binLo = histPM->GetYaxis()->FindBin(xLo);
322 Int_t binHi = histPM->GetYaxis()->FindBin(xHi);
324 Double_t fndLo = histPM->GetYaxis()->GetBinLowEdge(binLo);
325 Double_t fndHi = histPM->GetYaxis()->GetBinLowEdge(binHi + 1);
327 TObject::Info(
"Process",
328 "Bin %d: %.3f < %s < %.3f",
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"));
343 if (gErrorIgnoreLevel < kWarning) sig->
Print(
"");
350 fExt->setup = fInputKeys[
i];
353 fExt->var = (fndHi - fndLo) / 2 + fndLo;
354 fExt->varE = (fndHi - fndLo) / 2;
359 fExt->sb = sig->
GetSB();
365 fExt->eff = (histMC ? histMC->GetBinContent(bin + 1) : 1.);
366 fExt->effE = (histMC ? histMC->GetBinError(bin + 1) : 0.);
368 fExt->sref = (histMCtruth ? histMCtruth->GetBinContent(bin + 1) : 0.);
369 fExt->srefE = (histMCtruth ? histMCtruth->GetBinError(bin + 1) : 0.);
386 if (histMC)
delete histMC;
387 if (histMCtruth)
delete histMCtruth;
397 void PairAnalysisSpectrum::DrawSpectrum(
const char* varexp,
398 const char* selection,
416 TString optString(option);
418 printf(
"Plot spectrum: '%s' \t selection: '%s' \t options: '%s' \n",
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",
"");
440 TString ckey(varexp);
441 Int_t ndim = ckey.CountChar(
':') + 1;
444 TObjArray* carr = ckey.Tokenize(
":");
449 TObjArray* oarr = ckey.Tokenize(
"+-*/:");
451 TString fkey = ((TObjString*) oarr->At(0))->GetString();
455 ckey.ReplaceAll(
"/",
"#");
458 c = gPad->GetCanvas();
460 c = (TCanvas*) gROOT->FindObject(ckey.Data());
462 TObject::Info(
"Draw",
"create new canvas: '%s'", ckey.Data());
463 c =
new TCanvas(ckey.Data(), ckey.Data());
470 TList* prim = gPad->GetListOfPrimitives();
472 for (Int_t io = 0; io < prim->GetSize(); io++) {
474 if (obj->InheritsFrom(TGraph::Class()) && obj != prim->At(io + 1)) nobj++;
478 if ((optLeg && !nobj)) {
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"),
486 }
else if (optLeg && nobj) {
487 leg = (TLegend*) prim->FindObject(
"TPave");
490 Info(
"Draw",
"Basics: nobj: %d \t leg: %p", nobj, leg);
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",
"");
501 && !ckey.Contains(
"signal->")) {
504 TString varkey(varexp);
507 n = fTree->Draw(
">>elist", selection);
512 TEventList* elist = (TEventList*) gDirectory->Get(
"elist");
514 elist->SetReapplyCut(kTRUE);
516 fTree->SetEventList(elist);
521 if (varkey.Contains(
"setup")) {
524 hist->SetNameTitle(varexp, selection);
527 for (Long64_t
i = 0;
i < n;
i++) {
529 hist->Fill((fExt->setup).Data(), 1.);
534 hist->GetXaxis()->SetRange(1, n);
538 TString errkey = ((TObjString*) carr->At(0))->GetString();
539 errkey.ReplaceAll(fkey, fkey +
"E");
544 if (!fkey.Contains(
"E")) {
545 varkey.Append(errkey.Prepend(
":"));
547 " Appended '%s' by '%s' for error caluclation",
554 fTree->Draw(varkey, selection,
"goff");
557 Double_t* xval =
new Double_t[fTree->GetSelectedRows()];
558 for (Int_t ix = 0; ix < fTree->GetSelectedRows(); ix++)
562 TGraphErrors*
gr = 0x0;
564 gr =
new TGraphErrors(fTree->GetSelectedRows(),
570 gr =
new TGraphErrors(
571 fTree->GetSelectedRows(), fTree->GetV1(), xval, fTree->GetV2(), 0);
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()));
593 if (optSyst)
gr->SetName(GetTitle());
594 if (fkey.Contains(
"ref"))
gr->SetName(
"MC-truth");
598 TGraphErrors* grE = NULL;
599 TGraphErrors* grS = NULL;
600 TGraphErrors* grC = NULL;
603 grE =
new TGraphErrors();
604 grS =
new TGraphErrors();
605 grC =
new TGraphErrors();
606 Double_t* gx =
gr->GetX();
607 Double_t* gy =
gr->GetY();
608 Double_t* gye =
gr->GetEY();
613 while (first < gr->GetN()) {
620 Double_t xvar = gx[
first];
624 if ((gx[
i] - xvar) > 1.e-8)
break;
630 ysys /= (nsys ? nsys : 1);
635 for (Int_t
i = 0;
i < nsys;
i++) {
643 TMath::Abs(gye[j] * gye[j] - gye[
first] * gye[
first]));
646 esys = TMath::Max(esys, TMath::Abs(ysys - gy[j]) - 0.9 * uce);
649 esys = TMath::Max(esys, TMath::Abs(ysys - gy[j]));
651 case kSystRMS: esys += gy[j] * gy[j];
break;
659 case kSystMax:
break;
662 TMath::Sqrt(TMath::Abs(esys / (nsys ? nsys : 1) - ysys * ysys));
667 grE->SetPoint(ibin, xvar, ysys);
668 grE->SetPointError(ibin, 0.0, gye[
first]);
670 Double_t boxW = (fVarBinning->Max() - fVarBinning->Min())
671 / (fVarBinning->GetNrows() - 1);
672 grS->SetPoint(ibin, xvar, ysys);
673 grS->SetPointError(ibin, boxW * 0.35, esys);
676 grC->SetPoint(ibin, xvar, ysys);
677 grC->SetPointError(ibin,
679 TMath::Sqrt(esys * esys + gye[
first] * gye[
first]));
687 grC =
new TGraphErrors(*
gr);
690 if (optPrint) grC->Print();
692 Info(
"Draw",
" Draw object with options: '%s'", optString.Data());
694 gr->Draw((optString +
"A").Data());
696 gr->Draw((optString +
"same").Data());
699 Double_t* valE = grC->GetEY();
700 Double_t* val = grC->GetY();
701 Int_t npnts = grC->GetN();
703 TMath::Sort(npnts, val, idx, kTRUE);
706 (TMath::IsNaN(valE[idx[npnts - 1]]) ? 0. : valE[idx[npnts - 1]]);
707 Double_t
min = (val[idx[npnts - 1]] - errmin) * 0.9;
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;
715 (tmpmax >
max ? tmpmax :
max));
721 grE->Draw((optString +
"A").Data());
723 grS->SetFillColor(grS->GetLineColor());
724 grS->SetFillStyle(kFEmpty);
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",
"");
736 TString legkey =
gr->GetName();
737 if (leg) leg->AddEntry(
gr,
gr->GetName(), legOpt.Data());
744 fTree->Draw(varexp, selection, optString.Data());
745 fprintf(stderr,
"use plain TTree::Draw command \n");
755 TObjArray* arr = var.Tokenize(
":");
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"))
772 else if (xt.Contains(
"var")) {
778 if (arr->GetEntriesFast() < 2) {
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"))
796 else if (xt.Contains(
"var")) {
810 if (carr)
delete carr;
815 fVarBinning->Min(), fVarBinning->Max(),
"X");
823 if (!nobj) leg->Draw();
828 fTree->SetEventList(0);
832 void PairAnalysisSpectrum::Write() {
836 TFile* fout =
new TFile(
"test.root",
"RECREATE");
844 void PairAnalysisSpectrum::Fit(TString drawoption) {
853 drawoption.ToLower();
854 Bool_t optLeg = drawoption.Contains(
"leg");
855 drawoption.ReplaceAll(
"leg",
"");
859 Info(
"Fit",
"Spectrum fit method: %s", fFuncSigBack->GetName());
860 Int_t fitResult = fSignal->Fit(fFuncSigBack, (fFitOpt +
"EX0").Data());
863 if (fitResult != 0) {
864 Error(
"Fit",
"fit has error/issue (%d)", fitResult);
869 fFuncSigBack->SetLineColor(fSignal->GetLineColor());
873 TF1* fit = fFuncSigBack->DrawCopy((drawoption +
"same").Data());
876 fDof = fFuncSigBack->GetNDF();
877 if (fDof) fChi2Dof = fFuncSigBack->GetChisquare() / fFuncSigBack->GetNDF();
881 TList* prim = gPad->GetListOfPrimitives();
882 TLegend* leg = (TLegend*) prim->FindObject(
"TPave");
883 TString legkey = fFuncSigBack->GetName();
886 leg->AddEntry(fit, fFuncSigBack->GetName(), drawoption.Data());