diff --git a/AnalysisHowTo/PIDstudy/DSelector_p2gamma.C b/AnalysisHowTo/PIDstudy/DSelector_p2gamma.C new file mode 100644 index 00000000..82467d09 --- /dev/null +++ b/AnalysisHowTo/PIDstudy/DSelector_p2gamma.C @@ -0,0 +1,324 @@ +#include "DSelector_p2gamma.h" + +void DSelector_p2gamma::Init(TTree *locTree) +{ + // USERS: IN THIS FUNCTION, ONLY MODIFY SECTIONS WITH A "USER" OR "EXAMPLE" LABEL. LEAVE THE REST ALONE. + + // The Init() function is called when the selector needs to initialize a new tree or chain. + // Typically here the branch addresses and branch pointers of the tree will be set. + // Init() will be called many times when running on PROOF (once per file to be processed). + + //USERS: SET OUTPUT FILE NAME //can be overriden by user in PROOF + dOutputFileName = "p2gamma.root"; //"" for none + dOutputTreeFileName = ""; //"" for none + dFlatTreeFileName = ""; //output flat tree (one combo per tree entry), "" for none + dFlatTreeName = ""; //if blank, default name will be chosen + + //Because this function gets called for each TTree in the TChain, we must be careful: + //We need to re-initialize the tree interface & branch wrappers, but don't want to recreate histograms + bool locInitializedPriorFlag = dInitializedFlag; //save whether have been initialized previously + DSelector::Init(locTree); //This must be called to initialize wrappers for each new TTree + //gDirectory now points to the output file with name dOutputFileName (if any) + if(locInitializedPriorFlag) + return; //have already created histograms, etc. below: exit + + Get_ComboWrappers(); + dPreviousRunNumber = 0; + + /*********************************** EXAMPLE USER INITIALIZATION: ANALYSIS ACTIONS **********************************/ + + // EXAMPLE: Create deque for histogramming particle masses: + // // For histogramming the phi mass in phi -> K+ K- + // // Be sure to change this and dAnalyzeCutActions to match reaction + std::deque MyPi0; + MyPi0.push_back(Gamma); MyPi0.push_back(Gamma); + + //KINEMATICS (INITIAL) + dAnalysisActions.push_back(new DHistogramAction_ParticleComboKinematics(dComboWrapper, false, "Initial")); + + //PID (INITIAL) + dAnalysisActions.push_back(new DHistogramAction_ParticleID(dComboWrapper, false, "Initial")); + + //MASSES + dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, false, 0, MyPi0, 110, 0.08, 0.19, "Pi0")); + dAnalysisActions.push_back(new DHistogramAction_MissingMassSquared(dComboWrapper, false, 1000, -0.25, 0.25)); + //KINFIT RESULTS + dAnalysisActions.push_back(new DHistogramAction_KinFitResults(dComboWrapper)); + + // ************* CROSS SECTION SYSTEMATIC CUTS ****************** // + + //CUT MISSING MASS + dAnalysisActions.push_back(new DCutAction_MissingMassSquared(dComboWrapper, false, -0.01, 0.01)); + + //NOMINAL PID TIMING CUTS (https://halldweb.jlab.org/wiki/index.php/PID_study_proposal) + dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 1.5, Gamma, SYS_BCAL)); + dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 1.5, Gamma, SYS_FCAL)); + dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 1.0, Proton, SYS_BCAL)); + dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 0.6, Proton, SYS_TOF)); + dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 2.0, Proton, SYS_FCAL)); + dAnalysisActions.push_back(new DCutAction_PIDDeltaT(dComboWrapper, false, 2.5, Proton, SYS_START)); + + //MASSES (FINAL) + dAnalysisActions.push_back(new DHistogramAction_InvariantMass(dComboWrapper, false, 0, MyPi0, 110, 0.08, 0.19, "Pi0Final")); + + //KINEMATICS (FINAL) + dAnalysisActions.push_back(new DHistogramAction_ParticleComboKinematics(dComboWrapper, false, "Final")); + + //PID (FINAL) + dAnalysisActions.push_back(new DHistogramAction_ParticleID(dComboWrapper, false, "Final")); + + // ANALYZE CUT ACTIONS + // Change MyPi0 to match reaction + dAnalyzeCutActions = new DHistogramAction_AnalyzeCutActions( dAnalysisActions, dComboWrapper, false, 0, MyPi0, 110, 0.08, 0.19, "CutActionEffect" ); + + //INITIALIZE ACTIONS + //If you create any actions that you want to run manually (i.e. don't add to dAnalysisActions), be sure to initialize them here as well + Initialize_Actions(); + dAnalyzeCutActions->Initialize(); // manual action, must call Initialize() + + /******************************** EXAMPLE USER INITIALIZATION: STAND-ALONE HISTOGRAMS *******************************/ + + //EXAMPLE MANUAL HISTOGRAMS: + dHist_MissingMassSquared = new TH1I("MissingMassSquared", ";Missing Mass Squared (GeV/c^{2})^{2}", 600, -0.06, 0.06); + dHist_BeamEnergy = new TH1I("BeamEnergy", ";Beam Energy (GeV)", 600, 0.0, 12.0); +} + +Bool_t DSelector_p2gamma::Process(Long64_t locEntry) +{ + // The Process() function is called for each entry in the tree. The entry argument + // specifies which entry in the currently loaded tree is to be processed. + // + // This function should contain the "body" of the analysis. It can contain + // simple or elaborate selection criteria, run algorithms on the data + // of the event and typically fill histograms. + // + // The processing can be stopped by calling Abort(). + // Use fStatus to set the return value of TTree::Process(). + // The return value is currently not used. + + //CALL THIS FIRST + DSelector::Process(locEntry); //Gets the data from the tree for the entry + //cout << "RUN " << Get_RunNumber() << ", EVENT " << Get_EventNumber() << endl; + //TLorentzVector locProductionX4 = Get_X4_Production(); + + /******************************************** GET POLARIZATION ORIENTATION ******************************************/ + + //Only if the run number changes + //RCDB environment must be setup in order for this to work! (Will return false otherwise) + UInt_t locRunNumber = Get_RunNumber(); + if(locRunNumber != dPreviousRunNumber) + { + dIsPolarizedFlag = dAnalysisUtilities.Get_IsPolarizedBeam(locRunNumber, dIsPARAFlag); + dPreviousRunNumber = locRunNumber; + } + + /********************************************* SETUP UNIQUENESS TRACKING ********************************************/ + + //ANALYSIS ACTIONS: Reset uniqueness tracking for each action + //For any actions that you are executing manually, be sure to call Reset_NewEvent() on them here + Reset_Actions_NewEvent(); + dAnalyzeCutActions->Reset_NewEvent(); // manual action, must call Reset_NewEvent() + + //PREVENT-DOUBLE COUNTING WHEN HISTOGRAMMING + //Sometimes, some content is the exact same between one combo and the next + //e.g. maybe two combos have different beam particles, but the same data for the final-state + //When histogramming, you don't want to double-count when this happens: artificially inflates your signal (or background) + //So, for each quantity you histogram, keep track of what particles you used (for a given combo) + //Then for each combo, just compare to what you used before, and make sure it's unique + + //EXAMPLE 1: Particle-specific info: + set locUsedSoFar_BeamEnergy; //Int_t: Unique ID for beam particles. set: easy to use, fast to search + + //EXAMPLE 2: Combo-specific info: + //In general: Could have multiple particles with the same PID: Use a set of Int_t's + //In general: Multiple PIDs, so multiple sets: Contain within a map + //Multiple combos: Contain maps within a set (easier, faster to search) + set > > locUsedSoFar_MissingMass; + + //INSERT USER ANALYSIS UNIQUENESS TRACKING HERE + + /************************************************* LOOP OVER COMBOS *************************************************/ + + //Loop over combos + for(UInt_t loc_i = 0; loc_i < Get_NumCombos(); ++loc_i) + { + //Set branch array indices for combo and all combo particles + dComboWrapper->Set_ComboIndex(loc_i); + + // Is used to indicate when combos have been cut + if(dComboWrapper->Get_IsComboCut()) // Is false when tree originally created + continue; // Combo has been cut previously + + /********************************************** GET PARTICLE INDICES *********************************************/ + + //Used for tracking uniqueness when filling histograms, and for determining unused particles + + //Step 0 + Int_t locBeamID = dComboBeamWrapper->Get_BeamID(); + Int_t locPhoton1NeutralID = dPhoton1Wrapper->Get_NeutralID(); + Int_t locPhoton2NeutralID = dPhoton2Wrapper->Get_NeutralID(); + Int_t locProtonTrackID = dProtonWrapper->Get_TrackID(); + + /*********************************************** GET FOUR-MOMENTUM **********************************************/ + + // Get P4's: //is kinfit if kinfit performed, else is measured + //dTargetP4 is target p4 + //Step 0 + TLorentzVector locBeamP4 = dComboBeamWrapper->Get_P4(); + TLorentzVector locPhoton1P4 = dPhoton1Wrapper->Get_P4(); + TLorentzVector locPhoton2P4 = dPhoton2Wrapper->Get_P4(); + TLorentzVector locProtonP4 = dProtonWrapper->Get_P4(); + + // Get Measured P4's: + //Step 0 + TLorentzVector locBeamP4_Measured = dComboBeamWrapper->Get_P4_Measured(); + TLorentzVector locPhoton1P4_Measured = dPhoton1Wrapper->Get_P4_Measured(); + TLorentzVector locPhoton2P4_Measured = dPhoton2Wrapper->Get_P4_Measured(); + TLorentzVector locProtonP4_Measured = dProtonWrapper->Get_P4_Measured(); + + /********************************************* GET COMBO RF TIMING INFO *****************************************/ + + TLorentzVector locBeamX4_Measured = dComboBeamWrapper->Get_X4_Measured(); + Double_t locBunchPeriod = dAnalysisUtilities.Get_BeamBunchPeriod(Get_RunNumber()); + Double_t locDeltaT_RF = dAnalysisUtilities.Get_DeltaT_RF(Get_RunNumber(), locBeamX4_Measured, dComboWrapper); + Int_t locRelBeamBucket = dAnalysisUtilities.Get_RelativeBeamBucket(Get_RunNumber(), locBeamX4_Measured, dComboWrapper); // 0 for in-time events, non-zero integer for out-of-time photons + Int_t locNumOutOfTimeBunchesInTree = 4; //YOU need to specify this number + //Number of out-of-time beam bunches in tree (on a single side, so that total number out-of-time bunches accepted is 2 times this number for left + right bunches) + + Bool_t locSkipNearestOutOfTimeBunch = true; // True: skip events from nearest out-of-time bunch on either side (recommended). + Int_t locNumOutOfTimeBunchesToUse = locSkipNearestOutOfTimeBunch ? locNumOutOfTimeBunchesInTree-1:locNumOutOfTimeBunchesInTree; + Double_t locAccidentalScalingFactor = dAnalysisUtilities.Get_AccidentalScalingFactor(Get_RunNumber(), locBeamP4.E()); // Ideal value would be 1, but deviations observed: need added factor. + Double_t locAccidentalScalingFactorError = dAnalysisUtilities.Get_AccidentalScalingFactorError(Get_RunNumber(), locBeamP4.E()); // Ideal value would be 1, but deviations observed, need added factor. + Double_t locHistAccidWeightFactor = locRelBeamBucket==0 ? 1 : -locAccidentalScalingFactor/(2*locNumOutOfTimeBunchesToUse) ; // Weight by 1 for in-time events, ScalingFactor*(1/NBunches) for out-of-time + if(locSkipNearestOutOfTimeBunch && abs(locRelBeamBucket)==1) continue; // Skip nearest out-of-time bunch: tails of in-time distribution also leak in + + /********************************************* COMBINE FOUR-MOMENTUM ********************************************/ + + // DO YOUR STUFF HERE + + // Combine 4-vectors + TLorentzVector locMissingP4_Measured = locBeamP4_Measured + dTargetP4; + locMissingP4_Measured -= locPhoton1P4_Measured + locPhoton2P4_Measured + locProtonP4_Measured; + + TLorentzVector locPhotonPairP4_Measured = locPhoton1P4_Measured + locPhoton2P4_Measured; + + // Combine 4-vectors for Missing momentum for gamma p -> Xp + TLorentzVector locMissingP4_OffProton = locBeamP4_Measured + dTargetP4; + locMissingP4_OffProton -= locProtonP4_Measured; + + /******************************************** EXECUTE YOUR NON-PID SELECTION CUTS HERE **********************************************/ + + /**************************************** PreSelection ******************************************/ + //The space-time information when the photons hit the detector + TLorentzVector locPhoton1_X4_Measured = dPhoton1Wrapper->Get_X4_Shower(); + TLorentzVector locPhoton2_X4_Measured = dPhoton2Wrapper->Get_X4_Shower(); + + //The space-time information when the photons hit the detector + TLorentzVector X4_Photon1 = locPhoton1_X4_Measured; + double x_g1 = X4_Photon1.X(); + double y_g1 = X4_Photon1.Y(); + double z_g1 = X4_Photon1.Z(); + double r_g1 = X4_Photon1.Perp(); + double t_g1 = X4_Photon1.T(); + + TLorentzVector X4_Photon2 = locPhoton2_X4_Measured; + double x_g2 = X4_Photon2.X(); + double y_g2 = X4_Photon2.Y(); + double z_g2 = X4_Photon2.Z(); + double r_g2 = X4_Photon2.Perp(); + double t_g2 = X4_Photon2.T(); + + //the distant between two photons in XY plane, used for splits cut in FCAL + double dis_g1_g2 = sqrt((x_g1-x_g2)*(x_g1-x_g2)+(y_g1-y_g2)*(y_g1-y_g2)); + + //define the detector photon angle + TVector3 X3_Photon1 = X4_Photon1.Vect(); + TVector3 X3_Photon2 = X4_Photon2.Vect(); + TVector3 X3_target(0.0, 0.0, 65.0); + + double theta_det_g1 = (X3_Photon1-X3_target).Theta(); + double theta_det_g2 = (X3_Photon2-X3_target).Theta(); + double theta_det_gmin = theta_det_g1; + if(theta_det_g1>theta_det_g2) + theta_det_gmin = theta_det_g2; + double theta_det_g_min = theta_det_gmin*180/TMath::Pi(); + + //vertex information + TLorentzVector locProtonX4_Measured = dProtonWrapper->Get_X4_Measured(); + double z = 0.0; + double r = 0.0; + z = locProtonX4_Measured.Z(); + r = locProtonX4_Measured.Perp(); + + if(locProtonP4_Measured.Vect().Mag() < 0.25) //cut in momentum of proton + continue; + + if(z<=51.0 || z>=76.0 || r>=1.0) //cut in vertex + continue; + + if(dis_g1_g2<=12.5) //cut for FCAL splits + continue; + + // replace with Mark's cuts... + if((180.0/TMath::Pi())*theta_det_g1<=2.5 || (180.0/TMath::Pi())*theta_det_g2<=2.5) continue; //cut for beam hole + + if(((180.0/TMath::Pi())*theta_det_g1<=11.5 && (180.0/TMath::Pi())*theta_det_g1>=10.3) || ((180.0/TMath::Pi())*theta_det_g2<=11.5 && (180.0/TMath::Pi())*theta_det_g2>=10.3)) continue; //cut for FCAL-BCAL gap + + //Missing mass off proton + double mmo = locMissingP4_OffProton.M(); + + //Missing Mass squared + double mmsq = locMissingP4_Measured.M2(); + + //Missing Energy + double me = locMissingP4_Measured.E(); + + //Delta_phi + double phi_2g = (180.0/TMath::Pi())*(locPhotonPairP4_Measured.Phi()); + double phi_p = (180/TMath::Pi())*(locProtonP4_Measured.Phi()); + double delta_phi = phi_2g - phi_p; + if(delta_phi>360.) delta_phi += -360.; + if(delta_phi<0.) delta_phi += 360.; + + // Nominal non-PID analysis cuts + if(locBeamP4.E()<6.5 || locBeamP4.E()>11.5 || fabs(delta_phi-180.0)>5.0 || me<-0.50 || me>0.70 || mmo>0.5) + continue; + + /******************************************** EXECUTE ANALYSIS ACTIONS *******************************************/ + + // Loop through the analysis actions, executing them in order for the active particle combo + dAnalyzeCutActions->Perform_ActionWeight(locHistAccidWeightFactor); + + // Note: Using accidental subtraction in AnalyzeCutAction requires version > 1.9.0 of gluex_root_analysis. To use without accidental subtraction and a previous version of the software, comment out the line above and uncomment next line + //dAnalyzeCutActions->Perform_Action(); + + if(!Execute_Actions()) //if the active combo fails a cut, IsComboCutFlag automatically set + continue; + + //if you manually execute any actions, and it fails a cut, be sure to call: + //dComboWrapper->Set_IsComboCut(true); + + // Fill histograms and execute analysis + + } // end of combo loop + + //FILL HISTOGRAMS: Num combos / events surviving actions + Fill_NumCombosSurvivedHists(); + + return kTRUE; +} + +void DSelector_p2gamma::Finalize(void) +{ + //Save anything to output here that you do not want to be in the default DSelector output ROOT file. + + //Otherwise, don't do anything else (especially if you are using PROOF). + //If you are using PROOF, this function is called on each thread, + //so anything you do will not have the combined information from the various threads. + //Besides, it is best-practice to do post-processing (e.g. fitting) separately, in case there is a problem. + + //DO YOUR STUFF HERE + + //CALL THIS LAST + DSelector::Finalize(); //Saves results to the output file +} diff --git a/AnalysisHowTo/PIDstudy/DSelector_p2gamma.h b/AnalysisHowTo/PIDstudy/DSelector_p2gamma.h new file mode 100644 index 00000000..0dbf2dfe --- /dev/null +++ b/AnalysisHowTo/PIDstudy/DSelector_p2gamma.h @@ -0,0 +1,64 @@ +#ifndef DSelector_p2gamma_h +#define DSelector_p2gamma_h + +#include + +#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(dStep0Wrapper->Get_InitialParticle()); + dPhoton1Wrapper = static_cast(dStep0Wrapper->Get_FinalParticle(0)); + dPhoton2Wrapper = static_cast(dStep0Wrapper->Get_FinalParticle(1)); + dProtonWrapper = static_cast(dStep0Wrapper->Get_FinalParticle(2)); +} + +#endif // DSelector_p2gamma_h diff --git a/AnalysisHowTo/PIDstudy/README b/AnalysisHowTo/PIDstudy/README new file mode 100644 index 00000000..e5349e29 --- /dev/null +++ b/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 + diff --git a/AnalysisHowTo/PIDstudy/hist_sum_30274_31057_sim_g4.root b/AnalysisHowTo/PIDstudy/hist_sum_30274_31057_sim_g4.root new file mode 100644 index 00000000..30af2f97 Binary files /dev/null and b/AnalysisHowTo/PIDstudy/hist_sum_30274_31057_sim_g4.root differ diff --git a/AnalysisHowTo/PIDstudy/hist_sum_30730_30788.root b/AnalysisHowTo/PIDstudy/hist_sum_30730_30788.root new file mode 100644 index 00000000..41624e74 Binary files /dev/null and b/AnalysisHowTo/PIDstudy/hist_sum_30730_30788.root differ diff --git a/AnalysisHowTo/PIDstudy/plotPIDdataMC.py b/AnalysisHowTo/PIDstudy/plotPIDdataMC.py new file mode 100644 index 00000000..c848e35a --- /dev/null +++ b/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("/","_"))) diff --git a/AnalysisHowTo/PIDstudy/plotPIDyields.py b/AnalysisHowTo/PIDstudy/plotPIDyields.py new file mode 100644 index 00000000..d05db0cc --- /dev/null +++ b/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) diff --git a/AnalysisHowTo/PIDstudy/runSelector.C b/AnalysisHowTo/PIDstudy/runSelector.C new file mode 100644 index 00000000..ff6a74c3 --- /dev/null +++ b/AnalysisHowTo/PIDstudy/runSelector.C @@ -0,0 +1,56 @@ +// macro to process analysis TTree with TSelector +#include + +#include "TFile.h" +#include "TTree.h" +#include "TString.h" +#include "TSystem.h" + +void runSelector(TString runNumber = "30730", TString myPath = "/cache/halld/RunPeriod-2017-01/analysis/ver38/tree_gg__B4/merged/") +{ + // Load DSelector library + gROOT->ProcessLine(".x $(ROOT_ANALYSIS_HOME)/scripts/Load_DSelector.C"); + int Proof_Nthreads = 8; + + // process signal + TString sampleDir = myPath; + cout<<"running selector on files in: "<GetName(); + if(fileName.Contains(runNumber)) { + cout< skipping"< skipping"<Add(sampleDir+fileName); + ifile++; + } + } + + cout<<"total entries in TChain = "<GetEntries()<<" from "<