Update repo: committing code updates for Lc->pKpi analysis

ComputeEffAccWeightedAvg.C

@@ -0,0 +1,92 @@
+//___________________________________________________________________________________//
+// Macro for computing the weighted average of Lc resonant channels                  //
+// Main Function: ComputeEffAccWeightedAvg                                           //
+//___________________________________________________________________________________//
+
+#if !defined (__CINT__) || defined (__CLING__)
+
+#include <string>
+#include <vector>
+
+#include "yaml-cpp/yaml.h"
+
+#include "TROOT.h"
+#include "Riostream.h"
+#include "TH1.h"
+#include "TF1.h"
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TFile.h"
+#include "TGaxis.h"
+#include "TDatabasePDG.h"
+#include "TMath.h"
+#include "AliHFInvMassFitter.h"
+#include "AliVertexingHFUtils.h"
+
+#endif
+
+using namespace std;
+
+int ComputeEffAccWeightedAvg(TString effdir, TString particle, TString cutset, TString outFileName){
+
+    if(particle=="LctopKpi"){
+        TString f[5] = {
+            Form("%s/Eff_times_Acc_%s_NonRes%s.root", effdir.Data(), particle.Data(), cutset.Data()), //NonResonant
+            Form("%s/Eff_times_Acc_%s_KStar%s.root", effdir.Data(), particle.Data(), cutset.Data()), //KStar
+            Form("%s/Eff_times_Acc_%s_Delta%s.root", effdir.Data(), particle.Data(), cutset.Data()), //Delta
+            Form("%s/Eff_times_Acc_%s_Lambda1520%s.root", effdir.Data(), particle.Data(), cutset.Data()), //Lambda1520
+        };
+
+        TFile *file[4];
+        for (int i = 0; i < 4; i++){
+            file[i] = TFile::Open(f[i].Data());
+        }
+
+        Double_t BR[4] = {
+            //6.28 * 1e-02,        //all
+            3.5 * 1e-02,         //non res
+            1.96 * 0.667 * 1e-02, //Kstar
+            1.08 * 1e-02,        //Delta
+            2.2 * 0.225 * 1e-02 //L1520
+        };
+
+        TH1D *effC[4];
+        TH1D *effB[4];
+        for(int i = 0; i < 4; i++){
+            effC[i] = (TH1D *)file[i]->Get("hAccEffPrompt");
+            effB[i] = (TH1D *)file[i]->Get("hAccEffFD");
+            cout << " n bins " << effC[i]->GetNbinsX() << endl;
+        }
+        TH1D *effCw = (TH1D *)effC[0]->Clone("hAccEffPrompt");
+        TH1D *effBw = (TH1D *)effB[0]->Clone("hAccEffFD");
+        effCw->Reset("icse"); 
+        effBw->Reset("icse");
+
+
+        for(int iPt=0; iPt<effC[0]->GetNbinsX(); iPt++) {
+            double weightEffC = 0., weightUncEffC = 0., weightEffB = 0., weightUncEffB = 0., sumOfBR = 0.;
+            for(int iCh=0; iCh<4; iCh++) {
+                weightEffC += effC[iCh]->GetBinContent(iPt + 1) * BR[iCh];
+                weightUncEffC += effC[iCh]->GetBinError(iPt + 1) * effC[iCh]->GetBinError(iPt + 1) * BR[iCh] * BR[iCh];
+                weightEffB += effB[iCh]->GetBinContent(iPt + 1) * BR[iCh];
+                weightUncEffB += effB[iCh]->GetBinError(iPt + 1) * effB[iCh]->GetBinError(iPt + 1) * BR[iCh] * BR[iCh];
+                sumOfBR += BR[iCh];
+            }
+            weightEffC /= sumOfBR;
+            weightUncEffC = TMath::Sqrt(weightUncEffC) / sumOfBR;
+            weightEffB /= sumOfBR;
+            weightUncEffB = TMath::Sqrt(weightUncEffB) / sumOfBR;
+            effC[iCh]->SetBinContent(iPt + 1, weightEffC);
+            effC[iCh]->SetBinError(iPt + 1, weightUncEffC);
+            effB[iCh]->SetBinContent(iPt + 1, weightEffB);
+            effB[iCh]->SetBinError(iPt + 1, weightUncEffB);
+        }
+
+        TFile outFile(Form("%s/%s", effdir.Data(),outFileName.Data()),"recreate");
+        effCw->Write("hAccEffPrompt");
+        effBw->Write("hAccEffFD");
+    }
+
+    return 0;
+
+}

ComputeEffAccWeightedAvg.py

@@ -0,0 +1,94 @@
+'''
+Script for the combination of efficiency and acceptance factors
+run: python ComputeEffAccWeightedAvg.py effFileAll.root effFileNonRes.root effFileKStar.root effFileDelta.root effFileLambda1520.root outFileName.root
+'''
+
+import sys
+import argparse
+import ctypes
+import numpy as np
+from ROOT import gROOT, TFile, TCanvas, TLegend  # pylint: disable=import-error,no-name-in-module
+from ROOT import kBlack, kRed, kAzure  # pylint: disable=import-error,no-name-in-module
+from ROOT import kFullDiamond, kOpenDiamond, kFullTriangleUp, kOpenTriangleUp, kFullStar, kOpenStar, kFullCircle, kOpenCircle, kFullSquare, kOpenSquare, kFullCross, kOpenCross # pylint: disable=import-error,no-name-in-module
+from utils.StyleFormatter import SetGlobalStyle, SetObjectStyle
+
+parser = argparse.ArgumentParser(description='Arguments')
+parser.add_argument('effFileNonRes', metavar='text', default='')
+parser.add_argument('effFileKStar', metavar='text', default='')
+parser.add_argument('effFileDelta', metavar='text', default='')
+parser.add_argument('effFileLambda1520', metavar='text', default='')
+parser.add_argument('outFileName', metavar='text', default='')
+parser.add_argument("--batch", help="suppress video output", action="store_true")
+args = parser.parse_args()
+
+gROOT.SetBatch(args.batch)
+SetGlobalStyle(padleftmargin=0.14, padbottommargin=0.12, titlesize=0.045, labelsize=0.04)
+
+inputFile = [args.effFileNonRes, args.effFileKStar, args.effFileDelta, args.effFileLambda1520]
+colorMarkerPrompt = [kBlack, kRed, kRed, kRed, kRed]
+colorMarkerFD = [kBlack, kAzure, kAzure, kAzure, kAzure]
+styleMarkerPrompt = [kFullDiamond, kFullSquare, kFullStar, kFullTriangleUp, kFullCross]
+styleMarkerFD = [kOpenDiamond, kOpenSquare, kOpenStar, kOpenTriangleUp, kOpenCross]
+hEffPrompt, hEffFD = [], []
+
+for iReso, fileName, in enumerate(inputFile):
+    infile = TFile.Open(fileName)
+    hEffPrompt.append(infile.Get('hAccEffPrompt'))
+    hEffFD.append(infile.Get('hAccEffFD'))
+    hEffPrompt[iReso].SetDirectory(0)
+    hEffFD[iReso].SetDirectory(0)
+    SetObjectStyle(
+        hEffPrompt[iReso],
+        color=colorMarkerPrompt[iReso],
+        markerstyle=styleMarkerPrompt[iReso]
+    )
+    SetObjectStyle(
+        hEffFD[iReso],
+        color=colorMarkerFD[iReso],
+        markerstyle=styleMarkerFD[iReso],
+        markersize=1.5, linewidh=2, linestyle=7
+    )
+
+ptMin = hEffPrompt[0].GetBinLowEdge(1)
+ptMax = hEffPrompt[0].GetBinLowEdge(hEffPrompt[0].GetNbinsX()) + hEffPrompt[0].GetBinWidth(hEffPrompt[0].GetNbinsX())
+
+cEffPrompt = TCanvas('cEffPrompt', '', 800, 800)
+cEffPrompt.DrawFrame(ptMin, 1.e-5, ptMax, 1.,
+               ';#it{p}_{T} (GeV/#it{c});(Acc #times #epsilon);')
+cEffPrompt.SetLogy()
+for histo in hEffPrompt:
+    histo.Draw('same')
+
+cEffFD = TCanvas('cEffFD', '', 800, 800)
+cEffFD.DrawFrame(ptMin, 1.e-5, ptMax, 1.,
+               ';#it{p}_{T} (GeV/#it{c});(Acc #times #epsilon);')
+cEffFD.SetLogy()
+for histo in hEffFD:
+    histo.Draw('same')
+
+outFileNamePDF = args.outFileName.replace('.root', '')
+cEffPrompt.SaveAs('%s_prompt.pdf' % (outFileNamePDF))
+cEffFD.SaveAs('%s_FD.pdf' % (outFileNamePDF))
+
+hEffCw = hEffPrompt[0].Clone("hAccEffPrompt")
+hEffBw = hEffFD[0].Clone("hAccEffFD")
+
+BR = [3.5 * 1e-02, 1.96 * 0.667 * 1e-02, 1.08 * 1e-02, 2.2 * 0.225 * 1e-02]
+nPtBins = hEffPrompt[0].GetNbinsX()
+for iPt in range(nPtBins):
+    effC, effB, uncEffC, uncEffB, sumOfW = (0. for _ in range(5))
+    for histo, br, in zip(hEffPrompt, BR):
+        effC += histo.GetBinContent(iPt+1) * br
+        effB += histo.GetBinContent(iPt+1) * br
+        uncEffC += histo.GetBinError(iPt+1)**2 * br**2
+        uncEffB += histo.GetBinError(iPt+1)**2 * br**2
+        sumOfW += br
+    hEffCw.SetBinContent(iPt+1, (effC/sumOfW))
+    hEffCw.SetBinError(iPt+1, np.sqrt(uncEffC)/sumOfW)
+    hEffBw.SetBinContent(iPt+1, (effB/sumOfW))
+    hEffBw.SetBinError(iPt+1, np.sqrt(uncEffB)/sumOfW)
+
+outFile = TFile(args.outFileName, 'recreate')
+hEffCw.Write()
+hEffBw.Write()
+outFile.Close()

GetRawYieldsDplusDs.C

@@ -28,7 +28,7 @@
 using namespace std;
 
 enum {k010, k3050, k6080, kpp5TeVPrompt, kpp5TeVFD, kpp13TeVPrompt, kpp13TeVFD};
-enum {kDplus, kD0, kDs, kLc, kDstar};
+enum {kDplus, kD0, kDs, kLctopKpi, kLctopK0s, kDstar};
 
 //__________________________________________________________________________________________________________________
 int GetRawYieldsDplusDs(int cent = k010, bool isMC = false, TString infilename = "InvMassSpectraDplus_010_PbPb2015cuts.root", TString refFileName = "Distr_D0_MC_D0_pp13TeV_FD0.root", TString cfgfilename = "Dplus/config_Dplus_Fit.yml", TString outFileName = "RawYieldsDplus_010_PbPb2015cuts.root");
@@ -68,8 +68,10 @@ int GetRawYieldsDplusDs(int cent, bool isMC, TString infilename, TString refFile
         particle=kDs;
     } else if (ParticleName=="D0"){
         particle=kD0;
-    } else if (ParticleName=="Lc"){
-        particle=kLc;
+    } else if (ParticleName=="LctopKpi"){
+        particle=kLctopKpi;
+    } else if (ParticleName=="LctopK0s"){
+        particle=kLctopK0s;
     } else if (ParticleName=="Dstar"){
         particle=kDstar;
     } else{
@@ -181,7 +183,8 @@ int GetRawYieldsDplusDs(int cent, bool isMC, TString infilename, TString refFile
     if(particle==kDplus) massaxistit = "#it{M}(K#pi#pi) (GeV/#it{c}^{2})";
     else if(particle==kD0) massaxistit = "#it{M}(K#pi) (GeV/#it{c}^{2})";
     else if(particle==kDs) massaxistit = "#it{M}(KK#pi) (GeV/#it{c}^{2})";
-    else if(particle==kLc) massaxistit = "#it{M}(pK^{0}_{s}) (GeV/#it{c}^{2})";
+    else if(particle==kLctopKpi) massaxistit = "#it{M}(pK#pi) (GeV/#it{c}^{2})";
+    else if(particle==kLctopK0s) massaxistit = "#it{M}(pK^{0}_{s}) (GeV/#it{c}^{2})";
     else if (particle==kDstar) massaxistit = "#it{M}(pi^{+}) (GeV/#it{c}^{2})";
 
     //load inv-mass histos
@@ -374,7 +377,7 @@ int GetRawYieldsDplusDs(int cent, bool isMC, TString infilename, TString refFile
     if(particle==kDplus) massForFit = massDplus;
     else if (particle==kDs) massForFit = massDs;
     else if (particle==kD0) massForFit = massD0;
-    else if (particle==kLc) massForFit = massLc;
+    else if (particle==kLctopKpi || particle==kLctopK0s) massForFit = massLc;
     else if (particle==kDstar) massForFit = dmassDstar;
 
     TH1F* hMassForFit[nPtBins];

GetRawYieldsDplusDs.py

@@ -110,7 +110,9 @@
     massAxisTit = '#it{M}(K#pi#pi) (GeV/#it{c}^{2})'
 elif particleName == 'Ds':
     massAxisTit = '#it{M}(KK#pi) (GeV/#it{c}^{2})'
-elif particleName == 'Lc':
+elif particleName == 'LctopKpi':
+    massAxisTit = '#it{M}(pK#pi) (GeV/#it{c}^{2})'
+elif particleName == 'LctopK0s':
     massAxisTit = '#it{M}(pK^{0}_{s}) (GeV/#it{c}^{2})'
 elif particleName == 'Dstar':
     massAxisTit = '#it{M}(K#pi#pi) - #it{M}(K#pi) (GeV/#it{c}^{2})'

ProjectDplusDsSparse.py

@@ -38,6 +38,8 @@
 parser.add_argument('--Bspeciesweights', type=float, nargs=5, required=False,
                     help='values of weights for the different hadron species '
                          '(B0weight, Bplusweight, Bsweight, Lbweight, Otherweight)')
+parser.add_argument('--LctopKpireso', type=int, required=False, default=None,
+                    help='values to project single LctopKpi resonant channel')
 args = parser.parse_args()
 
 #TODO: add support for application of 2D weights
@@ -151,6 +153,8 @@
         binMax = sparseReco[refSparse].GetAxis(axisNum).FindBin(cutVars[iVar]['max'][iPt] * 0.9999)
         if 'RecoAll' in sparseReco:
             sparseReco['RecoAll'].GetAxis(axisNum).SetRange(binMin, binMax)
+            if args.LctopKpireso:
+                sparseReco['RecoAll'].GetAxis(2).SetRange(args.LctopKpireso+1, args.LctopKpireso+1) # "+1" applied to fix the discrepancy between the reso channel and the filled bin
         if isMC:
             sparseReco['RecoPrompt'].GetAxis(axisNum).SetRange(binMin, binMax)
             sparseReco['RecoFD'].GetAxis(axisNum).SetRange(binMin, binMax)
@@ -159,6 +163,9 @@
             if enableSecPeak:
                 sparseReco['RecoSecPeakPrompt'].GetAxis(axisNum).SetRange(binMin, binMax)
                 sparseReco['RecoSecPeakFD'].GetAxis(axisNum).SetRange(binMin, binMax)
+            if args.LctopKpireso:
+                sparseReco['RecoPrompt'].GetAxis(2).SetRange(args.LctopKpireso+1, args.LctopKpireso+1) # "+1" applied to fix the discrepancy between the reso channel and the filled bin
+                sparseReco['RecoFD'].GetAxis(3).SetRange(args.LctopKpireso+1, args.LctopKpireso+1) # "+1" applied to fix the discrepancy between the reso channel and the filled bin
 
     for iVar in ('InvMass', 'Pt'):
         varName = 'Pt' if iVar == 'Pt' else 'Mass'
@@ -298,6 +305,10 @@
         binGenMax = sparseGen['GenPrompt'].GetAxis(0).FindBin(ptMax*0.9999)
         sparseGen['GenPrompt'].GetAxis(0).SetRange(binGenMin, binGenMax)
         sparseGen['GenFD'].GetAxis(0).SetRange(binGenMin, binGenMax)
+        if args.LctopKpireso:
+            sparseGen['GenPrompt'].GetAxis(2).SetRange(args.LctopKpireso+1, args.LctopKpireso+1) # "+1" applied to fix the discrepancy between the reso channel and the filled bin
+            sparseGen['GenFD'].GetAxis(3).SetRange(args.LctopKpireso+1, args.LctopKpireso+1) # "+1" applied to fix the discrepancy between the reso channel and the filled bin
+
         if not args.multweights:
             hGenPtPrompt = sparseGen['GenPrompt'].Projection(0)
             # apply pt weights

ProjectDplusDsTree.py

@@ -37,6 +37,8 @@
                     help='First path of the pT weights file, second name of the pT weights histogram')
 parser.add_argument('--multweights', metavar=('text', 'text'), nargs=2, required=False,
                     help='First path of the mult weights file, second name of the mult weights histogram')
+parser.add_argument('--LctopKpireso', type=int, required=False,
+                    help='values to project single LctopKpi resonant channel (1: NonRes; 2: KStar; 3: Delta; 4: Lambda1520)')
 parser.add_argument('--std', help='adapt to std. analysis cuts', action='store_true')
 args = parser.parse_args()
 
@@ -65,6 +67,9 @@
 bitPrompt = 2
 bitFD = 3
 bitRefl = 4
+#define specific bits for Lc->pKpi resonant channel
+bitsLcResonance = [9, 10, 11, 12] # ['NonRes', 'Lambda1520', 'KStar', 'Delta']
+
 # define mass binning
 particle = inputCfg['tree']['particle']
 if particle == 'Ds':
@@ -147,13 +152,17 @@
     if 'cand_type' in dataFramePrompt.columns: #if not filtered tree, select only prompt and not reflected
         dataFramePrompt = FilterBitDf(dataFramePrompt, 'cand_type', [bitSignal, bitPrompt], 'and')
         dataFramePrompt = FilterBitDf(dataFramePrompt, 'cand_type', [bitRefl], 'not')
+        if args.LctopKpireso in range(1, 5):
+            dataFramePrompt = FilterBitDf(dataFramePrompt, 'cand_type', [bitsLcResonance[args.LctopKpireso]], 'and')
     dataFramePrompt.reset_index(inplace=True)
 
     dataFrameFD = LoadDfFromRootOrParquet(inputCfg['tree']['filenameFD'], inputCfg['tree']['dirname'],
                                           inputCfg['tree']['treename'])
     if 'cand_type' in dataFrameFD.columns: #if not filtered tree, select only FD and not reflected
         dataFrameFD = FilterBitDf(dataFrameFD, 'cand_type', [bitSignal, bitFD], 'and')
         dataFrameFD = FilterBitDf(dataFrameFD, 'cand_type', [bitRefl], 'not')
+        if args.LctopKpireso in range(1, 5):
+            dataFrameFD = FilterBitDf(dataFrameFD, 'cand_type', [bitsLcResonance[args.LctopKpireso]], 'and')
     dataFrameFD.reset_index(inplace=True)
 
     # compute pt weights
@@ -200,6 +209,9 @@
         binGenMax = sparseGen['GenPrompt'].GetAxis(0).FindBin(ptMax * 0.9999)
         sparseGen['GenPrompt'].GetAxis(0).SetRange(binGenMin, binGenMax)
         sparseGen['GenFD'].GetAxis(0).SetRange(binGenMin, binGenMax)
+        if args.LctopKpireso:
+            sparseGen['GenPrompt'].GetAxis(2).SetRange(args.LctopKpireso+1, args.LctopKpireso+1) # "+1" applied to fix the discrepancy between the reso channel and the filled bin
+            sparseGen['GenFD'].GetAxis(3).SetRange(args.LctopKpireso+1, args.LctopKpireso+1) # "+1" applied to fix the discrepancy between the reso channel and the filled bin
 
         if args.multweights:
             hMultVsGenPtPrompt = sparseGen['GenPrompt'].Projection(4, 0)