#if !defined(__CINT__) || defined(__MAKECINT__) #include "TList.h" #include "TFile.h" #include "TStyle.h" #include "TH1F.h" #include "TH2F.h" #include "THnSparse.h" #include "TLegend.h" #include "TSystem.h" #include "TMath.h" #include "TCanvas.h" #include "TLegend.h" #include "TLatex.h" #include "TF1.h" #include "TLine.h" #include "TPaletteAxis.h" // // #endif #include "/home/shahoian/ALICE/mymacros/SaveCanvas.h" const char kHStatName[]="hStat"; double kEps = 1e-6; // double kdPhiBgTailMin = 0.1; // lower limit of dphi tail to use for bg normalization double kdPhiBgTailMax = 0.3; // upper limit of dphi tail to use for bg normalization // double kWDistBgTailMin = 5.; // lower limit of wgh.distance tail to use for bg normalization double kWDistBgTailMax = 25.; // upper limit of wgh.distance tail to use for bg normalization double kdPhiSgCut = 0.06; // cut in dphi-bent used to extract the signal, extracted from stat histo double kWDistSgCut = 1.5; // cut in w.distance used to extract the signal, extracted from stat histo // enum { kNormShapeDist, // normalize bg tails usig weighted distance shape kNormShapeDPhi, // normalize bg tails usig dPhi-bend shape kNormShapes}; enum { kSclWghMean, // normalize bg tails to data using weighted mean of bin-by-bin ratios kSclIntegral, // normalize bg tails to data using integral kSclTypes}; // histograms to be used for bg nomalization for each of NormShapes used const char* kNormShapeH[kNormShapes] = { "EtaDist", // Weighted distance vs Eta "EtaDPhiS" // Dphi-bent vs distance }; const char* figDir = "fig161110"; TString useBgType = "Inj"; Int_t useShapeType = kNormShapeDist; // which distribution to use for bg normalization Bool_t useMCLB = 0;//kFALSE; // use Comb MC Labels as a template for Bg. Int_t useScaleType = kSclIntegral;//kSclWghMean; // which type of tails normalization to use Double_t useEtaCut = 1.; // cut on eta Double_t useZvMin = -5.; // cut on Z vertex Double_t useZvMax = 5.; // cut on Z vertex Int_t useBinGrEta = -1; // for bg fits group eta bins Int_t useBinGrZv = -1; // for bg fits group Zv bins enum {kBitNormPerEvent=BIT(14)}; // bins for saved parameters in the hStat histo enum {kDummyBin, kEvTot, // events read kEvProcData, // events with data mult.object (ESD or reco) kEvProcInj, // events Injected kEvProcRot, // events Rotated kEvProcMix, // events Mixed // kDPhi, // dphi window kDTht, // dtheta window kNStd, // N.standard deviations to keep kPhiShift, // bending shift kThtS2, // is dtheta scaled by 1/sin^2 kPhiOvl, // overlap params kZEtaOvl, // overlap params kNoOvl, // flag that overlap are suppressed // kPhiRot, // rotation phi kInjScl, // injection scaling kEtaCut, // eta cut kZVMin, // min ZVertex to process kZVMax, // max ZVertex to process kTrcMin, // min mult to process kTrcMax, // max mult to process // kOneUnit=49, // just 1 to track mergings kNWorkers=50, // n workers kNStatBins }; enum {kSigCorr,kMCPrim,kRawDtCut,kSignalEst,kSignalEstMC,kBgEst,k1MBeta,k1MBetaMC,kAlpha,kAlphaMC,kBgMC,kBgRescFc,kDataDist,kBgDist,kBgMCDist, kNHistos, kMCShift=20}; void CorrectSpectra(const char* flNameData, const char* flNameMC,const char* unique=""); void PrepareHistos(TList* lst, Bool_t isMC, TObjArray* outArray); void ProcessHistos(TObjArray* outArray); TH1* NormalizeBg(TH1* dataH, TH1* bgH, double &scl, double &scle); TObject* FindObject(const char* nameH, const TList* lst, Bool_t normPerEvent=kTRUE); TList* LoadList(const char* flName, const char* addPref, const char* nameL="clist"); void GetRatE(double x,double xe, double y,double ye, double &rat, double &rate); Int_t CheckStat(const TList* lst,const char* dtType); void Integrate(TH1* hist, double xmn,double xmx, double &val, double& err); void CropHisto(TH1* histo, int b00, int b01, int b10=-1, int b11=-1); void CropHisto(TH1* histo, double b00, double b01, double b10=-1, double b11=-1); void GetRealMinMax(TH1* h, double &vmn, double &vmx); const char* HName(const char* prefix,const char* htype); void ProjDZE(THnSparseF* hrawd, THnSparseF* hgenb, THnSparseF* hmcComb, TObjArray* res, int bstep0=-1,int bstep1=-1); void PlotResults(); void PlotDNDEta(); void PlotDist(); void PlotAlphaBeta(); void PlotSpecies(); TList *listDt=0, *listMC=0; TObjArray resArr; char outStr[1000]; char outTitle[1000]; TString uniqueName=""; void CorrectSpectra(const char* flNameData, const char* flNameMC, const char* uniqueNm) { // uniqueName = uniqueNm; listDt = LoadList(flNameData,"dt_"); listMC = LoadList(flNameMC,"mc_"); // resArr.Clear(); PrepareHistos(listMC, kTRUE, &resArr); PrepareHistos(listDt, kFALSE, &resArr); // ProcessHistos(&resArr); // sprintf(outStr,"CutEta%.1f_Zv%.1f_%.1f_bg%s_Shape_%s_mcLB%d_cutSig%.1f_cutBg%.1f",useEtaCut,useZvMin,useZvMax,useBgType.Data(), useShapeType==kNormShapeDist ? "wdst":"dphi", useMCLB, useShapeType==kNormShapeDist ? kWDistSgCut:kdPhiSgCut, useShapeType==kNormShapeDist ? kWDistBgTailMin:kdPhiBgTailMin); sprintf(outTitle,"%s, |#eta|<%.1f, %.1fGetName(),useShapeType); if (!CheckStat(lst,useBgType.Data())) {printf("Bg of type %s is absent in list %s\n",useBgType.Data(),lst->GetName()); return;} // THnSparseF* hD = (THnSparseF*) FindObject(HName("Data",xxZvEta),lst); THnSparseF* hB = (THnSparseF*) FindObject(HName(useBgType.Data(),xxZvEta),lst); // // special feature: use MC Labels bg as a shape instead of generated bg if (useMCLB/* && !isMC*/) { TString nm = hB->GetName(); nm += "_MCLB"; TString tit = hB->GetTitle(); tit += "_MCLB"; THnSparseF* hBMCLB = (THnSparseF*) FindObject(HName("Comb",xxZvEta),listMC)->Clone(nm.Data()); hBMCLB->SetTitle(tit.Data()); hB = hBMCLB; } THnSparseF* hBmc = 0; if (isMC) hBmc = (THnSparseF*) FindObject(HName("Comb",xxZvEta),lst); // ProjDZE(hD,hB,hBmc, outArray, useBinGrEta, useBinGrZv); // if (!isMC) return; // // prepare MC primary signal histo TH2F* mcPrim = (TH2F*)FindObject( "zvEtaPrimMC", lst ); mcPrim = (TH2F*) mcPrim->Clone("mcTrueSignal"); CropHisto(mcPrim,-useEtaCut,useEtaCut,useZvMin,useZvMax); outArray->AddAtAndExpand(mcPrim, kMCPrim + shift); // } //_____________________________________________________________________ void ProcessHistos(TObjArray* outArray) { // // build alpha matrix TH2* halp = (TH2*)outArray->At(kMCShift + kMCPrim); halp = (TH2*) halp->Clone("Alpha"); halp->SetTitle("#alpha"); halp->Divide( (TH2*)outArray->At(kMCShift + k1MBeta) ); halp->Divide( (TH2*)outArray->At(kMCShift + kRawDtCut) ); halp->SetMinimum(1.5); halp->SetMaximum(4.); outArray->AddAtAndExpand(halp, kAlpha); // // build alpha matrix with MC labels bg TH2* halpMC = (TH2*)outArray->At(kMCShift + kMCPrim); halpMC = (TH2*) halpMC->Clone("AlphaMC"); halpMC->SetTitle("#alpha MC labels"); halpMC->Divide( (TH2*)outArray->At(kMCShift + k1MBetaMC) ); halpMC->Divide( (TH2*)outArray->At(kMCShift + kRawDtCut) ); halpMC->SetMinimum(1.5); halpMC->SetMaximum(4.); outArray->AddAtAndExpand(halpMC, kAlphaMC); // // build corrected signal TH2* hsigCorr = (TH2*)outArray->At(kSignalEst); hsigCorr = (TH2*) hsigCorr->Clone("SignalEstCorr"); hsigCorr->SetTitle("Corrected Signal"); hsigCorr->Multiply( halp ); outArray->AddAtAndExpand(hsigCorr, kSigCorr); // } TCanvas *canvFin=0; TCanvas *canvAlp=0; TCanvas *canvBet=0; TCanvas* canvDst=0; TCanvas* canvSpec=0; void PlotResults() { PlotDist(); PlotAlphaBeta(); PlotSpecies(); PlotDNDEta(); } void PlotDNDEta() { // TObjArray *res = &resArr; char buff[1000]; // eta range double plEta = useEtaCut*1.1; gStyle->SetOptFit(0); gStyle->SetOptStat(0); gStyle->SetOptTitle(0); double mn = 1e6,mx = -1e6; canvFin = new TCanvas("canvFin", "canvFin",0,50,700,550); canvFin->SetLeftMargin(0.15); // canvFin->ToggleEventStatus(); // // corrected data TH1* hsigCorr = ((TH2F*)res->At(kSigCorr))->ProjectionX("DataCorrSignal"); SetHStyle(hsigCorr,kRed,20,1.0); hsigCorr->GetXaxis()->SetRangeUser(-plEta+kEps,plEta-kEps); hsigCorr->Scale(1./hsigCorr->GetBinWidth(1)); hsigCorr->Draw(); mx = TMath::Max(mx, hsigCorr->GetMaximum()); mn = TMath::Min(mn, hsigCorr->GetMinimum()); TF1* pl0 = new TF1("pl0","pol0"); pl0->SetParameter(0,hsigCorr->GetMinimum()); hsigCorr->Fit(pl0,"q0","",-0.5,.5); double fval = pl0->GetParameter(0); double ferr = pl0->GetParError(0); char ftres[1000]; sprintf(ftres,"dN/d#eta_{|#eta|<0.5} = %d #pm %d",int(fval),int(ferr)); printf("dN/d#eta_{|#eta|<0.5} = %.2f %.2f\n",fval,ferr); TLatex *txfit = new TLatex(-0.2,hsigCorr->GetMinimum()*0.9, ftres); txfit->SetTextSize(0.04); txfit->Draw(); // // raw data TH1* hraw = ((TH2F*)res->At(kRawDtCut))->ProjectionX("DataRaw"); SetHStyle(hraw,kRed,21,1.0); hraw->GetXaxis()->SetRangeUser(-plEta,plEta); hraw->Scale(1./hraw->GetBinWidth(1)); hraw->Draw("same"); mn = TMath::Min(mn, hraw->GetMinimum()); mx = TMath::Max(mx, hraw->GetMaximum()); // // raw data bg sub TH1* hraws = ((TH2F*)res->At(kSignalEst))->ProjectionX("DataRawSub"); SetHStyle(hraws,kRed,23,1.0); hraws->GetXaxis()->SetRangeUser(-plEta,plEta); hraws->Scale(1./hraw->GetBinWidth(1)); hraws->Draw("same"); mn = TMath::Min(mn, hraw->GetMinimum()); mx = TMath::Max(mx, hraw->GetMaximum()); // // bg TH1* hbg = ((TH2F*)res->At(kBgEst))->ProjectionX("BgEst"); SetHStyle(hbg,kMagenta,22,1.0); hbg->GetXaxis()->SetRangeUser(-plEta,plEta); hbg->Scale(1./hbg->GetBinWidth(1)); hbg->Draw("same"); mn = TMath::Min(mn, hbg->GetMinimum()); mx = TMath::Max(mx, hbg->GetMaximum()); // // mc part ---------------------------- // raw data TH1* hrawMC = ((TH2F*)res->At(kRawDtCut+kMCShift))->ProjectionX("DataRawMC"); SetHStyle(hrawMC,kBlue,24,1.0); hrawMC->GetXaxis()->SetRangeUser(-plEta,plEta); hrawMC->Scale(1./hrawMC->GetBinWidth(1)); hrawMC->Draw("same"); mn = TMath::Min(mn, hrawMC->GetMinimum()); mx = TMath::Max(mx, hrawMC->GetMaximum()); // // raw data bg sub TH1* hrawsMC = ((TH2F*)res->At(kSignalEst+kMCShift))->ProjectionX("DataRawSubMC"); SetHStyle(hrawsMC,kBlue,26,1.0); hrawsMC->GetXaxis()->SetRangeUser(-plEta,plEta); hrawsMC->Scale(1./hrawMC->GetBinWidth(1)); hrawsMC->Draw("same"); mn = TMath::Min(mn, hrawMC->GetMinimum()); mx = TMath::Max(mx, hrawMC->GetMaximum()); // // raw data bgMClabels sub TH1* hrawsMCLB = ((TH2F*)res->At(kSignalEstMC+kMCShift))->ProjectionX("DataRawSubMCLB"); SetHStyle(hrawsMCLB,kGreen+2,30,1.0); hrawsMCLB->GetXaxis()->SetRangeUser(-plEta,plEta); hrawsMCLB->Scale(1./hrawsMCLB->GetBinWidth(1)); hrawsMCLB->Draw("same"); mn = TMath::Min(mn, hrawsMCLB->GetMinimum()); mx = TMath::Max(mx, hrawsMCLB->GetMaximum()); // // bg est TH1* hbgMCEst = ((TH2F*)res->At(kBgEst+kMCShift))->ProjectionX("BgEstMC"); SetHStyle(hbgMCEst,kBlue,26,1.0); hbgMCEst->GetXaxis()->SetRangeUser(-plEta,plEta); hbgMCEst->Scale(1./hbgMCEst->GetBinWidth(1)); hbgMCEst->Draw("same"); mn = TMath::Min(mn, hbgMCEst->GetMinimum()); mx = TMath::Max(mx, hbgMCEst->GetMaximum()); // // bg MC TH1* hbgMC = ((TH2F*)res->At(kBgMC+kMCShift))->ProjectionX("BgMC"); SetHStyle(hbgMC,kGreen+2,25,1.0); hbgMC->GetXaxis()->SetRangeUser(-plEta,plEta); hbgMC->Scale(1./hbgMC->GetBinWidth(1)); hbgMC->Draw("same"); mn = TMath::Min(mn, hbgMC->GetMinimum()); mx = TMath::Max(mx, hbgMC->GetMaximum()); // mn = 0; hsigCorr->SetMinimum(mn); hsigCorr->SetMaximum(mx + 0.4*(mx-mn)); gPad->Modified(); // TLegend *legDnDeta = new TLegend(0.15,0.75, 0.45,0.95); legDnDeta->SetFillColor(kWhite); legDnDeta->SetHeader("Data"); legDnDeta->AddEntry(hsigCorr,"Corrected","pl"); legDnDeta->AddEntry(hraw, "Reconstructed","pl"); sprintf(buff,"Reconstructed - Bckg.%s.",useBgType.Data()); legDnDeta->AddEntry(hraws, buff,"pl"); sprintf(buff,"Background %s.",useBgType.Data()); legDnDeta->AddEntry(hbg, buff,"pl"); legDnDeta->Draw(); // TLegend *legDnDetaMC = new TLegend(0.60,0.72, 0.95,0.95); legDnDetaMC->SetFillColor(kWhite); legDnDetaMC->SetHeader("MC"); legDnDetaMC->AddEntry(hrawMC, "Reconstructed","pl"); sprintf(buff,"Reconstructed - Bckg.%s.",useBgType.Data()); legDnDetaMC->AddEntry(hrawsMC, buff,"pl"); sprintf(buff,"Reconstructed - Bckg.%s.","MC.Labels"); legDnDetaMC->AddEntry(hrawsMCLB, buff,"pl"); sprintf(buff,"Background %s.",useBgType.Data()); legDnDetaMC->AddEntry(hbgMCEst, buff,"pl"); sprintf(buff,"Background %s.","MC Labels"); legDnDetaMC->AddEntry(hbgMC, buff,"pl"); legDnDetaMC->Draw(); // gPad->SetGrid(1.1); gPad->Modified(); canvFin->cd(); AddLabel(outTitle,0.1,0.97, kBlack,0.025); sprintf(buff,"%s/%sdNdEta_%s",figDir,uniqueName.Data(),outStr); SaveCanvas(canvFin,buff,"cg"); // } // void PlotAlphaBeta() { char buff[1000]; gStyle->SetOptFit(0); gStyle->SetOptStat(0); gStyle->SetOptTitle(0); TObjArray* res = &resArr; //------------------------------------------------------ canvBet = new TCanvas("canvBet","canvBet",10,10,1000,400); canvBet->Divide(3,1,0.01,0.06); canvBet->cd(1); TH1* dtBet = (TH1*)res->At(k1MBeta); TH1* mcBet = (TH1*)res->At(k1MBeta+kMCShift); TH1* mcBetLB = (TH1*)res->At(k1MBetaMC+kMCShift); double mn,mx,mnt,mxt; GetRealMinMax(dtBet,mn,mx); GetRealMinMax(mcBet,mnt,mxt); if (mntmx) mx = mxt; GetRealMinMax(mcBetLB,mnt,mxt); if (mntmx) mx = mxt; // dtBet->SetMinimum(mn - 0.05*(mx-mn)); dtBet->SetMaximum(mx + 0.05*(mx-mn)); mcBet->SetMinimum(mn - 0.05*(mx-mn)); mcBet->SetMaximum(mx + 0.05*(mx-mn)); mcBetLB->SetMinimum(mn - 0.05*(mx-mn)); mcBetLB->SetMaximum(mx + 0.05*(mx-mn)); // canvBet->cd(1); gPad->SetRightMargin(0.15); dtBet->Draw("colz"); AddLabel("#beta Data",0.2,0.95,kBlack,0.05); gPad->Modified(); dtBet->GetYaxis()->SetTitleOffset(1.4); TPaletteAxis *p = (TPaletteAxis*)dtBet->FindObject("palette"); if (p) p->SetX1NDC(0.85); canvBet->cd(2); gPad->SetRightMargin(0.15); mcBet->Draw("colz"); AddLabel("#beta MC (bckg.estimated)",0.2,0.95,kBlack,0.05); gPad->Modified(); mcBet->GetYaxis()->SetTitleOffset(1.4); p = (TPaletteAxis*)mcBet->FindObject("palette"); if (p) p->SetX1NDC(0.85); canvBet->cd(3); gPad->SetRightMargin(0.15); mcBetLB->Draw("colz"); AddLabel("#beta MC (bckg.from MC labels)",0.2,0.95,kBlack,0.05); gPad->Modified(); mcBetLB->GetYaxis()->SetTitleOffset(1.4); p = (TPaletteAxis*)mcBetLB->FindObject("palette"); if (p) p->SetX1NDC(0.85); // canvBet->cd(); AddLabel(outTitle,0.1,0.97, kBlack, 0.025); // sprintf(buff,"%s/%sBeta_%s",figDir,uniqueName.Data(),outStr); SaveCanvas(canvBet,buff,"cg"); //------------------------------------------------------ canvAlp = new TCanvas("canvAlp","canvAlp",10,10,900,400); canvAlp->Divide(2,1,0.01,0.06); canvAlp->cd(1); TH1* dtAlp = (TH1*)res->At(kAlpha); TH1* mcAlp = (TH1*)res->At(kAlphaMC); GetRealMinMax(dtAlp,mn,mx); GetRealMinMax(mcAlp,mnt,mxt); if (mntmx) mx = mxt; dtAlp->SetMinimum(mn - 0.05*(mx-mn)); dtAlp->SetMaximum(mx + 0.05*(mx-mn)); mcAlp->SetMinimum(mn - 0.05*(mx-mn)); mcAlp->SetMaximum(mx + 0.05*(mx-mn)); // canvAlp->cd(1); gPad->SetRightMargin(0.15); dtAlp->Draw("colz"); AddLabel("#alpha (bckg.estimated)",0.2,0.95,kBlack,0.05); gPad->Modified(); dtAlp->GetYaxis()->SetTitleOffset(1.4); TPaletteAxis *pa = (TPaletteAxis*)dtBet->FindObject("palette"); if (pa) pa->SetX1NDC(0.85); canvAlp->cd(2); gPad->SetRightMargin(0.15); mcAlp->Draw("colz"); AddLabel("#alpha (bckg.from MC labels)",0.2,0.95,kBlack,0.05); gPad->Modified(); mcAlp->GetYaxis()->SetTitleOffset(1.4); pa = (TPaletteAxis*)mcBet->FindObject("palette"); if (pa) pa->SetX1NDC(0.85); gPad->Modified(); canvAlp->cd(); AddLabel(outTitle,0.1,0.97, kBlack, 0.025); // sprintf(buff,"%s/%sAlpha_%s",figDir,uniqueName.Data(),outStr); SaveCanvas(canvAlp,buff,"cg"); } void PlotSpecies() { char buff[1000]; gStyle->SetOptFit(0); gStyle->SetOptStat(0); gStyle->SetOptTitle(0); //------------------------------------------------------ TH2F* hSpecPrim = (TH2F*)FindObject( "pdgPrim", listMC,kFALSE); TH2F* hSpecSec = (TH2F*)FindObject( "pdgSec", listMC,kFALSE); TH2F* hSpecPrimP = (TH2F*)FindObject( "pdgPrimPar", listMC,kFALSE); TH2F* hSpecSecP = (TH2F*)FindObject( "pdgSecPar", listMC,kFALSE); int nbd = hSpecPrim->GetXaxis()->GetNbins(); // TH1* hSpecPrimAll = hSpecPrim->ProjectionX("specPrimAll",1,nbd+1,"e"); hSpecPrimAll->Scale(100./hSpecPrimAll->Integral()); hSpecPrimAll->GetYaxis()->SetTitle("Fraction,%"); hSpecPrimAll->GetXaxis()->SetLabelSize(0.06); hSpecPrimAll->GetXaxis()->LabelsOption("v"); // TH1* hSpecSecAll = hSpecSec->ProjectionX("specSecAll",1,nbd+1,"e"); hSpecSecAll->Scale(100./hSpecSecAll->Integral()); hSpecSecAll->GetYaxis()->SetTitle("Fraction,%"); hSpecSecAll->GetXaxis()->SetLabelSize(0.05); // TH1* hSpecPrimPAll = hSpecPrimP->ProjectionX("specPrimPAll",1,nbd+1,"e"); hSpecPrimPAll->Scale(100./hSpecPrimPAll->Integral()); hSpecPrimPAll->GetYaxis()->SetTitle("Fraction,%"); hSpecPrimPAll->GetXaxis()->SetLabelSize(0.06); hSpecPrimPAll->GetXaxis()->LabelsOption("v"); // TH1* hSpecSecPAll = hSpecSecP->ProjectionX("specSecPAll",1,nbd+1,"e"); hSpecSecPAll->Scale(100./hSpecSecPAll->Integral()); hSpecSecPAll->GetYaxis()->SetTitle("Fraction,%"); hSpecSecPAll->GetXaxis()->SetLabelSize(0.05); // int binCut = hSpecPrim->GetXaxis()->FindBin(kWDistSgCut-kEps); TH1* hSpecPrimSel = hSpecPrim->ProjectionX("specPrimSel",1,binCut,"e"); hSpecPrimSel->Scale(100./hSpecPrimSel->Integral()); hSpecPrimSel->GetYaxis()->SetTitle("Fraction,%"); hSpecPrimSel->GetXaxis()->SetLabelSize(0.05); // TH1* hSpecSecSel = hSpecSec->ProjectionX("specSecSel",1,binCut,"e"); hSpecSecSel->Scale(100./hSpecSecSel->Integral()); hSpecSecSel->GetYaxis()->SetTitle("Fraction,%"); hSpecSecSel->GetXaxis()->SetLabelSize(0.05); // TH1* hSpecPrimPSel = hSpecPrimP->ProjectionX("specPrimPSel",1,binCut,"e"); hSpecPrimPSel->Scale(100./hSpecPrimPSel->Integral()); hSpecPrimPSel->GetYaxis()->SetTitle("Fraction,%"); hSpecPrimPSel->GetXaxis()->SetLabelSize(0.05); // TH1* hSpecSecPSel = hSpecSecP->ProjectionX("specSecPSel",1,binCut,"e"); hSpecSecPSel->Scale(100./hSpecSecPSel->Integral()); hSpecSecPSel->GetYaxis()->SetTitle("Fraction,%"); hSpecSecPSel->GetXaxis()->SetLabelSize(0.05); canvSpec = new TCanvas("canvSpec","canvSpec",10,10,1100,800); canvSpec->Divide(1,2,0.01,0.01); canvSpec->cd(1); hSpecPrimAll->Draw(); SetHStyle(hSpecPrimAll,kBlue,25,1.1); hSpecPrimSel->Draw("same"); SetHStyle(hSpecPrimSel,kRed,20,1); // hSpecSecAll->Draw("same"); SetHStyle(hSpecSecAll,kGreen,32,1.1); hSpecSecSel->Draw("same"); SetHStyle(hSpecSecSel,kBlack,22,1); // TLegend *legPart = new TLegend(0.8,0.72, 0.999,0.999); legPart->SetFillColor(kWhite); legPart->SetHeader("Tracklet PDG"); // legPart->AddEntry(hSpecPrimAll, "Prim., before #Delta cut","pl"); legPart->AddEntry(hSpecPrimSel, "Prim., after #Delta cut","pl"); legPart->AddEntry(hSpecSecAll, "Sec., before #Delta cut","pl"); legPart->AddEntry(hSpecSecSel, "Sec., after #Delta cut","pl"); // legPart->Draw(); gPad->SetLogy(); gPad->SetGrid(1,1); gPad->Modified(); // canvSpec->cd(2); hSpecPrimPAll->Draw(); SetHStyle(hSpecPrimPAll,kBlue,25,1.1); hSpecPrimPSel->Draw("same"); SetHStyle(hSpecPrimPSel,kRed,20,1); // hSpecSecPAll->Draw("same"); SetHStyle(hSpecSecPAll,kGreen,32,1.1); hSpecSecPSel->Draw("same"); SetHStyle(hSpecSecPSel,kBlack,22,1); // TLegend *legPartP = new TLegend(0.8,0.72, 0.999,0.999); legPartP->SetFillColor(kWhite); legPartP->SetHeader("Tracklet Parents PDG"); // legPartP->AddEntry(hSpecPrimPAll, "Prim., before #Delta cut","pl"); legPartP->AddEntry(hSpecPrimPSel, "Prim., after #Delta cut","pl"); legPartP->AddEntry(hSpecSecPAll, "Sec., before #Delta cut","pl"); legPartP->AddEntry(hSpecSecPSel, "Sec., after #Delta cut","pl"); // legPartP->Draw(); gPad->SetLogy(); gPad->SetGrid(1,1); gPad->Modified(); // canvSpec->cd(1); AddLabel(outTitle,0.1,0.97, kBlack, 0.025); canvSpec->cd(); // sprintf(buff,"%s/%sSpecies_%s",figDir,uniqueName.Data(),outStr); SaveCanvas(canvSpec,buff,"cg"); } void PlotDist() { TObjArray* res = &resArr; char buff[1000]; canvDst = new TCanvas("canvDst","canvDst",10,10,700,500); TH1* mcdst = (TH1*)res->At(kDataDist+kMCShift); TH1* mcdstbg = (TH1*)res->At(kBgDist+kMCShift); TH1* mcdstbgLB = (TH1*)res->At(kBgMCDist+kMCShift); TH1* dtdst = (TH1*)res->At(kDataDist); TH1* dtdstbg = (TH1*)res->At(kBgDist); // TH2* mcDstZN = (TH2*)FindObject(useShapeType==kNormShapeDist ? "TrData_ZvDist":"TrData_ZvDPhiS", listMC ); TH2* mcDstZSec = (TH2*)FindObject(useShapeType==kNormShapeDist ? "TrSec_ZvDist":"TrSec_ZvDPhiS", listMC ); TH2* mcDstZCombU = (TH2*)FindObject(useShapeType==kNormShapeDist ? "TrCombU_ZvDist":"TrCombU_ZvDPhiS", listMC ); TH2* mcDstZCombC = (TH2*)FindObject(useShapeType==kNormShapeDist ? "TrComb_ZvDist":"TrComb_ZvDPhiS", listMC ); // mcDstZN->GetXaxis()->SetRangeUser(useZvMin+kEps,useZvMax-kEps); mcDstZSec->GetXaxis()->SetRangeUser(useZvMin+kEps,useZvMax-kEps); mcDstZCombU->GetXaxis()->SetRangeUser(useZvMin+kEps,useZvMax-kEps); mcDstZCombC->GetXaxis()->SetRangeUser(useZvMin+kEps,useZvMax-kEps); // TH1* mcDstN = mcDstZN->ProjectionY("mcDstN"); TH1* mcDstSec = mcDstZSec->ProjectionY("mcDstSec"); TH1* mcDstCombU = mcDstZCombU->ProjectionY("mcDstCombU"); TH1* mcDstCombC = mcDstZCombC->ProjectionY("mcDstCombC"); // double scl,sclE; mcDstN = NormalizeBg(mcdst,mcDstN,scl,sclE); mcDstSec->Scale(scl); mcDstCombU->Scale(scl); mcDstCombC->Scale(scl); mcDstCombC->Add(mcDstCombU,-1); // dtdst->Draw(""); gPad->Modified(); dtdst->GetXaxis()->SetLabelSize(0.03); dtdst->GetXaxis()->SetTitleSize(0.03); dtdst->GetXaxis()->SetTitleOffset(2); dtdstbg->Draw("same"); mcdst->Draw("same"); mcDstSec->Draw("same"); mcdstbgLB->Draw("same"); mcdstbg->Draw("same"); mcDstCombC->Draw("same"); // SetHStyle(mcdst,kBlue, 25,0.7); SetHStyle(mcdstbgLB,kGreen, 7/*32*/,0.5); SetHStyle(mcdstbg,kCyan, 7/*26*/,0.5); SetHStyle(mcDstCombC,kGreen+2, 21,0.7); SetHStyle(mcDstSec,kBlue+2, 22,0.7); // SetHStyle(dtdst,kRed, 20,0.7); SetHStyle(dtdstbg,kBlue, 34,0.7); // double vmcTot,vmcTotE; double vmcSec,vmcSecE, ratSec,ratSecE; double vmcCmbEst,vmcCmbEstE, ratCmbEst,ratCmbEstE; double vmcCmb,vmcCmbE, ratCmb,ratCmbE; double vmcCmbC,vmcCmbCE, ratCmbC,ratCmbCE; double cutSgMin,cutSgMax; double cutBgMin,cutBgMax; if (useShapeType==kNormShapeDist) { cutSgMin = 0; cutSgMax = kWDistSgCut; cutBgMin = kWDistBgTailMin; cutBgMax = kWDistBgTailMax; } else { cutSgMin = 0; cutSgMax = kdPhiSgCut; cutBgMin = kdPhiBgTailMin; cutBgMax = kdPhiBgTailMax; } Integrate(mcdst, cutSgMin,cutSgMax, vmcTot,vmcTotE); Integrate(mcDstSec, cutSgMin,cutSgMax, vmcSec,vmcSecE); GetRatE(vmcSec,vmcSecE, vmcTot,vmcTotE, ratSec,ratSecE); // Integrate(mcdstbgLB, cutSgMin,cutSgMax, vmcCmb,vmcCmbE); GetRatE(vmcCmb,vmcCmbE, vmcTot,vmcTotE, ratCmb,ratCmbE); // Integrate(mcdstbg, cutSgMin,cutSgMax, vmcCmbEst,vmcCmbEstE); GetRatE(vmcCmbEst,vmcCmbEstE, vmcTot,vmcTotE, ratCmbEst,ratCmbEstE); // Integrate(mcDstCombC, cutSgMin,cutSgMax, vmcCmbC,vmcCmbCE); GetRatE(vmcCmbC,vmcCmbCE, vmcTot,vmcTotE, ratCmbC,ratCmbCE); // double vdtTot,vdtTotE; double vdtBg,vdtBgE, ratdtBg,ratdtBgE; // Integrate(dtdst, cutSgMin,cutSgMax, vdtTot,vdtTotE); Integrate(dtdstbg, cutSgMin,cutSgMax, vdtBg,vdtBgE); GetRatE( vdtBg,vdtBgE, vdtTot,vdtTotE, ratdtBg,ratdtBgE); // // double dmn = mcdst->GetMinimum(); double dmx = mcdst->GetMaximum(); TLine *ln = new TLine(cutSgMax, dmn, cutSgMax, dmx); ln->SetLineColor(kBlack); ln->Draw(); TLine *lnc = new TLine(cutBgMin, dmn, cutBgMin, dmx); lnc->SetLineColor(kRed); lnc->Draw(); if (useShapeType==kNormShapeDPhi) { ln = new TLine(-cutSgMax, dmn, -cutSgMax, dmx); ln->SetLineColor(kBlack); ln->Draw(); // lnc = new TLine(-cutBgMin, dmn, -cutBgMin, dmx); lnc->SetLineColor(kRed); lnc->Draw(); } // TLegend *legDstMC1 = new TLegend(0.60,0.72, 0.95,0.95); legDstMC1->SetFillColor(kWhite); // legDstMC1->AddEntry(dtdst, "Data raw","pl"); sprintf(buff,"Data Comb. %s. : %.1f%%",useBgType.Data(),ratdtBg*100); legDstMC1->AddEntry(dtdstbg, buff,"pl"); // legDstMC1->AddEntry(mcdst, "MC raw","pl"); sprintf(buff,"MC Comb. %s. : %.1f%%",useBgType.Data(),ratCmbEst*100); legDstMC1->AddEntry(mcdstbg, buff,"pl"); // sprintf(buff,"MC Comb. %s. : %.1f%%","MC Lbl.",ratCmb*100); legDstMC1->AddEntry(mcdstbgLB, buff,"pl"); sprintf(buff,"MC Comb.Uncorr %s. : %.1f%%","MC Lbl.",ratCmbC*100); legDstMC1->AddEntry(mcDstCombC, buff,"pl"); sprintf(buff,"MC Sec. : %.1f%%",ratSec*100); legDstMC1->AddEntry(mcDstSec, buff,"pl"); legDstMC1->Draw(); if (useShapeType==kNormShapeDist) gPad->SetLogx(); gPad->SetLogy(); gPad->SetGrid(1,1); gPad->Modified(); // canvDst->cd(); AddLabel(outTitle,0.1,0.97, kBlack, 0.025); // sprintf(buff,"%s/%sDst_%s",figDir,uniqueName.Data(),outStr); SaveCanvas(canvDst,buff,"cg"); // } //______________________________________________________________________ void CropHisto(TH1* histo, int bx0, int bx1, int by0, int by1) { // fill 0 to all bins outside defined range TAxis *xax = histo->GetXaxis(); int nbx = xax->GetNbins(); double vmn=1e16,vmx=-1e16; if (histo->InheritsFrom(TH2::Class())) { TAxis *yax = histo->GetYaxis(); int nby = yax->GetNbins(); for (int ix=nbx+2;ix--;) { for (int iy=nby+2;iy--;) { if ((ixbx1)||(iyby1)) { histo->SetBinContent(ix,iy,0); histo->SetBinError(ix,iy,0); } else { double vl = histo->GetBinContent(ix,iy); if (vlvmx) vmx = vl; } } } } else { for (int ix=nbx+2;ix--;) { if ((ixbx1)) { histo->SetBinContent(ix,0); histo->SetBinError(ix,0); } else { double vl = histo->GetBinContent(ix); if (vlvmx) vmx = vl; } } } // if (vmn==vmx) { vmn = 0.95*vmn; vmx = 1.05*vmx; } histo->SetMaximum(vmx); histo->SetMinimum(vmn); } //______________________________________________________________________ void CropHisto(TH1* histo, double vx0, double vx1, double vy0, double vy1) { // fill 0 to all bins outside defined range TAxis *xax = histo->GetXaxis(); int bx0,bx1,by0=-1,by1=-1; bx0 = xax->FindBin(vx0+kEps); bx1 = xax->FindBin(vx1-kEps); if (histo->InheritsFrom(TH2::Class())) { TAxis *yax = histo->GetYaxis(); by0 = yax->FindBin(vy0+kEps); by1 = yax->FindBin(vy1-kEps); } CropHisto(histo,bx0,bx1,by0,by1); } //______________________________________________________________________ TH1* NormalizeBg(TH1* dataH, TH1* bgH, double &scl, double &sclE) { // match generated bg and data tails, calculate normalization, return normalized bg copy // TAxis* xax = dataH->GetXaxis(); int nbtot = xax->GetNbins(); int bgBins[2][2] = {{0}}; // limiting bins for tails integration Int_t ntails; // 0 for dphi plot, 1 for weighted dist plot if (useShapeType == kNormShapeDist) { // only positive tail bgBins[0][0] = xax->FindBin(kWDistBgTailMin+kEps); // positive tail min bin bgBins[0][1] = xax->FindBin(kWDistBgTailMax-kEps); // positive tail max bin ntails = 1; } else if (useShapeType == kNormShapeDPhi) { // both tails bgBins[0][0] = xax->FindBin( kdPhiBgTailMin+kEps); // positive tail min bin bgBins[0][1] = xax->FindBin( kdPhiBgTailMax-kEps); // positive tail max bin bgBins[1][0] = xax->FindBin(-kdPhiBgTailMax+kEps); // negative tail min bin bgBins[1][1] = xax->FindBin(-kdPhiBgTailMin-kEps); // positive tail max bin ntails = 2; } else {printf("NormalizeBg: unknown shape type %d\n",useShapeType);exit(1);} printf("NormalizeBg: bins for tails: right: %d:%d / left: %d:%d\n",bgBins[0][0],bgBins[0][1],bgBins[1][0],bgBins[1][1]); // double meanR=0,meanRE=0,meanRE2=0; double meanD=0,meanDE2=0; double meanB=0,meanBE2=0; double meanRI=0,meanRIE=0; for (int itp=0;itp<=ntails;itp++) { for (int ib=bgBins[itp][0];ib<=bgBins[itp][1];ib++) { if (ib<1||ib>nbtot) continue; double vD = dataH->GetBinContent(ib); double vB = bgH->GetBinContent(ib); double eD = dataH->GetBinError(ib); double eB = bgH->GetBinError(ib); meanD += vD; meanDE2 += eD*eD; meanB += vB; meanBE2 += eB*eB; if (vD<=0 || vB<=0 || eD<=0 || eB<=0) continue; double rat = vD/vB; double ratE2 = rat*rat*(eD*eD/vD/vD + eB*eB/vB/vB); meanR += rat/ratE2; meanRE2 += 1.0/ratE2; } } // if (meanRE2>0) { meanR /= meanRE2; meanRE2 = 1./meanRE2; meanRE = TMath::Sqrt(meanRE2); } if (meanDE2>0 && meanBE2>0) { meanRI = meanD/meanB; meanRIE = meanRI*TMath::Sqrt(meanDE2/meanD/meanD + meanBE2/meanB/meanB); } printf("NormalizeBg: Tails scaling %s wrt %s: Wgh.Mean:%.4f(%.4f) / Integral:%.4f(%.4f)\n", bgH->GetName(),dataH->GetName(), meanR,meanRE, meanRI,meanRIE); printf("NormalizeBg: Select scaling type %s\n",useScaleType==kSclWghMean ? "Wgh.Mean":"Integral"); // scl = useScaleType==kSclWghMean ? meanR : meanRI; sclE = useScaleType==kSclWghMean ? meanRE : meanRIE; // // rescaled bg char buff[1000]; sprintf(buff,"%s_bgNorm",bgH->GetName()); bgH = (TH1*)bgH->Clone(buff); sprintf(buff,"%s bgNorm%d %.4f+-%.4f",bgH->GetName(),useScaleType,scl,sclE); TH1* dumH = (TH1*)bgH->Clone("dummySCL$"); dumH->Reset(); for (int i=1;i<=nbtot;i++) { dumH->SetBinContent(i,scl); dumH->SetBinError(i,sclE); } bgH->Multiply(dumH); delete dumH; return bgH; } //______________________________________________________________________ TObject* FindObject(const char* nameH, const TList* lst, Bool_t normPerEvent) { // get histo, optionally normalizing it per processed event if (!lst) {printf("FindObject %s: No list provided\n",nameH); exit(1);} int nent = lst->GetEntries(); TString nm; TObject *hst = 0; for (int i=nent;i--;) { nm = lst->At(i)->GetName(); if (nm.EndsWith(nameH)) {hst = lst->At(i); break;} } if (!hst) {printf("FindObject: No %s histo in list %s\n",nameH,lst->GetName()); exit(1);} if (!normPerEvent || hst->TestBit(kBitNormPerEvent)) return hst; // already normalized TString nameHS = nameH; if (nameHS==kHStatName) return hst; // never scale stat. histo // TH1* hstat = (TH1*)FindObject(kHStatName,lst,kFALSE); double nrm = hstat->GetBinContent(kEvProcData); if (nrm<1) {printf("FindObject: Anomaluous %d number of events processed in list %p\n",int(nrm),lst); exit(1);} // // account for eventual cut in Z TH1* hzv = (TH1*)FindObject("zv",lst,kFALSE); int nbz = hzv->GetNbinsX(); double zvTot = hzv->Integral(1,nbz); int zb0 = hzv->FindBin( useZvMin+kEps); if (zb0<1) zb0 = 1; int zb1 = hzv->FindBin( useZvMax-kEps); if (zb1>nbz) zb1 = nbz; double zvSel = hzv->Integral(zb0,zb1); if (zvTot<1 || zvSel<1) {printf("No statistics: NzvTot: %.1f NzvSel:%.f\n",zvTot,zvSel); return 0;} else printf("%f fraction of selected events is used with current Zv cut %.1f:%.1f\n",zvSel/zvTot,useZvMin,useZvMax); nrm *= zvSel/zvTot; // if (hst->InheritsFrom(TH1::Class())) ((TH1*)hst)->Scale(1./nrm); else if (hst->InheritsFrom(THnSparse::Class())) { THnSparse* spr = (THnSparse*) hst; spr->Sumw2(); int coord[3] = {0,0,0}; for (Long64_t i = 0; i < spr->GetNbins(); ++i) { // Get the content of the bin from the current histogram Double_t v = spr->GetBinContent(i, coord); spr->SetBinContent(coord, v/nrm); spr->SetBinError(coord,TMath::Sqrt(v)/nrm); } } // hst->SetBit(kBitNormPerEvent); return hst; } //______________________________________________________________________ TList* LoadList(const char* flName, const char* addPref, const char* nameL) { // load list with histos TString nms = flName; gSystem->ExpandPathName(nms); TFile* fl = TFile::Open(nms.Data()); if (!fl) {printf("LoadList: No file %s\n",nms.Data()); exit(1);} TList* lst = (TList*)fl->Get(nameL); if (!lst) {printf("LoadList: No list %s in file %s\n",nameL,nms.Data()); exit(1);} lst->SetName(flName); // int nEnt = lst->GetSize(); TString nm; for (int i=0;iAt(i); nm = addPref; nm += ob->GetName(); ob->SetName(nm.Data()); } // return lst; } //____________________________________________________________________________ void GetRatE(double x,double xe, double y,double ye, double &rat, double &rate) { rat = 0; rate = 0; if (TMath::Abs(y)<1e-16 || TMath::Abs(x)<1e-16) return; rat = x/y; rate = rat*TMath::Sqrt( xe*xe/(x*x) + ye*ye/(y*y)); } //____________________________________________________________________________ void Integrate(TH1* hist, double xmn,double xmx, double &val, double& err) { // integrate 1d histo within given limits TAxis* xax = hist->GetXaxis(); int bmn = xax->FindBin(xmn+kEps); if (bmn<1) bmn = 0; // include int bmx = xax->FindBin(xmx-kEps); val = hist->IntegralAndError(bmn,bmx,err); // is this histo with symmetric axis ? then integrate also negative half axis if (TMath::Abs( xax->GetXmin() + xax->GetXmax() )<1e-6) { bmn = xax->FindBin(-xmx+kEps); bmx = xax->FindBin(-xmn-kEps); double errn; val += hist->IntegralAndError(bmn,bmx,errn); err = TMath::Sqrt(err*err + errn*errn); } } //____________________________________________________________________________ const char* HName(const char* prefix,const char* htype) { // compose the name of histo in the clist static TString strh; strh = "Tr"; strh += prefix; strh += "_"; strh += htype; return strh.Data(); } //____________________________________________________________________________ Int_t CheckStat(const TList* lst, const char* dtType) { // check if needed bg was generated TH1* hstat = (TH1*)FindObject(kHStatName,lst); TString dts = dtType; if (dts=="Data") return int( hstat->GetBinContent(kEvProcData) ); if (dts=="Mix") return int( hstat->GetBinContent(kEvProcMix) ); if (dts=="Inj") return int( hstat->GetBinContent(kEvProcInj) ); if (dts=="Rot") return int( hstat->GetBinContent(kEvProcRot) ); printf("Unknown process %s statistics is checked. Alowed: Data,Mix,Inj,Rot",dtType); return 0; } //____________________________________________________________________________ void ProjDZE(THnSparseF* hrawd, THnSparseF* hgenb, THnSparseF* hmcComb, TObjArray* res, int bStepEta,int bStepZv) { // project 3d histo of Dist vs Zv vs Eta to Zv vs Eta with all cuts int shift = hmcComb ? kMCShift : 0; // is this mc? // // determine boundaries for zv and eta cuts // double cutDstMin,cutDstMax; double cutBgMin,cutBgMax; double cutSgMin,cutSgMax; if (useShapeType==kNormShapeDist) { cutDstMin = 0; cutDstMax = kWDistBgTailMax; // cutBgMin = kWDistBgTailMin; cutBgMax = kWDistBgTailMax; // cutSgMin = 0; cutSgMax = kWDistSgCut; } else { cutDstMin = -kdPhiBgTailMax; cutDstMax = kdPhiBgTailMax; // cutBgMin = kdPhiBgTailMin; cutBgMax = kdPhiBgTailMax; // cutSgMin = 0; cutSgMax = kdPhiSgCut; } // int bn0[3] = {0}; // 1st bin to count int bn1[3] = {0}; // last bin to count // eta TAxis* axEta = hrawd->GetAxis(0); bn0[0] = axEta->FindBin(-useEtaCut+kEps); bn1[0] = axEta->FindBin( useEtaCut-kEps); // Zv TAxis* axZv = hrawd->GetAxis(1); bn0[1] = axZv->FindBin( useZvMin+kEps); bn1[1] = axZv->FindBin( useZvMax-kEps); // W.dist TAxis* axDst = hrawd->GetAxis(2); bn0[2] = axDst->FindBin( cutDstMin + kEps); bn1[2] = axDst->FindBin( cutDstMax - kEps); // // int nb[3] = { bn1[0]-bn0[0]+1, bn1[1]-bn0[1]+1, bn1[2]-bn0[2]+1}; // number of bins to count int nbTot[3] = {axEta->GetNbins(), axZv->GetNbins(), axDst->GetNbins() }; // total bins // if (bStepEta<1 || bStepEta>nb[0]) bStepEta = nb[0]; if (bStepZv<1 || bStepEta>nb[1]) bStepZv = nb[1]; // for (int i=3;i--;) { // set default axis range hrawd->GetAxis(i)->SetRange(1,nbTot[i]); hgenb->GetAxis(i)->SetRange(1,nbTot[i]); if (hmcComb) hmcComb->GetAxis(i)->SetRange(1,nbTot[i]); } // char grpTit[100]; sprintf(grpTit,"grp_eta%d_zv%d",bStepEta,bStepZv); // TString pref = hmcComb ? "mc" : "dt"; // // "Data" histo with cut on tails where we look for signal hrawd->GetAxis(2)->SetRangeUser(cutDstMin+kEps,cutDstMax-kEps); TH2* hRawDtCut = hrawd->Projection(1,0,"e"); hrawd->GetAxis(2)->SetRange(bn0[2],bn1[2]); hRawDtCut->SetName(pref+"_RawWithCut"); hRawDtCut->SetTitle(pref+" Raw with cut on tracklets"); res->AddAtAndExpand(hRawDtCut, kRawDtCut+shift); // // "Data - Est.Bg" histo with cut on tails where we look for signal hrawd->GetAxis(2)->SetRangeUser(cutDstMin+kEps,cutDstMax-kEps); TH2* hSignalEst = hrawd->Projection(1,0,"e"); hrawd->GetAxis(2)->SetRange(bn0[2],bn1[2]); hSignalEst->SetName(pref+"_SignalCut_"+grpTit); hSignalEst->SetTitle(pref+" Signal (raw-bg) with cut on tracklets "+grpTit); res->AddAtAndExpand(hSignalEst, kSignalEst+shift); // // "Data - MC.Bg" histo with cut on tails where we look for signal TH2* hSignalEstMC = 0; if (hmcComb) { hrawd->GetAxis(2)->SetRangeUser(cutDstMin+kEps,cutDstMax-kEps); hSignalEstMC = hrawd->Projection(1,0,"e"); hrawd->GetAxis(2)->SetRange(bn0[2],bn1[2]); hSignalEstMC->SetName(pref+"_SignalCut_bgMCLabels_"+grpTit); hSignalEstMC->SetTitle(pref+" Signal (raw-bg_MCLabels) with cut on tracklets "+grpTit); res->AddAtAndExpand(hSignalEstMC, kSignalEstMC+shift); } // // Estimated background in the cut range hgenb->GetAxis(2)->SetRangeUser(cutDstMin+kEps,cutDstMax-kEps); TH2* hBgEst = hgenb->Projection(1,0,"e"); hgenb->GetAxis(2)->SetRange(bn0[2],bn1[2]); hBgEst->SetName(pref+"_BgEst_"+grpTit); hBgEst->SetTitle(pref+" Estimated Bg "+grpTit); res->AddAtAndExpand(hBgEst, kBgEst+shift); // // 1-beta for "data" = (Data_cut - Bg_cut) / Data_cut TH2* h1mBeta = hrawd->Projection(1,0,"e"); h1mBeta->Reset(); h1mBeta->SetName(pref+"_h1mBeta_"+grpTit); h1mBeta->SetTitle(pref+" 1-#beta with gen.bg. "+grpTit); res->AddAtAndExpand(h1mBeta, k1MBeta+shift); // // If MC labels info is provided TH2* h1mBetaMC = 0; // 1-beta for MC with bg from labels TH2* hBgMC = 0; // bg from MC labels if (hmcComb) { hmcComb->GetAxis(2)->SetRangeUser(cutDstMin+kEps,cutDstMax-kEps); h1mBetaMC = hmcComb->Projection(1,0,"e"); h1mBetaMC->SetName(pref+"_h1mBetaMC"); h1mBetaMC->SetTitle(pref+" 1-#beta with bg. from MC labels"); h1mBetaMC->Divide(hRawDtCut); res->AddAtAndExpand(h1mBetaMC, k1MBetaMC+shift); for (int ib0=1;ib0<=nbTot[0];ib0++) for (int ib1=1;ib1<=nbTot[1];ib1++) h1mBetaMC->SetBinContent(ib0,ib1, 1.- h1mBetaMC->GetBinContent(ib0,ib1)); // hBgMC = hmcComb->Projection(1,0,"e"); hBgMC->SetName(pref+"_Bg_MClab"); hBgMC->SetTitle(pref+" Bg from MC labels"); res->AddAtAndExpand(hBgMC, kBgMC+shift); // // finalize estimated signal with bg from MC labels hSignalEstMC->Add(hBgMC,-1); // hmcComb->GetAxis(2)->SetRange(bn0[2],bn1[2]); } // // rescaling factors for generated bg TH2* hBgRescFc = hrawd->Projection(1,0,"e"); hBgRescFc->Reset(); hBgRescFc->SetName(pref+"_hBgRescFactors_"+grpTit); hBgRescFc->SetTitle(pref+" Scale.factor for gen.bg. "+grpTit); res->AddAtAndExpand(hBgRescFc, kBgRescFc+shift); // int nbint = bStepEta*bStepZv; float nbinstsq = TMath::Sqrt(nbint); hrawd->GetAxis(0)->SetRange(bn0[0],bn1[0]); hrawd->GetAxis(1)->SetRange(bn0[1],bn1[1]); TH1* hDstDt = hrawd->Projection(2,"e"); hrawd->GetAxis(0)->SetRange(1,nbTot[0]); hrawd->GetAxis(1)->SetRange(1,nbTot[1]); hDstDt->SetName(pref+"_DistRaw_"+grpTit); // "Data" projection on tracklet quality (w.distance) axis to check the tails hDstDt->SetTitle(pref+" DistRaw "+grpTit); double nrmDst,dumErr = 0; Integrate(hDstDt, cutBgMin,cutBgMax, nrmDst, dumErr); hDstDt->Scale(1./nrmDst); res->AddAtAndExpand(hDstDt, kDataDist+shift); // hgenb->GetAxis(0)->SetRange(bn0[0],bn1[0]); hgenb->GetAxis(1)->SetRange(bn0[1],bn1[1]); TH1* hDstBg = hgenb->Projection(2,"e"); hDstBg->Reset(); hgenb->GetAxis(0)->SetRange(1,nbTot[0]); hgenb->GetAxis(1)->SetRange(1,nbTot[1]); hDstBg->SetName(pref+"_DistBgNorm_"+grpTit); // Gen.Bg projection on tracklet quality (w.distance) axis to check the tails hDstBg->SetTitle(pref+" DistBgNorm "+grpTit); // Gen.Bg projection on tracklet quality (w.distance) axis to check the tails res->AddAtAndExpand(hDstBg, kBgDist+shift); TH1* hDstBgMC = 0; if (hmcComb) { hmcComb->GetAxis(0)->SetRange(bn0[0],bn1[0]); hmcComb->GetAxis(1)->SetRange(bn0[1],bn1[1]); hDstBgMC = hmcComb->Projection(2,"e"); hmcComb->GetAxis(0)->SetRange(1,nbTot[0]); hmcComb->GetAxis(1)->SetRange(1,nbTot[1]); hDstBgMC->SetName(pref+"_DistBgMC"); hDstBgMC->SetTitle(pref+" Bg. Distance from MC labels"); hDstBgMC->Scale(1./nrmDst); res->AddAtAndExpand(hDstBgMC, kBgMCDist+shift); } // // fill 1-beta matrix for (int ib0=bn0[0];ib0<=bn1[0];ib0+=bStepEta) { // eta hrawd->GetAxis(0)->SetRange(ib0,ib0+bStepEta-1); hgenb->GetAxis(0)->SetRange(ib0,ib0+bStepEta-1); for (int ib1=bn0[1];ib1<=bn1[1];ib1+=bStepZv) { // zv hrawd->GetAxis(1)->SetRange(ib1,ib1+bStepZv-1); hgenb->GetAxis(1)->SetRange(ib1,ib1+bStepZv-1); // TH1D* dstD = hrawd->Projection(2,"e"); // data "qaulity" for given eta:zv bin TH1D* dstB = hgenb->Projection(2,"e"); // data "qaulity" for given eta:zv bin double scl,sclE; TH1* nrmB = NormalizeBg(dstD,dstB,scl,sclE); // get rescaling factor for bg. from tails comparison double bgVal,bgErr; double dtVal,dtErr; // integral in the range where we look for signal Integrate(nrmB, cutSgMin, cutSgMax, bgVal, bgErr); Integrate(dstD, cutSgMin, cutSgMax, dtVal, dtErr); double beta,betaErr; GetRatE(bgVal,bgErr, dtVal, dtErr,beta,betaErr); // betaErr*=nbinstsq; // ??? RS for (int i=ib0;iSetBinContent(i,j, scl); hBgRescFc->SetBinError(i,j, sclE); } } hDstBg->Add(nrmB); delete dstD; delete dstB; delete nrmB; // } } hDstBg->Scale(1./nrmDst); // // finalize estimated bg and signal matrices hBgEst->Multiply(hBgRescFc); hSignalEst->Add(hBgEst,-1); // // finalize 1-beta for (int ib0=bn0[0];ib0<=bn1[0];ib0++) { // eta for (int ib1=bn0[1];ib1<=bn1[1];ib1++) { // zv // printf("Bin %d %d\n",ib0,ib1); double bg = hBgEst->GetBinContent(ib0,ib1); double bgE = hBgEst->GetBinError(ib0,ib1); double dt = hRawDtCut->GetBinContent(ib0,ib1); double dtE = hRawDtCut->GetBinError(ib0,ib1); double beta,betaE; GetRatE(bg,bgE,dt,dtE, beta,betaE ); h1mBeta->SetBinContent(ib0,ib1,1.-beta); h1mBeta->SetBinError(ib0,ib1,betaE); // } } // // crop to needed range for (int i=shift;iAt(i); if (!hist || !hist->InheritsFrom(TH2::Class())) continue; CropHisto(hist, bn0[0],bn1[0], bn0[1],bn1[1]); } // h1mBeta->SetMinimum(0.6); h1mBeta->SetMaximum(0.85); if (hmcComb) { h1mBetaMC->SetMinimum(0.6); h1mBetaMC->SetMaximum(0.85); } // // restore for (int i=3;i--;) { // set default axis range hrawd->GetAxis(i)->SetRange(1,nbTot[i]); hgenb->GetAxis(i)->SetRange(1,nbTot[i]); if (hmcComb) hmcComb->GetAxis(i)->SetRange(1,nbTot[i]); } // } void GetRealMinMax(TH1* histo, double &vmn, double &vmx) { TAxis *xax = histo->GetXaxis(); int nbx = xax->GetNbins(); vmn=1e6, vmx=-1e6; if (histo->InheritsFrom(TH2::Class())) { TAxis *yax = histo->GetYaxis(); int nby = yax->GetNbins(); for (int ix=nbx+2;ix--;) { for (int iy=nby+2;iy--;) { double vl = histo->GetBinContent(ix,iy); if (vlvmx) vmx = vl; } } } // else { for (int ix=nbx+2;ix--;) { double vl = histo->GetBinContent(ix); if (vlvmx) vmx = vl; } } // }