// $Id$ // // File: DCustomAction_ParticleStudies.cc // Created: Thu May 1 21:18:23 EDT 2014 // Creator: pmatt (on Darwin pmattLaptop 10.8.0 i386) // #include "DCustomAction_ParticleStudies.h" void DCustomAction_ParticleStudies::Initialize(JEventLoop* locEventLoop) { //Optional: Create histograms and/or modify member variables. //Create any histograms/trees/etc. within a ROOT lock. //This is so that when running multithreaded, only one thread is writing to the ROOT file at a time. //When creating a reaction-independent action, only modify member variables within a ROOT lock. //Objects created within a plugin (such as reaction-independent actions) can be accessed by many threads simultaneously. japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { // Optional: Useful utility functions. // locEventLoop->GetSingle(dAnalysisUtilities); //Required: Create a folder in the ROOT output file that will contain all of the output ROOT objects (if any) for this action. //If another thread has already created the folder, it just changes to it. CreateAndChangeTo_ActionDirectory(); CreateAndChangeTo_Directory("ChargedTracks", "ChargedTracks"); { //Thrown if(gDirectory->Get("NumThrown") == NULL) //check to see if already created by another thread dHist_NumThrownQTracks = new TH1D("NumThrown", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHist_NumThrownQTracks = static_cast(gDirectory->Get("NumThrown")); //TrueRecon if(gDirectory->Get("NumTrueRecon") == NULL) //check to see if already created by another thread dHist_NumReconTrueQTracks = new TH1D("NumTrueRecon", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHist_NumReconTrueQTracks = static_cast(gDirectory->Get("NumTrueRecon")); //NumDCHitsVsTheta_True if(gDirectory->Get("NumDCHitsVsTheta_True") == NULL) //check to see if already created by another thread dHist_NumDCHitsVsTheta_True = new TH2D("NumDCHitsVsTheta_True", ";#theta#circ;# DC Hits", 280, 0.0, 140.0, 50, 0.5, 50.5); else //already created by another thread dHist_NumDCHitsVsTheta_True = static_cast(gDirectory->Get("NumDCHitsVsTheta_True")); //NumTrueTracking5Sigma if(gDirectory->Get("NumTrueTracking5Sigma") == NULL) //check to see if already created by another thread dHist_NumReconTrueQTracks_Tracking5Sigma = new TH1D("NumTrueTracking5Sigma", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHist_NumReconTrueQTracks_Tracking5Sigma = static_cast(gDirectory->Get("NumTrueTracking5Sigma")); //NumTrueTracking3Sigma if(gDirectory->Get("NumTrueTracking3Sigma") == NULL) //check to see if already created by another thread dHist_NumReconTrueQTracks_Tracking3Sigma = new TH1D("NumTrueTracking3Sigma", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHist_NumReconTrueQTracks_Tracking3Sigma = static_cast(gDirectory->Get("NumTrueTracking3Sigma")); //NumTrueHasDetectorHit if(gDirectory->Get("NumTrueHasDetectorHit") == NULL) //check to see if already created by another thread dHist_NumReconTrueQTracks_HasDetectorHit = new TH1D("NumTrueHasDetectorHit", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHist_NumReconTrueQTracks_HasDetectorHit = static_cast(gDirectory->Get("NumTrueHasDetectorHit")); //NumAllRecon if(gDirectory->Get("NumAllRecon") == NULL) //check to see if already created by another thread dHist_NumReconQTracks = new TH1D("NumAllRecon", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHist_NumReconQTracks = static_cast(gDirectory->Get("NumAllRecon")); //NumAllDCHitsVsTheta if(gDirectory->Get("NumDCHitsVsTheta") == NULL) //check to see if already created by another thread dHist_NumDCHitsVsTheta = new TH2D("NumDCHitsVsTheta", ";#theta#circ;# DC Hits", 280, 0.0, 140.0, 50, 0.5, 50.5); else //already created by another thread dHist_NumDCHitsVsTheta = static_cast(gDirectory->Get("NumDCHitsVsTheta")); //NumAllTracking5Sigma if(gDirectory->Get("NumAllTracking5Sigma") == NULL) //check to see if already created by another thread dHist_NumReconQTracks_Tracking5Sigma = new TH1D("NumAllTracking5Sigma", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHist_NumReconQTracks_Tracking5Sigma = static_cast(gDirectory->Get("NumAllTracking5Sigma")); //NumAllTracking3Sigma if(gDirectory->Get("NumAllTracking3Sigma") == NULL) //check to see if already created by another thread dHist_NumReconQTracks_Tracking3Sigma = new TH1D("NumAllTracking3Sigma", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHist_NumReconQTracks_Tracking3Sigma = static_cast(gDirectory->Get("NumAllTracking3Sigma")); //NumAllHasDetectorHit if(gDirectory->Get("NumAllHasDetectorHit") == NULL) //check to see if already created by another thread dHist_NumReconQTracks_HasDetectorHit = new TH1D("NumAllHasDetectorHit", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHist_NumReconQTracks_HasDetectorHit = static_cast(gDirectory->Get("NumAllHasDetectorHit")); gDirectory->cd(".."); //return to the action directory } for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i) { Particle_t locPID = dFinalStatePIDs[loc_i]; string locParticleName = ParticleType(locPID); string locParticleROOTName = ParticleName_ROOT(locPID); CreateAndChangeTo_Directory(locParticleName, locParticleName); //Thrown if(gDirectory->Get("NumThrown") == NULL) //check to see if already created by another thread dHistMap_NumTracksByPID_Thrown[locPID] = new TH1D("NumThrown", ";# Thrown", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHistMap_NumTracksByPID_Thrown[locPID] = static_cast(gDirectory->Get("NumThrown")); //TrueRecon if(gDirectory->Get("NumTrueRecon") == NULL) //check to see if already created by another thread dHistMap_NumTracksByPID_TrueRecon[locPID] = new TH1D("NumTrueRecon", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHistMap_NumTracksByPID_TrueRecon[locPID] = static_cast(gDirectory->Get("NumTrueRecon")); //TruePID5Sigma if(gDirectory->Get("NumTruePID5Sigma") == NULL) //check to see if already created by another thread dHistMap_NumTracksByPID_TruePID5Sigma[locPID] = new TH1D("NumTruePID5Sigma", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHistMap_NumTracksByPID_TruePID5Sigma[locPID] = static_cast(gDirectory->Get("NumTruePID5Sigma")); //TrueTracking5Sigma //includes previous cuts if(gDirectory->Get("NumTrueTracking5Sigma") == NULL) //check to see if already created by another thread dHistMap_NumTracksByPID_TrueTracking5Sigma[locPID] = new TH1D("NumTrueTracking5Sigma", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHistMap_NumTracksByPID_TrueTracking5Sigma[locPID] = static_cast(gDirectory->Get("NumTrueTracking5Sigma")); //TrueTracking3Sigma //includes previous cuts if(gDirectory->Get("NumTrueTracking3Sigma") == NULL) //check to see if already created by another thread dHistMap_NumTracksByPID_TrueTracking3Sigma[locPID] = new TH1D("NumTrueTracking3Sigma", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHistMap_NumTracksByPID_TrueTracking3Sigma[locPID] = static_cast(gDirectory->Get("NumTrueTracking3Sigma")); //TrueHasDetectorHit if(gDirectory->Get("NumTrueHasDetectorHit") == NULL) //check to see if already created by another thread dHistMap_NumTracksByPID_TrueHasDetectorHit[locPID] = new TH1D("NumTrueHasDetectorHit", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHistMap_NumTracksByPID_TrueHasDetectorHit[locPID] = static_cast(gDirectory->Get("NumTrueHasDetectorHit")); //TruePID3Sigma //includes previous cuts if(gDirectory->Get("NumTruePID3Sigma") == NULL) //check to see if already created by another thread dHistMap_NumTracksByPID_TruePID3Sigma[locPID] = new TH1D("NumTruePID3Sigma", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHistMap_NumTracksByPID_TruePID3Sigma[locPID] = static_cast(gDirectory->Get("NumTruePID3Sigma")); //AllPID5Sigma if(gDirectory->Get("NumAllPID5Sigma") == NULL) //check to see if already created by another thread dHistMap_NumTracksByPID_AllPID5Sigma[locPID] = new TH1D("NumAllPID5Sigma", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHistMap_NumTracksByPID_AllPID5Sigma[locPID] = static_cast(gDirectory->Get("NumAllPID5Sigma")); //AllTracking5Sigma //includes previous cuts if(gDirectory->Get("NumAllTracking5Sigma") == NULL) //check to see if already created by another thread dHistMap_NumTracksByPID_AllTracking5Sigma[locPID] = new TH1D("NumAllTracking5Sigma", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHistMap_NumTracksByPID_AllTracking5Sigma[locPID] = static_cast(gDirectory->Get("NumAllTracking5Sigma")); //AllTracking3Sigma //includes previous cuts if(gDirectory->Get("NumAllTracking3Sigma") == NULL) //check to see if already created by another thread dHistMap_NumTracksByPID_AllTracking3Sigma[locPID] = new TH1D("NumAllTracking3Sigma", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHistMap_NumTracksByPID_AllTracking3Sigma[locPID] = static_cast(gDirectory->Get("NumAllTracking3Sigma")); //AllHasDetectorHit if(gDirectory->Get("NumAllHasDetectorHit") == NULL) //check to see if already created by another thread dHistMap_NumTracksByPID_AllHasDetectorHit[locPID] = new TH1D("NumAllHasDetectorHit", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHistMap_NumTracksByPID_AllHasDetectorHit[locPID] = static_cast(gDirectory->Get("NumAllHasDetectorHit")); //AllPID3Sigma //includes previous cuts if(gDirectory->Get("NumAllPID3Sigma") == NULL) //check to see if already created by another thread dHistMap_NumTracksByPID_AllPID3Sigma[locPID] = new TH1D("NumAllPID3Sigma", ";# Charged Tracks", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHistMap_NumTracksByPID_AllPID3Sigma[locPID] = static_cast(gDirectory->Get("NumAllPID3Sigma")); gDirectory->cd(".."); //return to the action directory } CreateAndChangeTo_Directory("Neutrals", "Neutrals"); { //Thrown if(gDirectory->Get("NumThrown") == NULL) //check to see if already created by another thread dHist_NumThrownNeutralShowers = new TH1D("NumThrown", ";# Neutral Showers", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHist_NumThrownNeutralShowers = static_cast(gDirectory->Get("NumThrown")); //TrueRecon if(gDirectory->Get("NumTrueRecon") == NULL) //check to see if already created by another thread dHist_NumNeutralShowers_True = new TH1D("NumTrueRecon", ";# Neutral Showers", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHist_NumNeutralShowers_True = static_cast(gDirectory->Get("NumTrueRecon")); //NumAllRecon if(gDirectory->Get("NumAllRecon") == NULL) //check to see if already created by another thread dHist_NumNeutralShowers = new TH1D("NumAllRecon", ";# Neutral Showers", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHist_NumNeutralShowers = static_cast(gDirectory->Get("NumAllRecon")); //AllShowerTime if(gDirectory->Get("AllShowerTime") == NULL) //check to see if already created by another thread dHist_NeutralShowerTime_All = new TH1D("AllShowerTime", ";Shower Time (ns)", 1600, -800.0, 800.0); else //already created by another thread dHist_NeutralShowerTime_All = static_cast(gDirectory->Get("AllShowerTime")); //TrueShowerTime if(gDirectory->Get("TrueShowerTime") == NULL) //check to see if already created by another thread dHist_NeutralShowerTime_True = new TH1D("TrueShowerTime", ";Shower Time (ns)", 1600, -800.0, 800.0); else //already created by another thread dHist_NeutralShowerTime_True = static_cast(gDirectory->Get("TrueShowerTime")); //NumTruePID5Sigma if(gDirectory->Get("NumTruePID5Sigma") == NULL) //check to see if already created by another thread dHist_NumNeutralShowers_TruePID5Sigma = new TH1D("NumTruePID5Sigma", ";# Neutral Showers", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHist_NumNeutralShowers_TruePID5Sigma = static_cast(gDirectory->Get("NumTruePID5Sigma")); //NumAllPID5Sigma if(gDirectory->Get("NumAllPID5Sigma") == NULL) //check to see if already created by another thread dHist_NumNeutralShowers_AllPID5Sigma = new TH1D("NumAllPID5Sigma", ";# Neutral Showers", dNumParticles + 1, -0.5, (float)dNumParticles + 0.5); else //already created by another thread dHist_NumNeutralShowers_AllPID5Sigma = static_cast(gDirectory->Get("NumAllPID5Sigma")); //AllShowerEnergy if(gDirectory->Get("AllShowerEnergy") == NULL) //check to see if already created by another thread dHist_NeutralShowerEnergy_All = new TH1D("AllShowerEnergy", ";Shower Energy (GeV)", 1600, 0.0, 8.0); else //already created by another thread dHist_NeutralShowerEnergy_All = static_cast(gDirectory->Get("AllShowerEnergy")); //TrueShowerEnergy if(gDirectory->Get("TrueShowerEnergy") == NULL) //check to see if already created by another thread dHist_NeutralShowerEnergy_True = new TH1D("TrueShowerEnergy", ";Shower Energy (GeV)", 1600, 0.0, 8.0); else //already created by another thread dHist_NeutralShowerEnergy_True = static_cast(gDirectory->Get("TrueShowerEnergy")); //AllShowerEnergy - BCAL if(gDirectory->Get("AllShowerEnergy_BCAL") == NULL) //check to see if already created by another thread dHist_NeutralShowerEnergy_All_BCAL = new TH1D("AllShowerEnergy_BCAL", ";BCAL Shower Energy (GeV)", 1600, 0.0, 8.0); else //already created by another thread dHist_NeutralShowerEnergy_All_BCAL = static_cast(gDirectory->Get("AllShowerEnergy_BCAL")); //TrueShowerEnergy - BCAL if(gDirectory->Get("TrueShowerEnergy_BCAL") == NULL) //check to see if already created by another thread dHist_NeutralShowerEnergy_True_BCAL = new TH1D("TrueShowerEnergy_BCAL", ";BCAL Shower Energy (GeV)", 1600, 0.0, 8.0); else //already created by another thread dHist_NeutralShowerEnergy_True_BCAL = static_cast(gDirectory->Get("TrueShowerEnergy_BCAL")); //AllShowerEnergy - FCAL if(gDirectory->Get("AllShowerEnergy_FCAL") == NULL) //check to see if already created by another thread dHist_NeutralShowerEnergy_All_FCAL = new TH1D("AllShowerEnergy_FCAL", ";FCAL Shower Energy (GeV)", 1600, 0.0, 8.0); else //already created by another thread dHist_NeutralShowerEnergy_All_FCAL = static_cast(gDirectory->Get("AllShowerEnergy_FCAL")); //TrueShowerEnergy - FCAL if(gDirectory->Get("TrueShowerEnergy_FCAL") == NULL) //check to see if already created by another thread dHist_NeutralShowerEnergy_True_FCAL = new TH1D("TrueShowerEnergy_FCAL", ";FCAL Shower Energy (GeV)", 1600, 0.0, 8.0); else //already created by another thread dHist_NeutralShowerEnergy_True_FCAL = static_cast(gDirectory->Get("TrueShowerEnergy_FCAL")); //DistanceToNearestTrack - BCAL if(gDirectory->Get("AllDistanceToNearestTrack_BCAL") == NULL) //check to see if already created by another thread dHist_DistanceToNearestTrack_All_BCAL = new TH1D("AllDistanceToNearestTrack_BCAL", ";BCAL Distance to Nearest Track (cm)", 1600, 0.0, 400.0); else //already created by another thread dHist_DistanceToNearestTrack_All_BCAL = static_cast(gDirectory->Get("AllDistanceToNearestTrack_BCAL")); //DistanceToNearestTrackVsTheta - BCAL if(gDirectory->Get("AllDistanceToNearestTrackVsTheta_BCAL") == NULL) //check to see if already created by another thread dHist_DistanceToNearestTrackVsTheta_All_BCAL = new TH2D("AllDistanceToNearestTrackVsTheta_BCAL", ";#theta#circ;BCAL Distance to Nearest Track (cm)", 280, 0.0, 140.0, 1600, 0.0, 400.0); else //already created by another thread dHist_DistanceToNearestTrackVsTheta_All_BCAL = static_cast(gDirectory->Get("AllDistanceToNearestTrackVsTheta_BCAL")); //DistanceToNearestTrack - FCAL if(gDirectory->Get("AllDistanceToNearestTrack_FCAL") == NULL) //check to see if already created by another thread dHist_DistanceToNearestTrack_All_FCAL = new TH1D("AllDistanceToNearestTrack_FCAL", ";FCAL Distance to Nearest Track (cm)", 1600, 0.0, 400.0); else //already created by another thread dHist_DistanceToNearestTrack_All_FCAL = static_cast(gDirectory->Get("AllDistanceToNearestTrack_FCAL")); //TrueDistanceToNearestTrack - BCAL if(gDirectory->Get("TrueDistanceToNearestTrack_BCAL") == NULL) //check to see if already created by another thread dHist_DistanceToNearestTrack_True_BCAL = new TH1D("TrueDistanceToNearestTrack_BCAL", ";BCAL Distance to Nearest Track (cm)", 1600, 0.0, 400.0); else //already created by another thread dHist_DistanceToNearestTrack_True_BCAL = static_cast(gDirectory->Get("TrueDistanceToNearestTrack_BCAL")); //TrueDistanceToNearestTrackVsTheta - BCAL if(gDirectory->Get("TrueDistanceToNearestTrackVsTheta_BCAL") == NULL) //check to see if already created by another thread dHist_DistanceToNearestTrackVsTheta_True_BCAL = new TH2D("TrueDistanceToNearestTrackVsTheta_BCAL", ";#theta#circ;BCAL Distance to Nearest Track (cm)", 280, 0.0, 140.0, 1600, 0.0, 400.0); else //already created by another thread dHist_DistanceToNearestTrackVsTheta_True_BCAL = static_cast(gDirectory->Get("TrueDistanceToNearestTrackVsTheta_BCAL")); //TrueDistanceToNearestTrack - FCAL if(gDirectory->Get("TrueDistanceToNearestTrack_FCAL") == NULL) //check to see if already created by another thread dHist_DistanceToNearestTrack_True_FCAL = new TH1D("TrueDistanceToNearestTrack_FCAL", ";FCAL Distance to Nearest Track (cm)", 1600, 0.0, 400.0); else //already created by another thread dHist_DistanceToNearestTrack_True_FCAL = static_cast(gDirectory->Get("TrueDistanceToNearestTrack_FCAL")); gDirectory->cd(".."); //return to the action directory } //NumQTracks2D if(gDirectory->Get("NumQTracks2D") == NULL) //check to see if already created by another thread dHist_NumQTracks2D = new TH2D("NumQTracks2D", ";#geq # q^{+};#geq # q^{-}", 11, -0.5, 10.5, 11, -0.5, 10.5); else //already created by another thread dHist_NumQTracks2D = static_cast(gDirectory->Get("NumQTracks2D")); //NumTrueQTracks2D if(gDirectory->Get("NumTrueQTracks2D") == NULL) //check to see if already created by another thread dHist_NumTrueQTracks2D = new TH2D("NumTrueQTracks2D", ";#geq # q^{+};#geq # q^{-}", 11, -0.5, 10.5, 11, -0.5, 10.5); else //already created by another thread dHist_NumTrueQTracks2D = static_cast(gDirectory->Get("NumTrueQTracks2D")); } japp->RootUnLock(); //RELEASE ROOT LOCK!! } bool DCustomAction_ParticleStudies::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) { //Expect locParticleCombo to be NULL since this is a reaction-independent action. //Optional: Quit the action if it has already been executed this event (else may result in double-counting when filling histograms) if(Get_NumPreviousParticleCombos() != 0) return true; const DDetectorMatches* locDetectorMatches = NULL; locEventLoop->GetSingle(locDetectorMatches); vector locNeutralParticles; locEventLoop->Get(locNeutralParticles); vector locNeutralShowers; locEventLoop->Get(locNeutralShowers); vector locChargedTracks; locEventLoop->Get(locChargedTracks); vector locMCThrowns; locEventLoop->Get(locMCThrowns, "FinalState"); const DMCThrownMatching* locMCThrownMatching = NULL; locEventLoop->GetSingle(locMCThrownMatching); size_t locNumThrownCharged = 0, locNumThrownNeutral = 0; map locNumThrownByPID; for(size_t loc_i = 0; loc_i < locMCThrowns.size(); ++loc_i) { Particle_t locPID = locMCThrowns[loc_i]->PID(); if(ParticleCharge(locPID) == 0) ++locNumThrownNeutral; if(locNumThrownByPID.find(locPID) == locNumThrownByPID.end()) locNumThrownByPID[locPID] = 1; else ++(locNumThrownByPID[locPID]); } locNumThrownCharged = locMCThrowns.size() - locNumThrownNeutral; //charged true recon map locNumTracksByPID_TrueRecon; map locNumTracksByPID_TruePID5Sigma; map locNumTracksByPID_TrueTracking5Sigma; map locNumTracksByPID_TrueTracking3Sigma; map locNumTracksByPID_TrueHasDetectorHit; map locNumTracksByPID_TruePID3Sigma; unsigned int locNumReconTrueQTracks = 0, locNumReconTrueQTracks_Tracking5Sigma = 0, locNumReconTrueQTracks_Tracking3Sigma = 0, locNumReconTrueQTracks_HasDetectorHit = 0; unsigned int locNumTrueTracksDefinitelyQPos = 0, locNumTrueTracksDefinitelyQNeg = 0; for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i) { if(locMCThrownMatching == NULL) break; double locMatchFOM = 0.0; const DMCThrown* locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locChargedTracks[loc_i], locMatchFOM); if((locMCThrown == NULL) || (locMatchFOM < 5.73303E-7)) continue; Particle_t locPID = (Particle_t)locMCThrown->type; const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_Hypothesis(locPID); if(locChargedTrackHypothesis == NULL) locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM(); //finally, have true recon ++locNumReconTrueQTracks; bool locHasDetectorHit = true; do { if(locChargedTrackHypothesis->dSCHitMatchParams.dTrackTimeBased != NULL) break; if(locChargedTrackHypothesis->dTOFHitMatchParams.dTrackTimeBased != NULL) break; if(locChargedTrackHypothesis->dBCALShowerMatchParams.dTrackTimeBased != NULL) break; if(locChargedTrackHypothesis->dFCALShowerMatchParams.dTrackTimeBased != NULL) break; locHasDetectorHit = false; } while(false); double locTrackingFOM = TMath::Prob(locChargedTrackHypothesis->dChiSq_Track, locChargedTrackHypothesis->dNDF_Track); if(locTrackingFOM >= 5.73303E-7) ++locNumReconTrueQTracks_Tracking5Sigma; if(locTrackingFOM >= 0.0027) { ++locNumReconTrueQTracks_Tracking3Sigma; if(locHasDetectorHit) { ++locNumReconTrueQTracks_HasDetectorHit; if(ParticleCharge(locPID) > 0.0) ++locNumTrueTracksDefinitelyQPos; else ++locNumTrueTracksDefinitelyQNeg; } } if(locNumTracksByPID_TrueRecon.find(locPID) == locNumTracksByPID_TrueRecon.end()) locNumTracksByPID_TrueRecon[locPID] = 1; else ++(locNumTracksByPID_TrueRecon[locPID]); //pid 5 sigma if(!Cut_PID(locChargedTrackHypothesis, 5.73303E-7)) continue; if(locNumTracksByPID_TruePID5Sigma.find(locPID) == locNumTracksByPID_TruePID5Sigma.end()) locNumTracksByPID_TruePID5Sigma[locPID] = 1; else ++(locNumTracksByPID_TruePID5Sigma[locPID]); //tracking 5 sigma if(locTrackingFOM < 5.73303E-7) continue; if(locNumTracksByPID_TrueTracking5Sigma.find(locPID) == locNumTracksByPID_TrueTracking5Sigma.end()) locNumTracksByPID_TrueTracking5Sigma[locPID] = 1; else ++(locNumTracksByPID_TrueTracking5Sigma[locPID]); //tracking 3 sigma if(locTrackingFOM < 0.0027) continue; if(locNumTracksByPID_TrueTracking3Sigma.find(locPID) == locNumTracksByPID_TrueTracking3Sigma.end()) locNumTracksByPID_TrueTracking3Sigma[locPID] = 1; else ++(locNumTracksByPID_TrueTracking3Sigma[locPID]); //has detector hit if(!locHasDetectorHit) continue; if(locNumTracksByPID_TrueHasDetectorHit.find(locPID) == locNumTracksByPID_TrueHasDetectorHit.end()) locNumTracksByPID_TrueHasDetectorHit[locPID] = 1; else ++(locNumTracksByPID_TrueHasDetectorHit[locPID]); //pid 3 sigma if(!Cut_PID(locChargedTrackHypothesis, 0.0027)) continue; if(locNumTracksByPID_TruePID3Sigma.find(locPID) == locNumTracksByPID_TruePID3Sigma.end()) locNumTracksByPID_TruePID3Sigma[locPID] = 1; else ++(locNumTracksByPID_TruePID3Sigma[locPID]); } //charged all recon map locNumTracksByPID_AllPID5Sigma; map locNumTracksByPID_AllTracking5Sigma; map locNumTracksByPID_AllTracking3Sigma; map locNumTracksByPID_AllHasDetectorHit; map locNumTracksByPID_AllPID3Sigma; unsigned int locNumReconQTracks_Tracking5Sigma = 0, locNumReconQTracks_Tracking3Sigma = 0, locNumReconQTracks_HasDetectorHit = 0; unsigned int locNumAllTracksDefinitelyQPos = 0, locNumAllTracksDefinitelyQNeg = 0, locNumAllTracksUnknownQ = 0; for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i) { bool locCouldBeQPos = false, locCouldBeQNeg = false; for(size_t loc_j = 0; loc_j < locChargedTracks[loc_i]->dChargedTrackHypotheses.size(); ++loc_j) { const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->dChargedTrackHypotheses[loc_j]; if(locChargedTrackHypothesis->charge() > 0.0) locCouldBeQPos = true; else locCouldBeQNeg = true; } bool locTracking5SigmaFoundFlag = false, locTracking3SigmaFoundFlag = false; for(size_t loc_j = 0; loc_j < locChargedTracks[loc_i]->dChargedTrackHypotheses.size(); ++loc_j) { const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->dChargedTrackHypotheses[loc_j]; double locTrackingFOM = TMath::Prob(locChargedTrackHypothesis->dChiSq_Track, locChargedTrackHypothesis->dNDF_Track); if((locTrackingFOM >= 5.73303E-7) && (!locTracking5SigmaFoundFlag)) { locTracking5SigmaFoundFlag = true; ++locNumReconQTracks_Tracking5Sigma; } if((locTrackingFOM >= 0.0027) && (!locTracking3SigmaFoundFlag)) { locTracking3SigmaFoundFlag = true; ++locNumReconQTracks_Tracking3Sigma; } bool locHasDetectorHit = true; do { if(locChargedTrackHypothesis->dSCHitMatchParams.dTrackTimeBased != NULL) break; if(locChargedTrackHypothesis->dTOFHitMatchParams.dTrackTimeBased != NULL) break; if(locChargedTrackHypothesis->dBCALShowerMatchParams.dTrackTimeBased != NULL) break; if(locChargedTrackHypothesis->dFCALShowerMatchParams.dTrackTimeBased != NULL) break; locHasDetectorHit = false; } while(false); if(locHasDetectorHit && (locTrackingFOM >= 0.0027)) { ++locNumReconQTracks_HasDetectorHit; if(locCouldBeQPos && locCouldBeQNeg) ++locNumAllTracksUnknownQ; else if(locCouldBeQPos) ++locNumAllTracksDefinitelyQPos; else ++locNumAllTracksDefinitelyQNeg; break; //no need to keep searching } } for(size_t loc_j = 0; loc_j < dFinalStatePIDs.size(); ++loc_j) { Particle_t locPID = dFinalStatePIDs[loc_j]; const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_Hypothesis(locPID); if(locChargedTrackHypothesis == NULL) continue; //pid 5 sigma if(!Cut_PID(locChargedTrackHypothesis, 5.73303E-7)) continue; if(locNumTracksByPID_AllPID5Sigma.find(locPID) == locNumTracksByPID_AllPID5Sigma.end()) locNumTracksByPID_AllPID5Sigma[locPID] = 1; else ++(locNumTracksByPID_AllPID5Sigma[locPID]); //tracking 5 sigma double locTrackingFOM = TMath::Prob(locChargedTrackHypothesis->dChiSq_Track, locChargedTrackHypothesis->dNDF_Track); if(locTrackingFOM < 5.73303E-7) continue; if(locNumTracksByPID_AllTracking5Sigma.find(locPID) == locNumTracksByPID_AllTracking5Sigma.end()) locNumTracksByPID_AllTracking5Sigma[locPID] = 1; else ++(locNumTracksByPID_AllTracking5Sigma[locPID]); //tracking 3 sigma if(locTrackingFOM < 0.0027) continue; if(locNumTracksByPID_AllTracking3Sigma.find(locPID) == locNumTracksByPID_AllTracking3Sigma.end()) locNumTracksByPID_AllTracking3Sigma[locPID] = 1; else ++(locNumTracksByPID_AllTracking3Sigma[locPID]); //has detector hit bool locHasDetectorHit = true; do { if(locChargedTrackHypothesis->dSCHitMatchParams.dTrackTimeBased != NULL) break; if(locChargedTrackHypothesis->dTOFHitMatchParams.dTrackTimeBased != NULL) break; if(locChargedTrackHypothesis->dBCALShowerMatchParams.dTrackTimeBased != NULL) break; if(locChargedTrackHypothesis->dFCALShowerMatchParams.dTrackTimeBased != NULL) break; locHasDetectorHit = false; } while(false); if(!locHasDetectorHit) continue; if(locNumTracksByPID_AllHasDetectorHit.find(locPID) == locNumTracksByPID_AllHasDetectorHit.end()) locNumTracksByPID_AllHasDetectorHit[locPID] = 1; else ++(locNumTracksByPID_AllHasDetectorHit[locPID]); //pid 3 sigma if(!Cut_PID(locChargedTrackHypothesis, 0.0027)) continue; if(locNumTracksByPID_AllPID3Sigma.find(locPID) == locNumTracksByPID_AllPID3Sigma.end()) locNumTracksByPID_AllPID3Sigma[locPID] = 1; else ++(locNumTracksByPID_AllPID3Sigma[locPID]); } } //neutral true recon unsigned int locNumReconTrueShowers = 0; for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i) { if(locMCThrownMatching == NULL) break; double locMatchFOM = 0.0; const DMCThrown* locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locNeutralParticles[loc_i], locMatchFOM); if((locMCThrown == NULL) || (locMatchFOM < 5.73303E-7)) continue; Particle_t locPID = (Particle_t)locMCThrown->type; const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_Hypothesis(locPID); if(locNeutralParticleHypothesis == NULL) locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_BestFOM(); //finally, have true recon ++locNumReconTrueShowers; } japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { // Fill any histograms here dHist_NumThrownQTracks->Fill(locNumThrownCharged); dHist_NumReconTrueQTracks->Fill(locNumReconTrueQTracks); dHist_NumReconTrueQTracks_Tracking5Sigma->Fill(locNumReconTrueQTracks_Tracking5Sigma); dHist_NumReconTrueQTracks_Tracking3Sigma->Fill(locNumReconTrueQTracks_Tracking3Sigma); dHist_NumReconTrueQTracks_HasDetectorHit->Fill(locNumReconTrueQTracks_HasDetectorHit); dHist_NumReconQTracks->Fill(locChargedTracks.size()); dHist_NumReconQTracks_Tracking5Sigma->Fill(locNumReconQTracks_Tracking5Sigma); dHist_NumReconQTracks_Tracking3Sigma->Fill(locNumReconQTracks_Tracking3Sigma); dHist_NumReconQTracks_HasDetectorHit->Fill(locNumReconQTracks_HasDetectorHit); for(size_t loc_i = 0; loc_i < dFinalStatePIDs.size(); ++loc_i) { Particle_t locPID = dFinalStatePIDs[loc_i]; if(locNumThrownByPID.find(locPID) == locNumThrownByPID.end()) dHistMap_NumTracksByPID_Thrown[locPID]->Fill(0); else dHistMap_NumTracksByPID_Thrown[locPID]->Fill(locNumThrownByPID[locPID]); //true recons if(locNumTracksByPID_TrueRecon.find(locPID) == locNumTracksByPID_TrueRecon.end()) dHistMap_NumTracksByPID_TrueRecon[locPID]->Fill(0); else dHistMap_NumTracksByPID_TrueRecon[locPID]->Fill(locNumTracksByPID_TrueRecon[locPID]); if(locNumTracksByPID_TruePID5Sigma.find(locPID) == locNumTracksByPID_TruePID5Sigma.end()) dHistMap_NumTracksByPID_TruePID5Sigma[locPID]->Fill(0); else dHistMap_NumTracksByPID_TruePID5Sigma[locPID]->Fill(locNumTracksByPID_TruePID5Sigma[locPID]); if(locNumTracksByPID_TrueTracking5Sigma.find(locPID) == locNumTracksByPID_TrueTracking5Sigma.end()) dHistMap_NumTracksByPID_TrueTracking5Sigma[locPID]->Fill(0); else dHistMap_NumTracksByPID_TrueTracking5Sigma[locPID]->Fill(locNumTracksByPID_TrueTracking5Sigma[locPID]); if(locNumTracksByPID_TrueTracking3Sigma.find(locPID) == locNumTracksByPID_TrueTracking3Sigma.end()) dHistMap_NumTracksByPID_TrueTracking3Sigma[locPID]->Fill(0); else dHistMap_NumTracksByPID_TrueTracking3Sigma[locPID]->Fill(locNumTracksByPID_TrueTracking3Sigma[locPID]); if(locNumTracksByPID_TrueHasDetectorHit.find(locPID) == locNumTracksByPID_TrueHasDetectorHit.end()) dHistMap_NumTracksByPID_TrueHasDetectorHit[locPID]->Fill(0); else dHistMap_NumTracksByPID_TrueHasDetectorHit[locPID]->Fill(locNumTracksByPID_TrueHasDetectorHit[locPID]); if(locNumTracksByPID_TruePID3Sigma.find(locPID) == locNumTracksByPID_TruePID3Sigma.end()) dHistMap_NumTracksByPID_TruePID3Sigma[locPID]->Fill(0); else dHistMap_NumTracksByPID_TruePID3Sigma[locPID]->Fill(locNumTracksByPID_TruePID3Sigma[locPID]); //all if(locNumTracksByPID_AllPID5Sigma.find(locPID) == locNumTracksByPID_AllPID5Sigma.end()) dHistMap_NumTracksByPID_AllPID5Sigma[locPID]->Fill(0); else dHistMap_NumTracksByPID_AllPID5Sigma[locPID]->Fill(locNumTracksByPID_AllPID5Sigma[locPID]); if(locNumTracksByPID_AllTracking5Sigma.find(locPID) == locNumTracksByPID_AllTracking5Sigma.end()) dHistMap_NumTracksByPID_AllTracking5Sigma[locPID]->Fill(0); else dHistMap_NumTracksByPID_AllTracking5Sigma[locPID]->Fill(locNumTracksByPID_AllTracking5Sigma[locPID]); if(locNumTracksByPID_AllTracking3Sigma.find(locPID) == locNumTracksByPID_AllTracking3Sigma.end()) dHistMap_NumTracksByPID_AllTracking3Sigma[locPID]->Fill(0); else dHistMap_NumTracksByPID_AllTracking3Sigma[locPID]->Fill(locNumTracksByPID_AllTracking3Sigma[locPID]); if(locNumTracksByPID_AllHasDetectorHit.find(locPID) == locNumTracksByPID_AllHasDetectorHit.end()) dHistMap_NumTracksByPID_AllHasDetectorHit[locPID]->Fill(0); else dHistMap_NumTracksByPID_AllHasDetectorHit[locPID]->Fill(locNumTracksByPID_AllHasDetectorHit[locPID]); if(locNumTracksByPID_AllPID3Sigma.find(locPID) == locNumTracksByPID_AllPID3Sigma.end()) dHistMap_NumTracksByPID_AllPID3Sigma[locPID]->Fill(0); else dHistMap_NumTracksByPID_AllPID3Sigma[locPID]->Fill(locNumTracksByPID_AllPID3Sigma[locPID]); } dHist_NumThrownNeutralShowers->Fill(locNumThrownNeutral); dHist_NumNeutralShowers->Fill(locNeutralShowers.size()); dHist_NumNeutralShowers_True->Fill(locNumReconTrueShowers); //track data for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i) { //all for(size_t loc_j = 0; loc_j < locChargedTracks[loc_i]->dChargedTrackHypotheses.size(); ++loc_j) { const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->dChargedTrackHypotheses[loc_j]; dHist_NumDCHitsVsTheta->Fill(locChargedTrackHypothesis->momentum().Theta()*180.0/TMath::Pi(), locChargedTrackHypothesis->dNDF_Track + 5); } //true if(locMCThrownMatching == NULL) continue; double locMatchFOM = 0.0; const DMCThrown* locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locChargedTracks[loc_i], locMatchFOM); if((locMCThrown == NULL) || (locMatchFOM < 5.73303E-7)) continue; Particle_t locPID = (Particle_t)locMCThrown->type; const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_Hypothesis(locPID); if(locChargedTrackHypothesis == NULL) locChargedTrackHypothesis = locChargedTracks[loc_i]->Get_BestFOM(); //finally, have true recon dHist_NumDCHitsVsTheta_True->Fill(locChargedTrackHypothesis->momentum().Theta()*180.0/TMath::Pi(), locChargedTrackHypothesis->dNDF_Track + 5); } //all shower data size_t locNumAllShowersPID5Sigma = 0; for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i) { const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_BestFOM(); dHist_NeutralShowerTime_All->Fill(locNeutralParticleHypothesis->t1()); //pid 5 sigma cut if((locNeutralParticleHypothesis->dNDF > 0) && (locNeutralParticleHypothesis->dFOM < 5.73303E-7)) continue; ++locNumAllShowersPID5Sigma; dHist_NeutralShowerEnergy_All->Fill(locNeutralParticles[loc_i]->dNeutralShower->dEnergy); if(locNeutralParticleHypothesis->t1_detector() == SYS_BCAL) { dHist_NeutralShowerEnergy_All_BCAL->Fill(locNeutralParticles[loc_i]->dNeutralShower->dEnergy); const DBCALShower* locBCALShower = NULL; locNeutralParticles[loc_i]->dNeutralShower->GetSingle(locBCALShower); double locDistance = 0.0; if(locDetectorMatches->Get_DistanceToNearestTrack(locBCALShower, locDistance)) { dHist_DistanceToNearestTrack_All_BCAL->Fill(locDistance); dHist_DistanceToNearestTrackVsTheta_All_BCAL->Fill(locNeutralParticleHypothesis->momentum().Theta()*180.0/TMath::Pi(), locDistance); } } else if(locNeutralParticleHypothesis->t1_detector() == SYS_FCAL) { dHist_NeutralShowerEnergy_All_FCAL->Fill(locNeutralParticles[loc_i]->dNeutralShower->dEnergy); const DFCALShower* locFCALShower = NULL; locNeutralParticles[loc_i]->dNeutralShower->GetSingle(locFCALShower); double locDistance = 0.0; if(locDetectorMatches->Get_DistanceToNearestTrack(locFCALShower, locDistance)) dHist_DistanceToNearestTrack_All_FCAL->Fill(locDistance); } } dHist_NumNeutralShowers_AllPID5Sigma->Fill(locNumAllShowersPID5Sigma); //true shower data size_t locNumTrueShowersPID5Sigma = 0; for(size_t loc_i = 0; loc_i < locNeutralParticles.size(); ++loc_i) { if(locMCThrownMatching == NULL) break; double locMatchFOM = 0.0; const DMCThrown* locMCThrown = locMCThrownMatching->Get_MatchingMCThrown(locNeutralParticles[loc_i], locMatchFOM); if((locMCThrown == NULL) || (locMatchFOM < 5.73303E-7)) continue; Particle_t locPID = (Particle_t)locMCThrown->type; const DNeutralParticleHypothesis* locNeutralParticleHypothesis = locNeutralParticles[loc_i]->Get_Hypothesis(locPID); if(locNeutralParticleHypothesis == NULL) continue; //finally, have a true shower dHist_NeutralShowerTime_True->Fill(locNeutralParticleHypothesis->t1()); //pid 5 sigma cut if((locNeutralParticleHypothesis->dNDF > 0) && (locNeutralParticleHypothesis->dFOM < 5.73303E-7)) continue; ++locNumTrueShowersPID5Sigma; dHist_NeutralShowerEnergy_True->Fill(locNeutralParticles[loc_i]->dNeutralShower->dEnergy); if(locNeutralParticleHypothesis->t1_detector() == SYS_BCAL) { dHist_NeutralShowerEnergy_True_BCAL->Fill(locNeutralParticles[loc_i]->dNeutralShower->dEnergy); const DBCALShower* locBCALShower = NULL; locNeutralParticles[loc_i]->dNeutralShower->GetSingle(locBCALShower); double locDistance = 0.0; if(locDetectorMatches->Get_DistanceToNearestTrack(locBCALShower, locDistance)) { dHist_DistanceToNearestTrack_True_BCAL->Fill(locDistance); dHist_DistanceToNearestTrackVsTheta_True_BCAL->Fill(locNeutralParticleHypothesis->momentum().Theta()*180.0/TMath::Pi(), locDistance); } } else if(locNeutralParticleHypothesis->t1_detector() == SYS_FCAL) { dHist_NeutralShowerEnergy_True_FCAL->Fill(locNeutralParticles[loc_i]->dNeutralShower->dEnergy); const DFCALShower* locFCALShower = NULL; locNeutralParticles[loc_i]->dNeutralShower->GetSingle(locFCALShower); double locDistance = 0.0; if(locDetectorMatches->Get_DistanceToNearestTrack(locFCALShower, locDistance)) dHist_DistanceToNearestTrack_True_FCAL->Fill(locDistance); } } dHist_NumNeutralShowers_TruePID5Sigma->Fill(locNumTrueShowersPID5Sigma); //recon tracks 2d for(int loc_i = 1; loc_i <= dHist_NumQTracks2D->GetNbinsX(); ++loc_i) { if((loc_i - 1) > (locNumAllTracksDefinitelyQPos + locNumAllTracksUnknownQ)) break; for(int loc_j = 1; loc_j <= dHist_NumQTracks2D->GetNbinsY(); ++loc_j) { if((loc_j - 1) > (locNumAllTracksDefinitelyQNeg + locNumAllTracksUnknownQ)) break; if((loc_i + loc_j - 2) > (locNumAllTracksDefinitelyQPos + locNumAllTracksDefinitelyQNeg + locNumAllTracksUnknownQ)) break; dHist_NumQTracks2D->Fill(loc_i - 1, loc_j - 1); } } //true tracks 2d for(int loc_i = 1; loc_i <= dHist_NumTrueQTracks2D->GetNbinsX(); ++loc_i) { if((loc_i - 1) > locNumTrueTracksDefinitelyQPos) break; for(int loc_j = 1; loc_j <= dHist_NumTrueQTracks2D->GetNbinsY(); ++loc_j) { if((loc_j - 1) > locNumTrueTracksDefinitelyQNeg) break; dHist_NumTrueQTracks2D->Fill(loc_i - 1, loc_j - 1); } } } japp->RootUnLock(); //RELEASE ROOT LOCK!! return true; //return false if you want to use this action to apply a cut (and it fails the cut!) } bool DCustomAction_ParticleStudies::Cut_PID(const DNeutralParticleHypothesis* locNeutralParticleHypothesis, double locMinimumConfidenceLevel) { return ((locNeutralParticleHypothesis->dNDF == 0) ? true : (TMath::Prob(locNeutralParticleHypothesis->dChiSq, locNeutralParticleHypothesis->dNDF) >= locMinimumConfidenceLevel)); } bool DCustomAction_ParticleStudies::Cut_PID(const DChargedTrackHypothesis* locChargedTrackHypothesis, double locMinimumConfidenceLevel) { if(locChargedTrackHypothesis->t1_detector() == SYS_BCAL) return true; //don't cut if bcal return ((locChargedTrackHypothesis->dNDF == 0) ? true : (TMath::Prob(locChargedTrackHypothesis->dChiSq, locChargedTrackHypothesis->dNDF) >= locMinimumConfidenceLevel)); }