// $Id$ // // File: DCustomAction_KsHunt_hists.cc // Created: Wed Jan 21 16:53:41 EST 2015 // Creator: jrsteven (on Linux ifarm1401 2.6.32-431.el6.x86_64 x86_64) // #include "DCustomAction_KsHunt_hists.h" void DCustomAction_KsHunt_hists::Initialize(JEventLoop* locEventLoop) { // get PID algos const DParticleID* locParticleID = NULL; locEventLoop->GetSingle(locParticleID); dParticleID = locParticleID; japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { //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(); if(gDirectory->Get("M2pi") == NULL) //check to see if already created by another thread dM2pi = new TH1D("M2pi", "#pi^{+} #pi^{-} Mass; Mass (GeV/c^{2})", 200, 0.0, 2.0); else //already created by another thread dM2pi = static_cast(gDirectory->Get("M2pi")); if(gDirectory->Get("MM_M2pi") == NULL) //check to see if already created by another thread dMM_M2pi = new TH2D("MM_M2pi", "MM off #pi^{+}#pi^{-} vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM", 200, 0.0, 2.0, 200, 0., 4.); else //already created by another thread dMM_M2pi = static_cast(gDirectory->Get("MM_M2pi")); if(gDirectory->Get("DeltaT2pi") == NULL) //check to see if already created by another thread dDeltaT2pi = new TH2D("DeltaT2pi", "Match SC-TAGGER #Delta t #pi_{1} vs #Delta t #pi_{2}; #Delta t #pi_{2}; #Delta t #pi_{1}", 100, -50.0, 50.0, 100, -50., 50.); else //already created by another thread dDeltaT2pi = static_cast(gDirectory->Get("DeltaT2pi")); if(gDirectory->Get("Egamma2pi") == NULL) //check to see if already created by another thread dEgamma2pi = new TH1D("Egamma2pi", "TAGGER photon energy; E_{gamma}", 320, 2., 10.); else //already created by another thread dEgamma2pi = static_cast(gDirectory->Get("Egamma2pi")); if(gDirectory->Get("MissingDeltaPT_DeltaPhi") == NULL) //check to see if already created by another thread dMissingDeltaPT_DeltaPhi = new TH2D("MissingDeltaPT_DeltaPhi", "Unused Track - MM: #Delta p_{T} vs #Delta t #phi; #Delta #phi; #Delta p_{T}", 180, -180.0, 180.0, 100, -5., 5.); else //already created by another thread dMissingDeltaPT_DeltaPhi = static_cast(gDirectory->Get("MissingDeltaPT_DeltaPhi")); if(gDirectory->Get("MM_M2pi_2SC") == NULL) //check to see if already created by another thread dMM_M2pi_2SC = new TH2D("MM_M2pi_2SC", "MM off #pi^{+}#pi^{-} vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM", 200, 0.0, 2.0, 200, 0., 4.); else //already created by another thread dMM_M2pi_2SC = static_cast(gDirectory->Get("MM_M2pi_2SC")); if(gDirectory->Get("MM_M2pi_Emeson") == NULL) //check to see if already created by another thread dMM_M2pi_Emeson = new TH2D("MM_M2pi_Emeson", "MM off #pi^{+}#pi^{-} vs M_{#pi^{+}#pi^{-}}; M_{#pi^{+}#pi^{-}}; MM", 200, 0.0, 2.0, 200, 0., 4.); else //already created by another thread dMM_M2pi_Emeson = static_cast(gDirectory->Get("MM_M2pi_Emeson")); } japp->RootUnLock(); //RELEASE ROOT LOCK!! } bool DCustomAction_KsHunt_hists::Perform_Action(JEventLoop* locEventLoop, const DParticleCombo* locParticleCombo) { const DDetectorMatches* locDetectorMatches = NULL; locEventLoop->GetSingle(locDetectorMatches); // for unused tracks analysis vector locChargedTracks; locEventLoop->Get(locChargedTracks, "PreSelect"); vector locUnusedChargedTrackHypotheses; for(size_t loc_i = 0; loc_i < locChargedTracks.size(); ++loc_i) locUnusedChargedTrackHypotheses.insert(locUnusedChargedTrackHypotheses.end(), locChargedTracks[loc_i]->dChargedTrackHypotheses.begin(), locChargedTracks[loc_i]->dChargedTrackHypotheses.end()); const DParticleComboStep* locParticleComboStep = locParticleCombo->Get_ParticleComboStep(0); if(locParticleComboStep->Get_InitialParticleID() != Gamma) return true; // get beam photon time const DKinematicData* locBeamPhoton = locParticleCombo->Get_ParticleComboStep(0)->Get_InitialParticle(); double locBeamPhotonTime = locBeamPhoton->time(); double locBeamPhotonEnergy = locBeamPhoton->energy(); if(locBeamPhotonEnergy < 6.8) return true; // cut low energy photons // get final state particles deque locParticles; if(!Get_UseKinFitResultsFlag()) //measured locParticleComboStep->Get_FinalParticles_Measured(locParticles); else locParticleComboStep->Get_FinalParticles(locParticles); // counter for deltaT int nFoundSC = 0; double locDeltaT[2] = {-999.,-999}; // sum 2pi mass DLorentzVector locP4_2pi; for(size_t loc_i = 0; loc_i < locParticles.size(); ++loc_i) { if(locParticleComboStep->Get_FinalParticleID(loc_i) != PiPlus && locParticleComboStep->Get_FinalParticleID(loc_i) != PiMinus) continue; locP4_2pi += locParticles[loc_i]->lorentzMomentum(); // get time based track const DChargedTrack* locChargedTrack = static_cast(locParticleComboStep->Get_FinalParticle_SourceObject(loc_i)); if(locChargedTrack == NULL) continue; // should never happen const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrack->Get_BestFOM(); const DTrackTimeBased* locTrackTimeBased = NULL; locChargedTrackHypothesis->GetSingleT(locTrackTimeBased); // get match to SC DSCHitMatchParams locSCHitMatchParams; bool foundSC = dParticleID->Get_BestSCMatchParams(locTrackTimeBased, locDetectorMatches, locSCHitMatchParams); if(foundSC){ double locSCTime = locSCHitMatchParams.dHitTime - locSCHitMatchParams.dFlightTime; locDeltaT[nFoundSC] = locSCTime - locBeamPhotonTime; nFoundSC++; } // elminate pi+/pi- from unused tracks for(size_t loc_j = 0; loc_j < locUnusedChargedTrackHypotheses.size(); ++loc_j) { if((locUnusedChargedTrackHypotheses[loc_j]->candidateid != locChargedTrackHypothesis->candidateid) || (locUnusedChargedTrackHypotheses[loc_j]->PID() != locParticles[loc_i]->PID())) continue; locUnusedChargedTrackHypotheses.erase(locUnusedChargedTrackHypotheses.begin() + loc_j); } } // calculate missing mass DLorentzVector locMissingP4 = dAnalysisUtilities->Calc_MissingP4(locParticleCombo,Get_UseKinFitResultsFlag()); japp->RootWriteLock(); //ACQUIRE ROOT LOCK!! { // Fill any histograms here dM2pi->Fill(locP4_2pi.M()); dMM_M2pi->Fill(locP4_2pi.M(), locMissingP4.M()); dEgamma2pi->Fill(locBeamPhotonEnergy); // count unused tracks and correlate with 3pi system (already has pre-select) int nUnusedChargedTrackHypotheses = 0; for(size_t loc_j = 0; loc_j < locUnusedChargedTrackHypotheses.size(); ++loc_j) { if(locUnusedChargedTrackHypotheses[loc_j]->PID() != Proton) continue; nUnusedChargedTrackHypotheses++; DLorentzVector locUnusedP4 = locUnusedChargedTrackHypotheses[loc_j]->lorentzMomentum(); double deltaPhi = (locUnusedP4.Vect().Phi() - locMissingP4.Vect().Phi())*180./TMath::Pi(); if(deltaPhi > 180.) deltaPhi -= 360.; else if(deltaPhi < -180.) deltaPhi += 360.; dMissingDeltaPT_DeltaPhi->Fill(deltaPhi, locUnusedP4.Vect().Perp()-locMissingP4.Vect().Perp()); } // 2pi energy cut (Dave M. suggestion) if(locP4_2pi.E() > 5) { dMM_M2pi_Emeson->Fill(locP4_2pi.M(), locMissingP4.M()); } // both tracks have SC hit if(nFoundSC == 2){ dDeltaT2pi->Fill(locDeltaT[0],locDeltaT[1]); dMM_M2pi_2SC->Fill(locP4_2pi.M(), locMissingP4.M()); } } 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!) }