Merge pull request #110 from JeffersonLab/pid_syst

Pid syst

AnalysisHowTo/PIDstudy/DSelector_p2gamma.h

@@ -0,0 +1,64 @@
+#ifndef DSelector_p2gamma_h
+#define DSelector_p2gamma_h
+
+#include <iostream>
+
+#include "DSelector/DSelector.h"
+#include "DSelector/DHistogramActions.h"
+#include "DSelector/DCutActions.h"
+
+#include "TH1I.h"
+#include "TH2I.h"
+
+class DSelector_p2gamma : public DSelector
+{
+	public:
+
+		DSelector_p2gamma(TTree* locTree = NULL) : DSelector(locTree){}
+		virtual ~DSelector_p2gamma(){}
+
+		void Init(TTree *tree);
+		Bool_t Process(Long64_t entry);
+
+	private:
+
+		void Get_ComboWrappers(void);
+		void Finalize(void);
+
+		// BEAM POLARIZATION INFORMATION
+		UInt_t dPreviousRunNumber;
+		bool dIsPolarizedFlag; //else is AMO
+		bool dIsPARAFlag; //else is PERP or AMO
+
+		// ANALYZE CUT ACTIONS
+		// // Automatically makes mass histograms where one cut is missing
+		DHistogramAction_AnalyzeCutActions* dAnalyzeCutActions;
+
+		//CREATE REACTION-SPECIFIC PARTICLE ARRAYS
+
+		//Step 0
+		DParticleComboStep* dStep0Wrapper;
+		DBeamParticle* dComboBeamWrapper;
+		DNeutralParticleHypothesis* dPhoton1Wrapper;
+		DNeutralParticleHypothesis* dPhoton2Wrapper;
+		DChargedTrackHypothesis* dProtonWrapper;
+
+		// DEFINE YOUR HISTOGRAMS HERE
+		// EXAMPLES:
+		TH1I* dHist_MissingMassSquared;
+		TH1I* dHist_BeamEnergy;
+
+	ClassDef(DSelector_p2gamma, 0);
+};
+
+void DSelector_p2gamma::Get_ComboWrappers(void)
+{
+	//Step 0
+	dStep0Wrapper = dComboWrapper->Get_ParticleComboStep(0);
+	dComboBeamWrapper = static_cast<DBeamParticle*>(dStep0Wrapper->Get_InitialParticle());
+	dPhoton1Wrapper = static_cast<DNeutralParticleHypothesis*>(dStep0Wrapper->Get_FinalParticle(0));
+	dPhoton2Wrapper = static_cast<DNeutralParticleHypothesis*>(dStep0Wrapper->Get_FinalParticle(1));
+	dProtonWrapper = static_cast<DChargedTrackHypothesis*>(dStep0Wrapper->Get_FinalParticle(2));
+}
+
+#endif // DSelector_p2gamma_h

AnalysisHowTo/PIDstudy/README

@@ -0,0 +1,17 @@
+# Example for PID cut study: Analysis How To's (JRS 6/1/20): jrsteven@jlab.org
+
+# This example shows how to use the AnalyzeCutAction tool to fill mass histograms for pi0 -> gg under differing PID requirements.  The yields can be extracted to show how they depend on the individual cuts.  Python scripts display the standard PID plots and make some data/MC comparisons
+
+# Step 0) Include your analysis channel in the launch with loosened PID cuts (https://halldweb.jlab.org/wiki-private/index.php/Spring_2017_Analysis_Launch#Version38), or run the analysis yourself over the REST files with those conditions.
+
+# Step 1) Produce Analysis TTrees for your MC sample under the same conditions as the data in Step 0): e.g. using the script run_MC_pid_syst.csh in this How To.  Note: this may require batch jobs for large samples.
+
+# Step 2) Open your reaction's ROOT tree and run the DSelector using the AnalyzeCutAction to compare the PID cuts (do this for data and MC).  
+root -l -b -q runSelector.C
+
+# Step 3) Use the python script to make standard PID plots (prior to default cuts) for all particles in your reaction.  Note: there are options to change the momentum slice for 1D plots, etc. in the script.  Feel free to change them or send requests for other options.
+python -b plotPIDdataMC.py 
+
+# Step 4) Use the python script to extract particle yields for each PID cut being applied.  Note: You'll need to define a fit range, background function, etc. for your specific reaction, but the procedure for looping over the different PID cut conditions should be the same for all reactions.
+python -b plotPIDyields.py
+

AnalysisHowTo/PIDstudy/plotPIDdataMC.py

@@ -0,0 +1,134 @@
+import os
+from ROOT import gDirectory,gStyle,TFile,TCanvas,TF1,TLegend
+
+# define function to get list of string keys in a directory
+def GetKeyNames( self, dir = "" ):
+    self.cd(dir)
+    return [key.GetName() for key in gDirectory.GetListOfKeys()]
+
+TFile.GetKeyNames = GetKeyNames
+
+# some basic setup for canvases and legends
+gStyle.SetOptStat(0)
+leg = TLegend(0.6,0.7,0.9,0.9)
+
+cc2D = TCanvas("cc2D","cc2D",800,400)
+cc2D.Divide(2,1)
+cc1D = TCanvas("cc1D","cc1D",600,400)
+
+plotDir = "plots/"
+if not os.path.exists(plotDir):
+    os.mkdir(plotDir)
+
+# default dE/dx cut for plotting on 2D distribution
+fMinProton_dEdx = TF1("fMinProton_dEdx", "exp(-1.*[0]*x + [1]) + [2]", 0., 10.) # cut for dEdx curve
+fMinProton_dEdx.SetParameters(5.2, 3.0, 0.9)
+fMinProton_dEdx.SetLineColor(2)
+
+
+# Parser needs
+## Files name and label from text input (defines canvas size and labels)
+## Momentum range for given particle?
+## DeltaT range to zoom in
+
+pidPath = "Hist_ParticleID_Initial"
+maxDeltaT = 2.5
+minP = 0.0
+maxP = 2.5
+minSliceP = 0.8
+maxSliceP = 1.0
+
+# Insert your input files here for Data and MC
+files = []
+files.append(TFile.Open("hist_sum_30730_30788.root")) # Data
+files.append(TFile.Open("hist_sum_30274_31057_sim_g4.root")) #MC
+
+# Open file and get list of keys
+f = files[0]
+stepKeyList = f.GetKeyNames(pidPath)
+
+# collect list of histogram names (full path)
+hists = []
+particles = []
+for step in stepKeyList:
+    print("Keys in file:", step)
+    particleKeyList = f.GetKeyNames(pidPath+"/"+step)
+
+    for particle in particleKeyList:
+        print("Particle in step:", particle)
+        histKeyList = f.GetKeyNames(pidPath+"/"+step+"/"+particle)
+
+        for hist in histKeyList:
+            print("Histogram for particle:", hist)
+            if "DIRC" in hist or "Beta" in hist: # skip DIRC and Best histograms for now
+                continue
+
+            particles.append(particle)
+            hists.append(pidPath+"/"+step+"/"+particle+"/"+hist)
+
+# make plots for each histogram in the list
+ihist = 0
+for hist,particle in zip(hists,particles):
+    print(hist)
+    ifile = 0
+    projMinBin = 0
+    projMaxBin = 999
+    norm = 0
+
+
+    for file in files:
+        h = file.Get(hist)
+
+        if "DeltaT" in hist:
+            h.GetYaxis().SetRangeUser(-1.*maxDeltaT,maxDeltaT)
+
+        # make 1D slices in momentum for comparison
+        plot1D = "p (GeV/c)" in h.GetXaxis().GetTitle()
+        if plot1D:
+
+            # 1D projection
+            cc1D.cd()
+            if ifile == 0: #Data specific values
+                projMinBin = h.ProjectionX().FindBin(minSliceP);
+                projMaxBin = h.ProjectionX().FindBin(maxSliceP);
+                h1D = h.ProjectionY("%s" % hist,projMinBin,projMaxBin)
+                h1D.SetTitle(h.GetTitle())
+                h1D.SetLineColor(ifile+1)
+                h1D.SetMarkerColor(ifile+1)
+                h1D.SetMarkerStyle(20)
+                norm = h1D.GetMaximum()
+                h1D.Draw("pe")
+
+                if ihist == 0:
+                    leg.AddEntry(h1D,"Data","p")
+
+            else: # MC specific values
+                h1D = h.ProjectionY("%d" % ifile,projMinBin,projMaxBin)
+                if h1D.GetMaximum() > 0:
+                    h1D.Scale(norm/h1D.GetMaximum())
+                h1D.SetLineColor(ifile+1)
+                h1D.SetMarkerColor(ifile+1)
+                h1D.SetMarkerStyle(20)
+                h1D.Draw("pe same")
+
+                if ihist == 0:
+                    leg.AddEntry(h1D,"MC","p")
+
+                leg.Draw("same")
+
+            if "Proton" in particle:
+                h.GetXaxis().SetRangeUser(minP,maxP)
+
+        cc2D.cd(ifile+1)
+        h.Draw("colz")
+
+        if "CDC dE/dx" in h.GetYaxis().GetTitle():
+            fMinProton_dEdx.Draw("same")
+
+        ifile += 1
+    ihist += 1
+
+    # print plots
+    cc2D.Print("%s%s.pdf" % (plotDir,hist.replace("/","_")))
+    if plot1D:
+        cc1D.Print("%spSlice_%s.pdf" % (plotDir,hist.replace("/","_")))

AnalysisHowTo/PIDstudy/plotPIDyields.py

@@ -0,0 +1,129 @@
+import os
+from ROOT import gDirectory,gPad,TFile,TCanvas,TF1,TH1F,TMath
+
+# define function to get list of string keys in a directory
+def GetKeyNames( self, dir = "" ):
+    self.cd(dir)
+    return [key.GetName() for key in gDirectory.GetListOfKeys()]
+
+TFile.GetKeyNames = GetKeyNames
+
+# Canvas for plotting results
+cc1D = TCanvas("cc1D","cc1D",800,400)
+cc1D.Divide(2,1)
+
+plotDir = "plots/"
+if not os.path.exists(plotDir):
+    os.mkdir(plotDir)
+
+# Define double gaussian fit function
+fPi0 = TF1("fitPi0", "gaus(0) + gaus(3) + pol1(6)", 0.08, 0.19)
+fPi0.SetLineColor(2)
+
+# Parser needs
+## Files name and label from text input (defines canvas size and labels)
+
+pidPath = "Hist_AnalyzeCutActions_CutActionEffect"
+
+# Insert your input files here for Data and MC
+files = []
+files.append(TFile.Open("hist_sum_30730_30788.root")) # Data
+files.append(TFile.Open("hist_sum_30274_31057_sim_g4.root")) # MC
+
+# Open data file file and get list of keys in AnalyzeCutActions directory
+f = files[0]
+cuts = f.GetKeyNames(pidPath)
+hists = []
+
+# collect list of histogram names (full path)
+for cut in cuts:
+    print("Cuts in file:", cut)
+    hists.append(pidPath+"/"+cut)
+#print(cuts)
+#print(hists)
+
+nominalYieldData = 0
+nominalYieldMC = 0
+
+# crate histograms to store yields
+nEmptyBins = 4
+yieldHistData = TH1F("yieldData",";", len(cuts)+nEmptyBins, -0.5, -0.5+len(cuts) +nEmptyBins)
+yieldHistMC = TH1F("yieldMC",";", len(cuts)+nEmptyBins, -0.5, -0.5+len(cuts) +nEmptyBins)
+for cut in cuts: # set histogram bin names
+    yieldHistData.Fill(cut, 0)
+    yieldHistMC.Fill(cut, 0)
+yieldRatioData = yieldHistData.Clone("yieldRatioData")
+yieldRatioMC = yieldHistMC.Clone("yieldRatioMC")
+
+# make plots for each histogram in the list
+ihist = 0
+for hist,cut in zip(hists,cuts):
+    ifile = 0
+
+    # fit yield in data and MC for each histogram
+    for file in files:
+        cc1D.cd(ifile+1)
+        h = file.Get(hist)
+
+        # setup simple histogram fit with double gaussian and polynomial
+        fPi0.SetParameters(h.GetMaximum(), 0.135, 0.007, h.GetMaximum()/10., 0.135, 0.015)
+
+        # draw histogram and fit
+        h.SetMarkerStyle(20)
+        h.Draw("pe")
+        fitResult = h.Fit(fPi0,"Q","",0.08,0.19)
+        fitYield = fPi0.Integral(0.11,0.16)/h.GetXaxis().GetBinWidth(1)
+        fitYieldError = fPi0.IntegralError(0.11,0.16)/h.GetXaxis().GetBinWidth(1) #,fitResult.GetParams(),fitResult.GetCovarianceMatrix().GetMatrixArray(),0.05
+
+        if ifile==0: #Data specific values
+            if ihist == 0:
+                nominalYieldData = fitYield
+                nominalYieldErrorData = fitYieldError
+            yieldHistData.SetBinContent(ihist+1,fitYield)
+            yieldHistData.SetBinError(ihist+1,fitYieldError)
+            yieldRatio = nominalYieldData/fitYield
+            yieldRatioRelError = TMath.Sqrt(pow(nominalYieldErrorData/nominalYieldData,2) + pow(fitYieldError/fitYield,2))
+            yieldRatioData.SetBinContent(ihist+1,nominalYieldData/fitYield)
+            yieldRatioData.SetBinError(ihist+1,yieldRatio*yieldRatioRelError)
+        else: # MC specific values
+            if ihist == 0:
+                nominalYieldMC = fitYield
+                nominalYieldErrorMC = fitYieldError
+            yieldHistMC.SetBinContent(ihist+1,fitYield)
+            yieldHistMC.SetBinError(ihist+1,fitYieldError)
+            yieldRatio = nominalYieldMC/fitYield
+            yieldRatioRelError = TMath.Sqrt(pow(nominalYieldErrorMC/nominalYieldMC,2) + pow(fitYieldError/fitYield,2))
+            yieldRatioMC.SetBinContent(ihist+1,yieldRatio)
+            yieldRatioMC.SetBinError(ihist+1,yieldRatio*yieldRatioRelError)
+        ifile += 1
+    ihist += 1
+
+    cc1D.Print("%sfit_%s.pdf" % (plotDir,cut.replace("/","_")))
+    cc1D.Print("%sfitSummary.pdf(" % plotDir)
+
+# plot yields and take ratio of data/MC efficiency ratio
+ccYieldRatio = TCanvas("ccYieldRatio","ccYieldRatio",900,400)
+ccYieldRatio.Divide(3,1)
+
+ccYieldRatio.cd(1);
+gPad.SetBottomMargin(0.15)
+yieldHistData.SetMarkerStyle(20)
+yieldHistData.Draw("pe")
+ccYieldRatio.cd(2);
+gPad.SetBottomMargin(0.15)
+yieldHistMC.SetMarkerStyle(20)
+yieldHistMC.Draw("pe")
+
+ccYieldRatio.cd(3)
+gPad.SetBottomMargin(0.15)
+gPad.SetLeftMargin(0.15)
+efficRatio = yieldRatioData.Clone("efficRatio")
+efficRatio.GetYaxis().SetTitle("Efficiency Ratio Data/MC")
+efficRatio.Divide(yieldRatioMC)
+efficRatio.SetMinimum(0.9)
+efficRatio.SetMaximum(1.1)
+efficRatio.SetMarkerStyle(20)
+efficRatio.Draw("pe")
+
+ccYieldRatio.Print("%sefficiency.pdf" % plotDir)
+ccYieldRatio.Print("%sfitSummary.pdf)" % plotDir)