From be537991fbd4408e38fd5a0ecde986c8b1956ed0 Mon Sep 17 00:00:00 2001 From: gxproj7 gxproj7 Date: Mon, 22 Apr 2019 10:00:49 -0400 Subject: [PATCH 01/19] change to dirc_tree plugin for to properly compute invariant mass --- .../Analysis/dirc_tree/DEventProcessor_dirc_tree.cc | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc b/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc index 952dcf1da..3134aa864 100644 --- a/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc +++ b/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc @@ -111,13 +111,21 @@ jerror_t DEventProcessor_dirc_tree::evnt(jana::JEventLoop* loop, uint64_t locEve for(size_t icombo = 0; icombo < locPassedParticleCombos.size(); ++icombo){ double chisq = locPassedParticleCombos[icombo]->Get_KinFitResults()->Get_ChiSq(); DLorentzVector locMissingP4 = fAnalysisUtilities->Calc_MissingP4(locReaction, locPassedParticleCombos[icombo], false); - DLorentzVector locInvP4; auto locParticleComboStep = locPassedParticleCombos[icombo]->Get_ParticleComboStep(0); + + // calculate pi+pi- and K+K- invariant mass by taking the final state and subtracting the proton + DLorentzVector locInvP4 = fAnalysisUtilities->Calc_FinalStateP4(locReaction, locPassedParticleCombos[icombo], 0, false); + for(size_t parti=0; partiGet_NumFinalParticles(); parti++){ + auto locParticle = locParticleComboStep->Get_FinalParticle(parti); + if(locParticle->PID() == Proton) locInvP4 -= locParticle->lorentzMomentum(); + } int numIt=0; for(size_t parti=0; partiGet_NumFinalParticles(); parti++){ auto locParticle = locParticleComboStep->Get_FinalParticle(parti); // Get_FinalParticle_Measured(parti); + if(locParticle->charge() == 0) continue; +/* // we expect only rho or phi events: g,p->pi+,pi-,p g,p->K+,K-,p if(locParticle->PID() == PiPlus || locParticle->PID() == PiMinus || locParticle->PID() == KPlus || locParticle->PID() == KMinus){ locInvP4 += locParticle->lorentzMomentum(); @@ -129,6 +137,7 @@ jerror_t DEventProcessor_dirc_tree::evnt(jana::JEventLoop* loop, uint64_t locEve previousInv.push_back(tm); } } +*/ // get track auto locChargedTrack = static_cast(locParticleComboStep->Get_FinalParticle_SourceObject(parti)); From 55416b0a07ab6fe83158a6b41bd083eebdd0c3e9 Mon Sep 17 00:00:00 2001 From: gxproj7 gxproj7 Date: Mon, 22 Apr 2019 15:41:13 -0400 Subject: [PATCH 02/19] fix dirc_tree InvMass calculation and add lines for loading library by Ahmed --- .../Analysis/dirc_tree/DEventProcessor_dirc_tree.cc | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc b/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc index 3134aa864..76ebd0b02 100644 --- a/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc +++ b/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc @@ -63,6 +63,14 @@ jerror_t DEventProcessor_dirc_tree::init(void) //------------------ jerror_t DEventProcessor_dirc_tree::brun(jana::JEventLoop* locEventLoop, int locRunNumber) { + ////////////////////////////////////////////////////////////////////////////// + // dapp and geom are not used but without it dirc_hits.so will no be loaded + // root [0] gSystem->Load("dirc_hits.so"); + // cling::DynamicLibraryManager::loadLibrary(): dirc_hits.so: undefined symbol: _ZTV25DEventSourceRESTGenerator + DApplication* dapp=dynamic_cast(locEventLoop->GetJApplication()); + DGeometry *geom = dapp->GetDGeometry(locRunNumber); + ////////////////////////////////////////////////////////////////////////////// + // This is called whenever the run number changes // get DIRC geometry @@ -114,7 +122,7 @@ jerror_t DEventProcessor_dirc_tree::evnt(jana::JEventLoop* loop, uint64_t locEve auto locParticleComboStep = locPassedParticleCombos[icombo]->Get_ParticleComboStep(0); // calculate pi+pi- and K+K- invariant mass by taking the final state and subtracting the proton - DLorentzVector locInvP4 = fAnalysisUtilities->Calc_FinalStateP4(locReaction, locPassedParticleCombos[icombo], 0, false); + DLorentzVector locInvP4 = fAnalysisUtilities->Calc_FinalStateP4(locReaction, locPassedParticleCombos[icombo], 0, true); for(size_t parti=0; partiGet_NumFinalParticles(); parti++){ auto locParticle = locParticleComboStep->Get_FinalParticle(parti); if(locParticle->PID() == Proton) locInvP4 -= locParticle->lorentzMomentum(); From 8349d4e4b2a3c11d10349da35141e7fcc8ac3cdf Mon Sep 17 00:00:00 2001 From: gxproj7 gxproj7 Date: Wed, 24 Apr 2019 15:35:37 -0400 Subject: [PATCH 03/19] eliminate duplicated reactions --- .../dirc_tree/DEventProcessor_dirc_tree.cc | 26 +++---------------- 1 file changed, 3 insertions(+), 23 deletions(-) diff --git a/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc b/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc index 76ebd0b02..a7ea51e29 100644 --- a/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc +++ b/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc @@ -63,14 +63,6 @@ jerror_t DEventProcessor_dirc_tree::init(void) //------------------ jerror_t DEventProcessor_dirc_tree::brun(jana::JEventLoop* locEventLoop, int locRunNumber) { - ////////////////////////////////////////////////////////////////////////////// - // dapp and geom are not used but without it dirc_hits.so will no be loaded - // root [0] gSystem->Load("dirc_hits.so"); - // cling::DynamicLibraryManager::loadLibrary(): dirc_hits.so: undefined symbol: _ZTV25DEventSourceRESTGenerator - DApplication* dapp=dynamic_cast(locEventLoop->GetJApplication()); - DGeometry *geom = dapp->GetDGeometry(locRunNumber); - ////////////////////////////////////////////////////////////////////////////// - // This is called whenever the run number changes // get DIRC geometry @@ -113,6 +105,9 @@ jerror_t DEventProcessor_dirc_tree::evnt(jana::JEventLoop* loop, uint64_t locEve deque locPassedParticleCombos; locAnalysisResultsVector[loc_i]->Get_PassedParticleCombos(locPassedParticleCombos); const DReaction* locReaction = locAnalysisResultsVector[loc_i]->Get_Reaction(); + if(locReaction->Get_ReactionName() != "p2pi_dirc_tree" && + locReaction->Get_ReactionName() != "p2k_dirc_tree") continue; + std::vector previousInv; // loop over combos @@ -128,25 +123,10 @@ jerror_t DEventProcessor_dirc_tree::evnt(jana::JEventLoop* loop, uint64_t locEve if(locParticle->PID() == Proton) locInvP4 -= locParticle->lorentzMomentum(); } - int numIt=0; for(size_t parti=0; partiGet_NumFinalParticles(); parti++){ auto locParticle = locParticleComboStep->Get_FinalParticle(parti); // Get_FinalParticle_Measured(parti); if(locParticle->charge() == 0) continue; -/* - // we expect only rho or phi events: g,p->pi+,pi-,p g,p->K+,K-,p - if(locParticle->PID() == PiPlus || locParticle->PID() == PiMinus || locParticle->PID() == KPlus || locParticle->PID() == KMinus){ - locInvP4 += locParticle->lorentzMomentum(); - numIt++; - if(numIt==2){ - double tm = locInvP4.M(); - // do not store doubles - if(std::find_if(previousInv.begin(), previousInv.end(), [&tm](double m){return fabs(m-tm)<0.000000001;}) != previousInv.end()) { continue; } - previousInv.push_back(tm); - } - } -*/ - // get track auto locChargedTrack = static_cast(locParticleComboStep->Get_FinalParticle_SourceObject(parti)); auto locChargedTrackHypothesis = locChargedTrack->Get_Hypothesis(locParticle->PID()); From 3c4e17009ffd2916b58c605d00c7fb64c744e6e8 Mon Sep 17 00:00:00 2001 From: gxproj7 gxproj7 Date: Thu, 2 May 2019 18:26:00 -0400 Subject: [PATCH 04/19] update dirc_tree plugin and access to DIRCLut per-photon calculations --- src/libraries/DIRC/DDIRCLut.cc | 43 ++++++++++------- src/libraries/DIRC/DDIRCLut.h | 11 +++-- src/libraries/HDDM/DEventSourceREST.cc | 6 ++- src/libraries/PID/DParticleID.cc | 2 +- .../dirc_hists/DEventProcessor_dirc_hists.cc | 8 ++-- .../DCustomAction_dirc_reactions.cc | 47 ++++++++++++------- .../DCustomAction_dirc_reactions.h | 5 +- .../DReaction_factory_dirc_reactions.cc | 2 + .../dirc_tree/DEventProcessor_dirc_tree.cc | 2 + src/plugins/Analysis/dirc_tree/DrcEvent.cc | 2 +- src/plugins/Analysis/dirc_tree/DrcEvent.h | 3 ++ .../JEventProcessor_HLDetectorTiming.cc | 4 +- 12 files changed, 88 insertions(+), 47 deletions(-) diff --git a/src/libraries/DIRC/DDIRCLut.cc b/src/libraries/DIRC/DDIRCLut.cc index 9f8255555..8b5a71ada 100644 --- a/src/libraries/DIRC/DDIRCLut.cc +++ b/src/libraries/DIRC/DDIRCLut.cc @@ -26,9 +26,9 @@ DDIRCLut::DDIRCLut() gPARMS->SetDefaultParameter("DIRC:TRUTH_PIXELTIME",DIRC_TRUTH_PIXELTIME); // timing cuts for photons - DIRC_CUT_TDIFFD = 2; // direct cut in ns + DIRC_CUT_TDIFFD = 3; // direct cut in ns gPARMS->SetDefaultParameter("DIRC:CUT_TDIFFD",DIRC_CUT_TDIFFD); - DIRC_CUT_TDIFFR = 3; // reflected cut in ns + DIRC_CUT_TDIFFR = 4; // reflected cut in ns gPARMS->SetDefaultParameter("DIRC:CUT_TDIFFR",DIRC_CUT_TDIFFR); // Gives DeltaT = 0, but shouldn't it be v=20.3767 [cm/ns] for 1.47125 @@ -36,7 +36,7 @@ DDIRCLut::DDIRCLut() gPARMS->SetDefaultParameter("DIRC:LIGHT_V",DIRC_LIGHT_V); // sigma (thetaC for single photon) in radiams - DIRC_SIGMA_THETAC = 0.0085; + DIRC_SIGMA_THETAC = 0.01; gPARMS->SetDefaultParameter("DIRC:SIGMA_THETAC",DIRC_SIGMA_THETAC); // set PID for different passes in debuging histograms @@ -50,6 +50,9 @@ DDIRCLut::DDIRCLut() dCriticalAngle = asin(1.00028/1.47125); // n_quarzt = 1.47125; //(1.47125 <==> 390nm) dIndex = 1.473; + vector new_thetac(108); + dThetaC_offsets.push_back(new_thetac); dThetaC_offsets.push_back(new_thetac); + if(DIRC_DEBUG_HISTS) CreateDebugHistograms(); } @@ -64,6 +67,12 @@ bool DDIRCLut::brun(JEventLoop *loop) { loop->Get(locDIRCGeometry); dDIRCGeometry = locDIRCGeometry[0]; + // load constant tables + if (loop->GetCalib("/DIRC/North/thetac_offsets", dThetaC_offsets[0])) + jout << "Error loading /DIRC/North/thetac_offsets !" << endl; + if (loop->GetCalib("/DIRC/South/thetac_offsets", dThetaC_offsets[1])) + jout << "Error loading /DIRC/South/thetac_offsets !" << endl; + return true; } @@ -119,10 +128,8 @@ bool DDIRCLut::CreateDebugHistograms() { return true; } -bool DDIRCLut::CalcLUT(TVector3 locProjPos, TVector3 locProjMom, const vector locDIRCHits, double locFlightTime, Particle_t locPID, shared_ptr& locDIRCMatchParams, const vector locDIRCBarHits, map, vector >& locDIRCTrackMatchParams) const +bool DDIRCLut::CalcLUT(TVector3 locProjPos, TVector3 locProjMom, const vector locDIRCHits, double locFlightTime, double locMass, shared_ptr& locDIRCMatchParams, const vector locDIRCBarHits, map, vector >& locDIRCTrackMatchParams) const { - double locMass = ParticleMass(locPID); - // get bar and track position/momentum from extrapolation TVector3 momInBar = locProjMom; TVector3 posInBar = locProjPos; @@ -156,13 +163,13 @@ bool DDIRCLut::CalcLUT(TVector3 locProjPos, TVector3 locProjMom, const vector(); // initialize variables for LUT summary - double locAngle = CalcAngle(momInBar, locMass); - map locExpectedAngle = CalcExpectedAngles(momInBar); + double locAngle = CalcAngle(momInBar.Mag(), locMass); + map locExpectedAngle = CalcExpectedAngles(momInBar.Mag()); map logLikelihoodSum; Particle_t locHypothesisPID = PiPlus; for(uint loc_i = 0; loc_i> DDIRCLut::CalcPhoton(const DDIRCPmtHit *locDIRCHit, int box = (bar < 24) ? 1 : 0; if((box == 0 && channel < dMaxChannels) || (box == 1 && channel >= dMaxChannels)) return locDIRCPhotons; + + int pmt = channel/64; + double thetac_offset = dThetaC_offsets[box][pmt]; // use hit time to determine if reflected or not double hitTime = locDIRCHit->t - locFlightTime; @@ -247,7 +257,7 @@ vector> DDIRCLut::CalcPhoton(const DDIRCPmtHit *locDIRCHit, } // needs to be X dependent choice for reflection cut (from CCDB?) - bool reflected = hitTime>0; // try all photons as reflected for now + bool reflected = hitTime>35; // try only some photons as reflected for now // get position along bar for calculated time double radiatorL = dDIRCGeometry->GetBarLength(bar); @@ -292,7 +302,8 @@ vector> DDIRCLut::CalcPhoton(const DDIRCPmtHit *locDIRCHit, luttheta = dir.Angle(TVector3(-1,0,0)); if(luttheta > TMath::PiOver2()) luttheta = TMath::Pi()-luttheta; - tangle = momInBar.Angle(dir);//-0.002; //correction + tangle = momInBar.Angle(dir); + tangle -= thetac_offset; //correction double bartime = lenz/cos(luttheta)/DIRC_LIGHT_V; double totalTime = bartime+evtime; @@ -336,7 +347,7 @@ vector> DDIRCLut::CalcPhoton(const DDIRCPmtHit *locDIRCHit, } // remove photon candidates not used in likelihood - if(fabs(tangle-0.5*(locExpectedAngle[PiPlus]+locExpectedAngle[KPlus]))>0.05) continue; + if(fabs(tangle-0.5*(locExpectedAngle[PiPlus]+locExpectedAngle[KPlus]))>0.03) continue; // save good photons to matched list isGood = true; @@ -376,15 +387,15 @@ double DDIRCLut::CalcLikelihood(double locExpectedThetaC, double locThetaC) cons return locLikelihood; } -double DDIRCLut::CalcAngle(TVector3 momInBar, double locMass) const { - return acos(sqrt(momInBar.Mag()*momInBar.Mag() + locMass*locMass)/momInBar.Mag()/dIndex); +double DDIRCLut::CalcAngle(double locP, double locMass) const { + return acos(sqrt(locP*locP + locMass*locMass)/locP/dIndex); } -map DDIRCLut::CalcExpectedAngles(TVector3 momInBar) const { +map DDIRCLut::CalcExpectedAngles(double locP) const { map locExpectedAngles; for(uint loc_i = 0; loc_i dirc_lut_constants_t; + class DDIRCLut: public JObject { public: @@ -34,14 +37,16 @@ class DDIRCLut: public JObject { DDIRCLut(); ~DDIRCLut(){}; + vector dThetaC_offsets; + bool brun(JEventLoop *loop); bool CreateDebugHistograms(); - bool CalcLUT(TVector3 locProjPos, TVector3 locProjMom, const vector locDIRCHits, double locFlightTime, Particle_t locPID, shared_ptr& locDIRCMatchParams, const vector locDIRCBarHits, map, vector >& locDIRCTrackMatchParams) const; + bool CalcLUT(TVector3 locProjPos, TVector3 locProjMom, const vector locDIRCHits, double locFlightTime, double locMass, shared_ptr& locDIRCMatchParams, const vector locDIRCBarHits, map, vector >& locDIRCTrackMatchParams) const; vector> CalcPhoton(const DDIRCPmtHit *locDIRCHit, double locFlightTime, TVector3 posInBar, TVector3 momInBar, map locExpectedAngle, double locAngle, Particle_t locPID, map &logLikelihoodSum, int &nPhotonsThetaC, double &meanThetaC, double &meanDeltaT, bool &isGood) const; vector> CalcPhoton(const DDIRCPmtHit *locDIRCHit, double locFlightTime, TVector3 posInBar, TVector3 momInBar, map locExpectedAngle, double locAngle, Particle_t locPID, map &logLikelihoodSum) const; double CalcLikelihood(double locExpectedThetaC, double locThetaC) const; - double CalcAngle(TVector3 momInBar, double locMass) const; - map CalcExpectedAngles(TVector3 momInBar) const; + double CalcAngle(double locP, double locMass) const; + map CalcExpectedAngles(double locP) const; private: DApplication *dapp; diff --git a/src/libraries/HDDM/DEventSourceREST.cc b/src/libraries/HDDM/DEventSourceREST.cc index 438772275..6c511f4b8 100644 --- a/src/libraries/HDDM/DEventSourceREST.cc +++ b/src/libraries/HDDM/DEventSourceREST.cc @@ -1250,6 +1250,10 @@ jerror_t DEventSourceREST::Extract_DDetectorMatches(JEventLoop* locEventLoop, hd for(; dircIter != dircList.end(); ++dircIter) { size_t locTrackIndex = dircIter->getTrack(); + if(locTrackIndex > locTrackTimeBasedVector.size()) continue; + + auto locTrackTimeBased = locTrackTimeBasedVector[locTrackIndex]; + if( !locTrackTimeBased ) continue; auto locDIRCMatchParams = std::make_shared(); map ,vector > locDIRCTrackMatchParams; @@ -1260,7 +1264,7 @@ jerror_t DEventSourceREST::Extract_DDetectorMatches(JEventLoop* locEventLoop, hd TVector3 locProjMom(dircIter->getPx(),dircIter->getPy(),dircIter->getPz()); double locFlightTime = dircIter->getT(); - if( locParticleID->Get_DIRCLut()->CalcLUT(locProjPos, locProjMom, locDIRCHits, locFlightTime, locTrackTimeBasedVector[locTrackIndex]->PID(), locDIRCMatchParams, locDIRCBarHits, locDIRCTrackMatchParams) ) + if( locParticleID->Get_DIRCLut()->CalcLUT(locProjPos, locProjMom, locDIRCHits, locFlightTime, locTrackTimeBased->mass(), locDIRCMatchParams, locDIRCBarHits, locDIRCTrackMatchParams) ) locDetectorMatches->Add_Match(locTrackTimeBasedVector[locTrackIndex], std::const_pointer_cast(locDIRCMatchParams)); } else { diff --git a/src/libraries/PID/DParticleID.cc b/src/libraries/PID/DParticleID.cc index e68d9342f..1b6e7ca90 100644 --- a/src/libraries/PID/DParticleID.cc +++ b/src/libraries/PID/DParticleID.cc @@ -1787,7 +1787,7 @@ bool DParticleID::Cut_MatchDIRC(const vector &ext } // Calculate DIRC LUT - return dDIRCLut->CalcLUT(locProjPos, locProjMom, locDIRCHits, locFlightTime, locPID, locDIRCMatchParams, locDIRCBarHits, locDIRCTrackMatchParams); + return dDIRCLut->CalcLUT(locProjPos, locProjMom, locDIRCHits, locFlightTime, ParticleMass(locPID), locDIRCMatchParams, locDIRCBarHits, locDIRCTrackMatchParams); } diff --git a/src/plugins/Analysis/dirc_hists/DEventProcessor_dirc_hists.cc b/src/plugins/Analysis/dirc_hists/DEventProcessor_dirc_hists.cc index 57637d843..8a4500c14 100644 --- a/src/plugins/Analysis/dirc_hists/DEventProcessor_dirc_hists.cc +++ b/src/plugins/Analysis/dirc_hists/DEventProcessor_dirc_hists.cc @@ -177,8 +177,8 @@ jerror_t DEventProcessor_dirc_hists::evnt(JEventLoop *loop, uint64_t eventnumber hExtrapolatedBarHitXY[locPID]->Fill(posInBar.X(), posInBar.Y()); japp->RootFillUnLock(this); //RELEASE ROOT FILL LOCK - double locAngle = dDIRCLut->CalcAngle(momInBar, locMass); - map locExpectedAngle = dDIRCLut->CalcExpectedAngles(momInBar); + double locAngle = dDIRCLut->CalcAngle(locP, locMass); + map locExpectedAngle = dDIRCLut->CalcExpectedAngles(locP); // get map of DIRCMatches to PMT hits map, vector > locDIRCTrackMatchParamsMap; @@ -192,7 +192,7 @@ jerror_t DEventProcessor_dirc_hists::evnt(JEventLoop *loop, uint64_t eventnumber vector> locDIRCPhotons = dDIRCLut->CalcPhoton(locDIRCPmtHits[loc_i], locExtrapolatedTime, posInBar, momInBar, locExpectedAngle, locAngle, locPID, logLikelihoodSum); double locHitTime = locDIRCPmtHits[loc_i]->t - locExtrapolatedTime; int locChannel = locDIRCPmtHits[loc_i]->ch%dMaxChannels; - if(locHitTime > 0 && locHitTime < 100) locPhotonInclusive++; + if(locHitTime > 0 && locHitTime < 150) locPhotonInclusive++; int pixel_row = dDIRCGeometry->GetPixelRow(locChannel); int pixel_col = dDIRCGeometry->GetPixelColumn(locChannel); @@ -230,7 +230,7 @@ jerror_t DEventProcessor_dirc_hists::evnt(JEventLoop *loop, uint64_t eventnumber } // fill histograms for candidate photons in timing cut - if(fabs(locDeltaT) < 100.0) { + if(fabs(locDeltaT) < 5.0) { hThetaC[locPID]->Fill(locThetaC); hDeltaThetaC[locPID]->Fill(locThetaC-locExpectedThetaC); hDeltaThetaCVsP[locPID]->Fill(momInBar.Mag(), locThetaC-locExpectedThetaC); diff --git a/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.cc b/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.cc index 27f4ac642..f7b749a05 100644 --- a/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.cc +++ b/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.cc @@ -61,6 +61,8 @@ void DCustomAction_dirc_reactions::Initialize(JEventLoop* locEventLoop) hThetaCVsP = GetOrCreate_Histogram(Form("hThetaCVsP_%s",locParticleName.data()), Form("cherenkov angle vs. momentum; p (GeV/c); %s #theta_{C} [rad]", locParticleROOTName.data()), 120, 0.0, 12.0, 250, 0.6, 1.0); hDeltaThetaCVsP = GetOrCreate_Histogram(Form("hDeltaThetaCVsP_%s",locParticleName.data()), Form("cherenkov angle vs. momentum; p (GeV/c); %s #Delta#theta_{C} [rad]", locParticleROOTName.data()), 120, 0.0, 12.0, 200,-0.2,0.2); + hDeltaThetaCVsChannel = GetOrCreate_Histogram(Form("hDeltaThetaCVsChannel_%s",locParticleName.data()), Form("cherenkov angle vs. channel; channel ID; %s #Delta#theta_{C} [rad]", locParticleROOTName.data()), 6912, 0, 6912, 200,-0.2,0.2); + hLikelihoodDiffVsP = GetOrCreate_Histogram(Form("hLikelihoodDiffVsP_%s",locParticleName.data()), Form("; p (GeV/c); %s", locLikelihoodName.data()), 120, 0.0, 12.0, 100, -200, 200); hReactionLikelihoodDiffVsP = GetOrCreate_Histogram(Form("hReactionLikelihoodDiffVsP_%s",locParticleName.data()), Form("; p (GeV/c); %s", locLikelihoodName.data()), 120, 0.0, 12.0, 100, -200, 200); @@ -82,19 +84,20 @@ void DCustomAction_dirc_reactions::Initialize(JEventLoop* locEventLoop) if(locBinVect.Theta()*180/TMath::Pi() > 12.) continue; - hDiffMap[locBar][locXbin] = GetOrCreate_Histogram(Form("hDiff_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; %s t_{calc}-t_{measured} [ns]; entries [#]", locBar,xbin_min,xbin_max,locParticleROOTName.data()), 200,-100,100); - hHitTimeMap[locBar][locXbin] = GetOrCreate_Histogram(Form("hHitTimeMap_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; %s t_{measured} [ns]; entries [#]", locBar,xbin_min,xbin_max,locParticleROOTName.data()), 200,0,200); + hDiffMap[locBar][locXbin] = GetOrCreate_Histogram(Form("hDiff_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; %s t_{calc}-t_{measured} [ns]; entries [#]", locBar,xbin_min,xbin_max,locParticleROOTName.data()), 100,-20,20); + hHitTimeMap[locBar][locXbin] = GetOrCreate_Histogram(Form("hHitTimeMap_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; %s t_{measured} [ns]; entries [#]", locBar,xbin_min,xbin_max,locParticleROOTName.data()), 100,0,150); hNphCMap[locBar][locXbin] = GetOrCreate_Histogram(Form("hNphC_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f] # photons; %s # photons", locBar,xbin_min,xbin_max,locParticleROOTName.data()), 80, 0, 80); hNphCMapSlot4[locBar][locXbin] = GetOrCreate_Histogram(Form("hNphCSlot4_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f] # photons; %s # photons", locBar,xbin_min,xbin_max,locParticleROOTName.data()), 80, 0, 80); hNphCMapSlot5[locBar][locXbin] = GetOrCreate_Histogram(Form("hNphCSlot5_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f] # photons; %s # photons", locBar,xbin_min,xbin_max,locParticleROOTName.data()), 80, 0, 80); + hThetaCVsPMap[locBar][locXbin] = GetOrCreate_Histogram(Form("hThetaCVsP_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f] cherenkov angle vs. momentum; p (GeV/c); %s #theta_{C} [rad]", locBar,xbin_min,xbin_max,locParticleROOTName.data()), 60, 0.0, 12.0, 250,0.6,1.0); hDeltaThetaCVsPMap[locBar][locXbin] = GetOrCreate_Histogram(Form("hDeltaThetaCVsP_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f] cherenkov angle vs. momentum; p (GeV/c); %s #Delta#theta_{C} [rad]", locBar,xbin_min,xbin_max,locParticleROOTName.data()), 60, 0.0, 12.0, 60,-0.15,0.15); hReactionLikelihoodDiffVsPMap[locBar][locXbin] = GetOrCreate_Histogram(Form("hReactionLikelihoodDiffVsP_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; p (GeV/c); %s", locBar,xbin_min,xbin_max,locLikelihoodName.data()), 60, 0.0, 12.0, 50, -200, 200); - + hPixelHitMap3D[locBar][locXbin] = GetOrCreate_Histogram(Form("hPixelHit3D_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; pixel rows; pixel columns; hit time", locBar,xbin_min,xbin_max), 144, -0.5, 143.5, 48, -0.5, 47.5, 50, 0, 100); hPixelHitMap[locBar][locXbin] = GetOrCreate_Histogram(Form("hPixelHit_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; pixel rows; pixel columns", locBar,xbin_min,xbin_max), 144, -0.5, 143.5, 48, -0.5, 47.5); hPixelHitMapReflected[locBar][locXbin] = GetOrCreate_Histogram(Form("hPixelHitReflected_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; pixel rows; pixel columns", locBar,xbin_min,xbin_max), 144, -0.5, 143.5, 48, -0.5, 47.5); - hHitTimeMapAll[locBar][locXbin] = GetOrCreate_Histogram(Form("hHitTimeMapAll_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; %s t_{measured} [ns]; entries [#]", locBar,xbin_min,xbin_max,locParticleROOTName.data()), 200,0,200); + hHitTimeMapAll[locBar][locXbin] = GetOrCreate_Histogram(Form("hHitTimeMapAll_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; %s t_{measured} [ns]; entries [#]", locBar,xbin_min,xbin_max,locParticleROOTName.data()), 100,0,150); hPixelHitMapAll[locBar][locXbin] = GetOrCreate_Histogram(Form("hPixelHitAll_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; pixel rows; pixel columns", locBar,xbin_min,xbin_max), 144, -0.5, 143.5, 48, -0.5, 47.5); hPixelHitMapAllReflected[locBar][locXbin] = GetOrCreate_Histogram(Form("hPixelHitAllReflected_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; pixel rows; pixel columns", locBar,xbin_min,xbin_max), 144, -0.5, 143.5, 48, -0.5, 47.5); //hPixelHitTimeMap[locBar][locXbin] = GetOrCreate_Histogram(Form("hPixelHitTime_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; Pixel Hit Channel; Pixel Hit t [ns]", locBar,xbin_min,xbin_max), 6912, 0, 6912, 50, 0, 100); @@ -178,8 +181,8 @@ bool DCustomAction_dirc_reactions::Perform_Action(JEventLoop* locEventLoop, cons hExtrapolatedBarHitXY->Fill(posInBar.X(), posInBar.Y()); Unlock_Action(); //RELEASE ROOT LOCK!! - double locAngle = dDIRCLut->CalcAngle(momInBar, locMass); - map locExpectedAngle = dDIRCLut->CalcExpectedAngles(momInBar); + double locAngle = dDIRCLut->CalcAngle(locP, locMass); + map locExpectedAngle = dDIRCLut->CalcExpectedAngles(locP); // get map of DIRCMatches to PMT hits map, vector > locDIRCTrackMatchParamsMap; @@ -208,7 +211,7 @@ bool DCustomAction_dirc_reactions::Perform_Action(JEventLoop* locEventLoop, cons int pixel_col = dDIRCGeometry->GetPixelColumn(locChannel); // fill histograms for candidate photons in timing cut - if(locHitTime > 0 && locHitTime < 100.0) { + if(locHitTime > 0 && locHitTime < 150.0) { if((pixel_row < 64 && pixel_col < 24) || (pixel_row < 32 && pixel_col > 23)) locNPhCSlot4++; else @@ -217,6 +220,7 @@ bool DCustomAction_dirc_reactions::Perform_Action(JEventLoop* locEventLoop, cons Lock_Action(); //ACQUIRE ROOT LOCK!! if(DIRC_FILL_BAR_MAP && locP > 4.) { //hPixelHitTimeMap[locBar][locXbin]->Fill(locChannel, locHitTime); + hPixelHitMap3D[locBar][locXbin]->Fill(pixel_row, pixel_col, locHitTime); if(locHitTime < 38.) hPixelHitMapAll[locBar][locXbin]->Fill(pixel_row, pixel_col); else @@ -234,21 +238,28 @@ bool DCustomAction_dirc_reactions::Perform_Action(JEventLoop* locEventLoop, cons double locDeltaT = locDIRCPhotons[loc_j].first - locHitTime; double locThetaC = locDIRCPhotons[loc_j].second; - Lock_Action(); //ACQUIRE ROOT LOCK!! - hDiff->Fill(locDeltaT); - if(DIRC_FILL_BAR_MAP) { - hDiffMap[locBar][locXbin]->Fill(locDeltaT); - hHitTimeMap[locBar][locXbin]->Fill(locHitTime); - } - Unlock_Action(); //RELEASE ROOT LOCK!! - + if(fabs(locThetaC-locExpectedThetaC)<0.02) { + Lock_Action(); //ACQUIRE ROOT LOCK!! + hDiff->Fill(locDeltaT); + if(DIRC_FILL_BAR_MAP) { + hDiffMap[locBar][locXbin]->Fill(locDeltaT); + hHitTimeMap[locBar][locXbin]->Fill(locHitTime); + } + Unlock_Action(); //RELEASE ROOT LOCK!! + } + // fill histograms for candidate photons in timing cut - if(fabs(locDeltaT) < 100.0) { + if(fabs(locDeltaT) < 2.0) { Lock_Action(); //ACQUIRE ROOT LOCK!! hThetaC->Fill(locThetaC); hDeltaThetaC->Fill(locThetaC-locExpectedThetaC); + hThetaCVsP->Fill(momInBar.Mag(), locThetaC); hDeltaThetaCVsP->Fill(momInBar.Mag(), locThetaC-locExpectedThetaC); - if(DIRC_FILL_BAR_MAP) hDeltaThetaCVsPMap[locBar][locXbin]->Fill(momInBar.Mag(), locThetaC-locExpectedThetaC); + hDeltaThetaCVsChannel->Fill(locChannel, locThetaC-locExpectedThetaC); + if(DIRC_FILL_BAR_MAP) { + hDeltaThetaCVsPMap[locBar][locXbin]->Fill(momInBar.Mag(), locThetaC-locExpectedThetaC); + hThetaCVsPMap[locBar][locXbin]->Fill(momInBar.Mag(), locThetaC); + } Unlock_Action(); //RELEASE ROOT LOCK!! if(std::find(locUsedPixel.begin(), locUsedPixel.end(), locChannel) != locUsedPixel.end()) @@ -272,7 +283,7 @@ bool DCustomAction_dirc_reactions::Perform_Action(JEventLoop* locEventLoop, cons // fill histograms with per-track quantities Lock_Action(); //ACQUIRE ROOT LOCK!! hNphC->Fill(locDIRCMatchParams->dNPhotons); - hThetaCVsP->Fill(momInBar.Mag(), locDIRCMatchParams->dThetaC); + //hThetaCVsP->Fill(momInBar.Mag(), locDIRCMatchParams->dThetaC); if(DIRC_FILL_BAR_MAP) { hNphCMap[locBar][locXbin]->Fill(locDIRCMatchParams->dNPhotons); hNphCMapSlot4[locBar][locXbin]->Fill(locNPhCSlot4); diff --git a/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.h b/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.h index 79bcda7c8..afc611515 100644 --- a/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.h +++ b/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.h @@ -10,6 +10,7 @@ #include #include "TH1.h" +#include "TH3.h" #include "TLorentzRotation.h" #include "JANA/JEventLoop.h" @@ -58,14 +59,16 @@ class DCustomAction_dirc_reactions : public DAnalysisAction //Store any histograms as member variables here TH1I *hDiff, *hNphC, *hThetaC, *hDeltaThetaC, *hLikelihood, *hLikelihoodDiff; TH2I *hThetaCVsP, *hDeltaThetaCVsP, *hLikelihoodDiffVsP, *hReactionLikelihoodDiffVsP; + TH2I *hDeltaThetaCVsChannel; TH2I *hExtrapolatedBarHitXY, *hExtrapolatedBarHitXY_PreCut; TH1I *hDiffMap[48][40], *hHitTimeMap[48][40], *hHitTimeMapAll[48][40]; TH1I *hNphCMap[48][40], *hNphCMapSlot4[48][40], *hNphCMapSlot5[48][40]; + TH3S *hPixelHitMap3D[48][40]; TH2S *hPixelHitMap[48][40], *hPixelHitMapReflected[48][40]; TH2S *hPixelHitMapAll[48][40], *hPixelHitMapAllReflected[48][40]; //TH2I *hPixelHitTimeMap[48][40]; - TH2I *hDeltaThetaCVsPMap[48][40], *hReactionLikelihoodDiffVsPMap[48][40]; + TH2I *hThetaCVsPMap[48][40], *hDeltaThetaCVsPMap[48][40], *hReactionLikelihoodDiffVsPMap[48][40]; }; #endif // _DCustomAction_dirc_reactions_ diff --git a/src/plugins/Analysis/dirc_reactions/DReaction_factory_dirc_reactions.cc b/src/plugins/Analysis/dirc_reactions/DReaction_factory_dirc_reactions.cc index 24d17a712..871fd253e 100644 --- a/src/plugins/Analysis/dirc_reactions/DReaction_factory_dirc_reactions.cc +++ b/src/plugins/Analysis/dirc_reactions/DReaction_factory_dirc_reactions.cc @@ -72,6 +72,7 @@ jerror_t DReaction_factory_dirc_reactions::evnt(JEventLoop* locEventLoop, uint64 locReaction->Add_AnalysisAction(new DHistogramAction_MissingMassSquared(locReaction, false, 1000, -0.1, 0.1, "PostKinFitCut")); locReaction->Add_AnalysisAction(new DHistogramAction_InvariantMass(locReaction, 0, locRhoPIDs, false, 900, 0.3, 1.2, "Rho_PostKinFitCut")); locReaction->Add_AnalysisAction(new DHistogramAction_InvariantMass(locReaction, 0, locRhoPIDs, true, 900, 0.3, 1.2, "Rho_KinFit_PostKinFitCut")); + locReaction->Add_AnalysisAction(new DCutAction_InvariantMass(locReaction, 0, locRhoPIDs, true, 0.5, 1.0)); // Custom histograms for DIRC //locReaction->Add_AnalysisAction(new DCustomAction_dirc_reactions(locReaction, false, PiPlus, 0, 1, "PiPlus_DIRC")); @@ -132,6 +133,7 @@ jerror_t DReaction_factory_dirc_reactions::evnt(JEventLoop* locEventLoop, uint64 locReaction->Add_AnalysisAction(new DHistogramAction_MissingMassSquared(locReaction, false, 1000, -0.1, 0.1, "PostKinFitCut")); locReaction->Add_AnalysisAction(new DHistogramAction_InvariantMass(locReaction, 0, locPhiPIDs, false, 500, 0.9, 1.4, "Phi_PostKinFitCut")); locReaction->Add_AnalysisAction(new DHistogramAction_InvariantMass(locReaction, 0, locPhiPIDs, true, 500, 0.9, 1.4, "Phi_KinFit_PostKinFitCut")); + locReaction->Add_AnalysisAction(new DCutAction_InvariantMass(locReaction, 0, locPhiPIDs, true, 0.95, 1.05)); // Custom histograms for DIRC //locReaction->Add_AnalysisAction(new DCustomAction_dirc_reactions(locReaction, false, KPlus, 0, 1, "KPlus_DIRC")); diff --git a/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc b/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc index a7ea51e29..7093f44e2 100644 --- a/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc +++ b/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc @@ -143,6 +143,7 @@ jerror_t DEventProcessor_dirc_tree::evnt(jana::JEventLoop* loop, uint64_t locEve bool foundTOF = dParticleID->Get_BestTOFMatchParams(locTrackTimeBased, locDetectorMatches, locTOFHitMatchParams); if(!foundTOF || locTOFHitMatchParams->Get_DistanceToTrack() > 20.0) continue; double toftrackdist = locTOFHitMatchParams->Get_DistanceToTrack(); + double toftrackdeltat = locChargedTrackHypothesis->Get_TimeAtPOCAToVertex() - locChargedTrackHypothesis->t0(); Particle_t locPID = locTrackTimeBased->PID(); // get DIRC match parameters (contains LUT information) @@ -165,6 +166,7 @@ jerror_t DEventProcessor_dirc_tree::evnt(jana::JEventLoop* loop, uint64_t locEve fEvent->SetPosition(TVector3(posInBar.X(), posInBar.Y(), posInBar.Z())); fEvent->SetDcHits(locDCHits); fEvent->SetTofTrackDist(toftrackdist); + fEvent->SetTofTrackDeltaT(toftrackdeltat); fEvent->SetInvMass(locInvP4.M()); fEvent->SetMissMass(locMissingP4.M2()); fEvent->SetChiSq(chisq); diff --git a/src/plugins/Analysis/dirc_tree/DrcEvent.cc b/src/plugins/Analysis/dirc_tree/DrcEvent.cc index 745062fc1..f280fb256 100644 --- a/src/plugins/Analysis/dirc_tree/DrcEvent.cc +++ b/src/plugins/Analysis/dirc_tree/DrcEvent.cc @@ -4,7 +4,7 @@ ClassImp(DrcEvent) // // ----- Default constructor ------------------------------------------- DrcEvent::DrcEvent(): fId(-1),fType(0),fPdg(0),fParent(0),fDcHits(0), - fTime(-100),fInvMass(-100),fMissMass(-100),fChiSq(-100),fTofTrackDist(-100),fHitSize(0),fTest1(0),fTest2(0){ + fTime(-100),fInvMass(-100),fMissMass(-100),fChiSq(-100),fTofTrackDist(-100),fTofTrackDeltaT(-100),fHitSize(0),fTest1(0),fTest2(0){ } void DrcEvent::AddHit(DrcHit hit){ diff --git a/src/plugins/Analysis/dirc_tree/DrcEvent.h b/src/plugins/Analysis/dirc_tree/DrcEvent.h index 64f5bb38e..5fd949c0a 100644 --- a/src/plugins/Analysis/dirc_tree/DrcEvent.h +++ b/src/plugins/Analysis/dirc_tree/DrcEvent.h @@ -33,6 +33,7 @@ class DrcEvent: public TObject { Double_t GetMissMass() const { return fMissMass; } Double_t GetChiSq() const { return fChiSq; } Double_t GetTofTrackDist() const { return fTofTrackDist; } + Double_t GetTofTrackDeltaT() const { return fTofTrackDeltaT; } Int_t GetPdg() const { return fPdg; } Int_t GetParent() const { return fParent; } @@ -51,6 +52,7 @@ class DrcEvent: public TObject { void SetMissMass(Double_t val){ fMissMass=val; } void SetChiSq(Double_t val) { fChiSq=val; } void SetTofTrackDist(Double_t val) { fTofTrackDist=val; } + void SetTofTrackDeltaT(Double_t val) { fTofTrackDeltaT=val; } void SetPdg(Int_t val) { fPdg = val; } void SetParent(Int_t val) { fParent = val; } @@ -71,6 +73,7 @@ class DrcEvent: public TObject { Double_t fMissMass; Double_t fChiSq; Double_t fTofTrackDist; + Double_t fTofTrackDeltaT; Int_t fHitSize; std::vector fHitArray; diff --git a/src/plugins/Calibration/HLDetectorTiming/JEventProcessor_HLDetectorTiming.cc b/src/plugins/Calibration/HLDetectorTiming/JEventProcessor_HLDetectorTiming.cc index 11669c751..e94934c13 100644 --- a/src/plugins/Calibration/HLDetectorTiming/JEventProcessor_HLDetectorTiming.cc +++ b/src/plugins/Calibration/HLDetectorTiming/JEventProcessor_HLDetectorTiming.cc @@ -1085,8 +1085,8 @@ jerror_t JEventProcessor_HLDetectorTiming::evnt(JEventLoop *loop, uint64_t event Particle_t locPID = locTrackTimeBased->PID(); double locMass = ParticleMass(locPID); - double locAngle = dDIRCLut->CalcAngle(momInBar, locMass); - map locExpectedAngle = dDIRCLut->CalcExpectedAngles(momInBar); + double locAngle = dDIRCLut->CalcAngle(momInBar.Mag(), locMass); + map locExpectedAngle = dDIRCLut->CalcExpectedAngles(momInBar.Mag()); // get map of DIRCMatches to PMT hits map, vector > locDIRCTrackMatchParamsMap; From d1efa4380abf7e3de613d203c3e7f8a712e23ade Mon Sep 17 00:00:00 2001 From: Justin Stevens Date: Tue, 7 May 2019 07:31:05 -0400 Subject: [PATCH 05/19] fix DIRC detector matches in DChargedTrackHypothesis --- src/libraries/ANALYSIS/DEventWriterROOT.cc | 21 +++++++++++----- src/libraries/DIRC/DDIRCLutReader.cc | 2 +- src/libraries/PID/DChargedTrackHypothesis.h | 1 + .../DCustomAction_dirc_reactions.cc | 24 ++++++++++++++++++- 4 files changed, 40 insertions(+), 8 deletions(-) diff --git a/src/libraries/ANALYSIS/DEventWriterROOT.cc b/src/libraries/ANALYSIS/DEventWriterROOT.cc index 510ab4aa2..9795bcbcd 100644 --- a/src/libraries/ANALYSIS/DEventWriterROOT.cc +++ b/src/libraries/ANALYSIS/DEventWriterROOT.cc @@ -1624,17 +1624,26 @@ void DEventWriterROOT::Fill_ChargedHypo(DTreeFillData* locTreeFillData, unsigned // DIRC if(DIRC_OUTPUT) { - int locDIRCNumPhotons = (locChargedTrackHypothesis->Get_DIRCMatchParams() != NULL) ? locChargedTrackHypothesis->Get_DIRCMatchParams()->dNPhotons : 0; + int locDIRCNumPhotons = 0; + double locDIRCThetaC = 999.; + double locDIRCLele = 999.; + double locDIRCLpi = 999.; + double locDIRCLk = 999.; + double locDIRCLp = 999.; + auto locDIRCMatchParams = locChargedTrackHypothesis->Get_DIRCMatchParams(); + if(locDIRCMatchParams != NULL) { + locDIRCNumPhotons = locDIRCMatchParams->dNPhotons; + locDIRCThetaC = locDIRCMatchParams->dThetaC; + locDIRCLele = locDIRCMatchParams->dLikelihoodElectron; + locDIRCLpi = locDIRCMatchParams->dLikelihoodPion; + locDIRCLk = locDIRCMatchParams->dLikelihoodKaon; + locDIRCLp = locDIRCMatchParams->dLikelihoodProton; + } locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "NumPhotons_DIRC"), locDIRCNumPhotons, locArrayIndex); - double locDIRCThetaC = (locChargedTrackHypothesis->Get_DIRCMatchParams() != NULL) ? locChargedTrackHypothesis->Get_DIRCMatchParams()->dThetaC : 0.0; locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "ThetaC_DIRC"), locDIRCThetaC, locArrayIndex); - double locDIRCLele = (locChargedTrackHypothesis->Get_DIRCMatchParams() != NULL) ? locChargedTrackHypothesis->Get_DIRCMatchParams()->dLikelihoodElectron : 0.0; locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "Lele_DIRC"), locDIRCLele, locArrayIndex); - double locDIRCLpi = (locChargedTrackHypothesis->Get_DIRCMatchParams() != NULL) ? locChargedTrackHypothesis->Get_DIRCMatchParams()->dLikelihoodPion : 0.0; locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "Lpi_DIRC"), locDIRCLpi, locArrayIndex); - double locDIRCLk = (locChargedTrackHypothesis->Get_DIRCMatchParams() != NULL) ? locChargedTrackHypothesis->Get_DIRCMatchParams()->dLikelihoodKaon : 0.0; locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "Lk_DIRC"), locDIRCLk, locArrayIndex); - double locDIRCLp = (locChargedTrackHypothesis->Get_DIRCMatchParams() != NULL) ? locChargedTrackHypothesis->Get_DIRCMatchParams()->dLikelihoodProton : 0.0; locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "Lp_DIRC"), locDIRCLp, locArrayIndex); } } diff --git a/src/libraries/DIRC/DDIRCLutReader.cc b/src/libraries/DIRC/DDIRCLutReader.cc index 5a4547e5e..d4ab6b709 100644 --- a/src/libraries/DIRC/DDIRCLutReader.cc +++ b/src/libraries/DIRC/DDIRCLutReader.cc @@ -29,7 +29,7 @@ DDIRCLutReader::DDIRCLutReader(JApplication *japp, unsigned int run_number) jcalib = japp->GetJCalibration(run_number); if(jcalib->GetCalib("/DIRC/LUT/lut_map", lut_map_name)) jout << "Can't find requested /DIRC/LUT/lut_map in CCDB for this run!" << endl; - else if(lut_map_name.find("map_name") != lut_map_name.end() && lut_map_name["map_name"] != "None") { + else if(lut_map_name.find("map_name") != lut_map_name.end() && lut_map_name["map_name"] != "None" && lut_file.empty()) { jresman = japp->GetJResourceManager(run_number); lut_file = jresman->GetResource(lut_map_name["map_name"]); } diff --git a/src/libraries/PID/DChargedTrackHypothesis.h b/src/libraries/PID/DChargedTrackHypothesis.h index d1d8221a1..10c6bcd79 100644 --- a/src/libraries/PID/DChargedTrackHypothesis.h +++ b/src/libraries/PID/DChargedTrackHypothesis.h @@ -338,6 +338,7 @@ inline void DChargedTrackHypothesis::DTrackingInfo::Reset(void) dTOFHitMatchParams = nullptr; dBCALShowerMatchParams = nullptr; dFCALShowerMatchParams = nullptr; + dDIRCMatchParams = nullptr; } #endif // _DChargedTrackHypothesis_ diff --git a/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.cc b/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.cc index f7b749a05..cfd7add08 100644 --- a/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.cc +++ b/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.cc @@ -167,7 +167,29 @@ bool DCustomAction_dirc_reactions::Perform_Action(JEventLoop* locEventLoop, cons DVector3 momInBar = locDIRCMatchParams->dExtrapolatedMom; double locExpectedThetaC = locDIRCMatchParams->dExpectedThetaC; double locExtrapolatedTime = locDIRCMatchParams->dExtrapolatedTime; - + + /////////////////////////////////////////// + // option to cheat and use truth bar hit // + /////////////////////////////////////////// + if(DIRC_TRUTH_BARHIT && locDIRCBarHits.size() > 0) { + + TVector3 bestMatchPos, bestMatchMom; + double bestMatchDist = 999.; + for(int i=0; i<(int)locDIRCBarHits.size(); i++) { + TVector3 locDIRCBarHitPos(locDIRCBarHits[i]->x, locDIRCBarHits[i]->y, locDIRCBarHits[i]->z); + TVector3 locDIRCBarHitMom(locDIRCBarHits[i]->px, locDIRCBarHits[i]->py, locDIRCBarHits[i]->pz); + if((posInBar - locDIRCBarHitPos).Mag() < bestMatchDist) { + bestMatchDist = (posInBar - locDIRCBarHitPos).Mag(); + bestMatchPos = locDIRCBarHitPos; + bestMatchMom = locDIRCBarHitMom; + } + } + + momInBar = bestMatchMom; + posInBar = bestMatchPos; + } + + // get binning for histograms int locBar = dDIRCGeometry->GetBar(posInBar.Y()); int locXbin = (int)(posInBar.X()/5.0) + 19; From 7c87aed18ab73ada19296aa33aeec29c6fe0a73d Mon Sep 17 00:00:00 2001 From: Justin Stevens Date: Wed, 9 Oct 2019 21:47:39 -0400 Subject: [PATCH 06/19] add wavelength plot to truth_dirc --- src/plugins/Analysis/truth_dirc/DEventProcessor_truth_dirc.cc | 4 ++++ src/plugins/Analysis/truth_dirc/DEventProcessor_truth_dirc.h | 1 + 2 files changed, 5 insertions(+) diff --git a/src/plugins/Analysis/truth_dirc/DEventProcessor_truth_dirc.cc b/src/plugins/Analysis/truth_dirc/DEventProcessor_truth_dirc.cc index ca847c7e0..4bf50869d 100644 --- a/src/plugins/Analysis/truth_dirc/DEventProcessor_truth_dirc.cc +++ b/src/plugins/Analysis/truth_dirc/DEventProcessor_truth_dirc.cc @@ -34,6 +34,8 @@ jerror_t DEventProcessor_truth_dirc::init(void) { TDirectory *dircDir = gDirectory->mkdir("DIRC_truth"); dircDir->cd(); + hTruthWavelength = new TH1F("hTruthWavelength", "; Wavelength (nm)", 300, 200, 800); + int nChannels = 108*64; hTruthPixelHitTime = new TH2F("hTruthPixelHitTime", "; Pixel Channel # ; #Delta t (ns)", nChannels, 0, nChannels, 200, -100, 100); @@ -104,6 +106,7 @@ jerror_t DEventProcessor_truth_dirc::evnt(JEventLoop *loop, uint64_t eventnumber double z = dircPmtHits[h]->z; double t = dircPmtHits[h]->t; double t_fixed = dircPmtHits[h]->t_fixed; + double lambda = 1240./dircPmtHits[h]->E/1e9; // get PMT labels int pmt_column = dDIRCGeometry->GetPmtColumn(ch); @@ -114,6 +117,7 @@ jerror_t DEventProcessor_truth_dirc::evnt(JEventLoop *loop, uint64_t eventnumber int pixel_col = dDIRCGeometry->GetPixelColumn(ch); japp->RootWriteLock(); //ACQUIRE ROOT LOCK + hTruthWavelength->Fill(lambda); hTruthPixelHitTime->Fill(ch, t-t_fixed); if(x < 0.) { hTruthPmtHitZY_South->Fill(z, y); diff --git a/src/plugins/Analysis/truth_dirc/DEventProcessor_truth_dirc.h b/src/plugins/Analysis/truth_dirc/DEventProcessor_truth_dirc.h index af8f8b0ff..6d941ff04 100644 --- a/src/plugins/Analysis/truth_dirc/DEventProcessor_truth_dirc.h +++ b/src/plugins/Analysis/truth_dirc/DEventProcessor_truth_dirc.h @@ -56,6 +56,7 @@ class DEventProcessor_truth_dirc: public JEventProcessor { const DDIRCGeometry* dDIRCGeometry; + TH1F *hTruthWavelength; TH1F *hTruthBarHitBar; TH2F *hTruthBarHitXY; TH2F *hTruthPmtHitZY_North, *hTruthPmtHitZY_South; From 3f7c7e697cf64a66e5d1d4f6da102f0ff8be0e82 Mon Sep 17 00:00:00 2001 From: Justin Stevens Date: Wed, 9 Oct 2019 21:48:29 -0400 Subject: [PATCH 07/19] add option for timewalk plots (not include by default) --- .../DIRC_online/JEventProcessor_DIRC_online.cc | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/plugins/monitoring/DIRC_online/JEventProcessor_DIRC_online.cc b/src/plugins/monitoring/DIRC_online/JEventProcessor_DIRC_online.cc index 0b1cbdd15..40a73824a 100644 --- a/src/plugins/monitoring/DIRC_online/JEventProcessor_DIRC_online.cc +++ b/src/plugins/monitoring/DIRC_online/JEventProcessor_DIRC_online.cc @@ -100,6 +100,9 @@ JEventProcessor_DIRC_online::~JEventProcessor_DIRC_online() { jerror_t JEventProcessor_DIRC_online::init(void) { + FillTimewalk = false; + gPARMS->SetDefaultParameter("DIRC:FILLTIMEWALK", FillTimewalk, "Fill timewalk histograms, default = false"); + // create root folder for psc and cd to it, store main dir TDirectory *mainDir = gDirectory; TDirectory *dircDir = gDirectory->mkdir("DIRC_online"); @@ -135,7 +138,7 @@ jerror_t JEventProcessor_DIRC_online::init(void) { hHit_pixelOccupancy[i][j] = new TH2I("Hit_PixelOccupancy"+strN,"DIRCPmtHit pixel occupancy "+strT+"; pixel rows; pixel columns",Npixelcolumns,-0.5,-0.5+Npixelcolumns,Npixelrows,-0.5,-0.5+Npixelrows); hHit_TimeOverThreshold[i][j] = new TH1I("Hit_TimeOverThreshold"+strN,"DIRCPmtHit time-over-threshold "+strT+"; time-over-threshold (ns); hits",100,0.0,100.); hHit_TimeOverThresholdVsChannel[i][j] = new TH2I("Hit_TimeOverThresholdVsChannel"+strN,"DIRCPmtHit time-over-threshold vs channel "+strT+"; channel; time-over-threshold [ns]",Nchannels,-0.5,-0.5+Nchannels,100,0.0,100.); - hHit_tdcTime[i][j] = new TH1I("Hit_Time"+strN,"DIRCPmtHit time "+strT+";time [ns]; hits",500,0.0,500.0); + hHit_tdcTime[i][j] = new TH1I("Hit_Time"+strN,"DIRCPmtHit time "+strT+";time [ns]; hits",6500,-500.0,6000.0); hHit_tdcTimeVsEvent[i][j] = new TH2I("Hit_TimeVsEvent"+strN,"DIRCPmtHit time "+strT+"; event #; time [ns]; hits",1000,0,100e6,500,0.0,500.0); hHit_tdcTimeVsChannel[i][j] = new TH2I("Hit_TimeVsChannel"+strN,"DIRCPmtHit time vs. channel "+strT+"; channel;time [ns]",Nchannels,-0.5,-0.5+Nchannels,500,0.0,500.0); } @@ -147,7 +150,7 @@ jerror_t JEventProcessor_DIRC_online::init(void) { hHit_TimeEventMeanVsLEDRef[i] = new TH2I("Hit_TimeEventMeanVsLEDRef","LED Time Event Mean DIRCPmtHit time vs. LED Reference time; LED reference time [ns] ; LED pixel event mean time [ns]", 100, 100, 150, 400, -10, 30); hHit_TimeDiffEventMeanLEDRefVsTimestamp[i] = new TH2I("Hit_TimeDiffeventMeanLEDRefVsTimestamp","LED Time Event Mean DIRCPmtHit time - LED reference time vs. event timestamp; event timestamp [ns?] ; time difference [ns]", 1000, 0, 1e10, 400, 100, 140); - if(i==1) { + if(FillTimewalk) { gDirectory->mkdir("Timewalk")->cd(); for(int j=0; jt - locRefTime; if(ledFiber[1]) locDeltaT -= 10.; if(ledFiber[2]) locDeltaT -= 20.; - hHit_Timewalk[box][channel]->Fill(hit->tot, locDeltaT); + if(FillTimewalk) hHit_Timewalk[box][channel]->Fill(hit->tot, locDeltaT); } } From 3c82c0ddb3f685281b523c6b42d53c3b7603e60a Mon Sep 17 00:00:00 2001 From: Justin Stevens Date: Wed, 9 Oct 2019 21:48:51 -0400 Subject: [PATCH 08/19] add missing header --- src/plugins/monitoring/DIRC_online/JEventProcessor_DIRC_online.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/plugins/monitoring/DIRC_online/JEventProcessor_DIRC_online.h b/src/plugins/monitoring/DIRC_online/JEventProcessor_DIRC_online.h index 8c20758fb..ec37f2a7c 100644 --- a/src/plugins/monitoring/DIRC_online/JEventProcessor_DIRC_online.h +++ b/src/plugins/monitoring/DIRC_online/JEventProcessor_DIRC_online.h @@ -24,6 +24,7 @@ class JEventProcessor_DIRC_online:public jana::JEventProcessor{ jerror_t fini(void); ///< Called after last event of last event source has been processed. const DDIRCGeometry* dDIRCGeometry; + bool FillTimewalk; }; #endif // _JEventProcessor_DIRC_online_ From 85c704c2c7a8d47c768a6da431f7450f8464ebd1 Mon Sep 17 00:00:00 2001 From: Justin Stevens Date: Wed, 9 Oct 2019 22:15:39 -0400 Subject: [PATCH 09/19] add north box to DIRC_occupancy macros --- src/plugins/monitoring/occupancy_online/DIRC_occupancy.C | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/plugins/monitoring/occupancy_online/DIRC_occupancy.C b/src/plugins/monitoring/occupancy_online/DIRC_occupancy.C index f85c06305..73c470d5f 100644 --- a/src/plugins/monitoring/occupancy_online/DIRC_occupancy.C +++ b/src/plugins/monitoring/occupancy_online/DIRC_occupancy.C @@ -3,9 +3,7 @@ // which histograms to fetch for the macro. // // hnamepath: /occupancy/dirc_num_events -// hnamepath: /occupancy/dirc_tdc_pixel_N_occ // hnamepath: /occupancy/dirc_tdc_pixel_S_occ -// hnamepath: /occupancy/dirc_tdc_pixel_N_occ_led // hnamepath: /occupancy/dirc_tdc_pixel_S_occ_led // // e-mail: davidl@jlab.org @@ -17,14 +15,12 @@ // RootSpy saves the current directory and style before // calling the macro and restores it after so it is OK to // change them and not change them back. - TDirectory *savedir = gDirectory; + TDirectory *savedir = gDirectory; TDirectory *dir = (TDirectory*)gDirectory->FindObjectAny("occupancy"); if(dir) dir->cd(); - TH2I *dirc_tdc_pixel_N_occ_led = (TH2I*)gDirectory->FindObjectAny("dirc_tdc_pixel_N_occ_led"); TH2I *dirc_tdc_pixel_S_occ_led = (TH2I*)gDirectory->FindObjectAny("dirc_tdc_pixel_S_occ_led"); - TH2I *dirc_tdc_pixel_N_occ = (TH2I*)gDirectory->FindObjectAny("dirc_tdc_pixel_N_occ"); - TH2I *dirc_tdc_pixel_S_occ = (TH2I*)gDirectory->FindObjectAny("dirc_tdc_pixel_S_occ"); + TH2I *dirc_tdc_pixel_S_occ = (TH2I*)gDirectory->FindObjectAny("dirc_tdc_pixel_S_occ"); TH1I *dirc_num_events = (TH1I*)gDirectory->FindObjectAny("dirc_num_events"); double Nevents = 1.0; From bac4b520f818d4342e4418cdf649024d3bbab56d Mon Sep 17 00:00:00 2001 From: Justin Stevens Date: Wed, 9 Oct 2019 22:18:34 -0400 Subject: [PATCH 10/19] add north box to DIRC_occupancy macros --- .../occupancy_online/DIRC_North_occupancy.C | 69 +++++++++++++++++++ 1 file changed, 69 insertions(+) create mode 100644 src/plugins/monitoring/occupancy_online/DIRC_North_occupancy.C diff --git a/src/plugins/monitoring/occupancy_online/DIRC_North_occupancy.C b/src/plugins/monitoring/occupancy_online/DIRC_North_occupancy.C new file mode 100644 index 000000000..dfba30e58 --- /dev/null +++ b/src/plugins/monitoring/occupancy_online/DIRC_North_occupancy.C @@ -0,0 +1,69 @@ + +// The following are special comments used by RootSpy to know +// which histograms to fetch for the macro. +// +// hnamepath: /occupancy/dirc_num_events +// hnamepath: /occupancy/dirc_tdc_pixel_N_occ +// hnamepath: /occupancy/dirc_tdc_pixel_N_occ_led +// +// e-mail: davidl@jlab.org +// e-mail: marki@jlab.org +// e-mail: tbritton@jlab.org +// + +{ + // RootSpy saves the current directory and style before + // calling the macro and restores it after so it is OK to + // change them and not change them back. + TDirectory *savedir = gDirectory; + TDirectory *dir = (TDirectory*)gDirectory->FindObjectAny("occupancy"); + if(dir) dir->cd(); + + TH2I *dirc_tdc_pixel_N_occ_led = (TH2I*)gDirectory->FindObjectAny("dirc_tdc_pixel_N_occ_led"); + TH2I *dirc_tdc_pixel_N_occ = (TH2I*)gDirectory->FindObjectAny("dirc_tdc_pixel_N_occ"); + TH1I *dirc_num_events = (TH1I*)gDirectory->FindObjectAny("dirc_num_events"); + + double Nevents = 1.0; + if(dirc_num_events) Nevents = (double)dirc_num_events->GetBinContent(1); + + // Just for testing + if(gPad == NULL){ + TCanvas *c1 = new TCanvas("c1","DIRC North Occupancy",600,400); + c1->cd(0); + c1->Draw(); + c1->Update(); + } + if(!gPad) {savedir->cd(); return;} + + TCanvas *c1 = gPad->GetCanvas(); + c1->cd(0); + c1->Clear(); + + TPad *p1 = (TPad*)gDirectory->FindObjectAny("dirc_occ_pad1"); + if(!p1)p1 = new TPad("dirc_occ_pad1","p1",0.,0.5,1.,1.); + p1->Draw(); + p1->cd(); + gStyle->SetOptStat(0); + if(dirc_tdc_pixel_N_occ_led) dirc_tdc_pixel_N_occ_led->DrawCopy("colz"); + + c1->cd(0); + TPad *p2 = (TPad*)gDirectory->FindObjectAny("dirc_occ_pad2"); + if(!p2)p2 = new TPad("dirc_occ_pad2","p1",0.,0.,1.,0.5); + p2->Draw(); + p2->cd(); + gStyle->SetOptStat(0); + if(dirc_tdc_pixel_N_occ) dirc_tdc_pixel_N_occ->DrawCopy("colz"); + +#ifdef ROOTSPY_MACROS + // ------ The following is used by RSAI -------- + if( rs_GetFlag("Is_RSAI")==1 ){ + auto min_events = rs_GetFlag("MIN_EVENTS_RSAI"); + if( min_events < 1 ) min_events = 1E4; + if( Nevents >= min_events ) { + cout << "DIRC Flagging AI check after " << Nevents << " events (>=" << min_events << ")" << endl; + rs_SavePad("DIRC_North_occupancy", 0); + rs_ResetAllMacroHistos("//DIRC_North_occupancy"); + } + } +#endif +} From 307c72cbf79c88035983d0bebbbe71efc2deaf9d Mon Sep 17 00:00:00 2001 From: Justin Stevens Date: Tue, 22 Oct 2019 11:35:05 -0400 Subject: [PATCH 11/19] fixes for North OB in online monitoring --- src/libraries/DIRC/DDIRCPmtHit_factory.cc | 3 +- .../monitoring/DIRC_online/DIRC_North_hit.C | 73 +++++++++++++++++++ .../monitoring/DIRC_online/DIRC_digihit.C | 21 ------ src/plugins/monitoring/DIRC_online/DIRC_hit.C | 70 +++--------------- .../JEventProcessor_DIRC_online.cc | 4 +- 5 files changed, 89 insertions(+), 82 deletions(-) create mode 100644 src/plugins/monitoring/DIRC_online/DIRC_North_hit.C diff --git a/src/libraries/DIRC/DDIRCPmtHit_factory.cc b/src/libraries/DIRC/DDIRCPmtHit_factory.cc index 0256bb19c..b3c7242ee 100644 --- a/src/libraries/DIRC/DDIRCPmtHit_factory.cc +++ b/src/libraries/DIRC/DDIRCPmtHit_factory.cc @@ -125,6 +125,7 @@ jerror_t DDIRCPmtHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber) if(digihit_lead->channel != digihit_trail->channel) continue; int channel = digihit_lead->channel; int box = (channel < DIRC_MAX_CHANNELS) ? 1 : 0; // North=0 and South=1 + if(box == 0) channel -= DIRC_MAX_CHANNELS; // box-local channel to index CCDB tables // get time-over-threshold timeOverThreshold = (double)digihit_trail->time - (double)digihit_lead->time; @@ -152,7 +153,7 @@ jerror_t DDIRCPmtHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber) double slope = 0.3; double timeOverThresholdPeak = 50; if(applyTimeOffset) { - hit->t = hit->t - time_offsets[box][channel] + t_base[box]; + hit->t = hit->t + t_base[box] - time_offsets[box][channel]; } if(applyTimewalk) { hit->t += slope*(timeOverThreshold - timeOverThresholdPeak); diff --git a/src/plugins/monitoring/DIRC_online/DIRC_North_hit.C b/src/plugins/monitoring/DIRC_online/DIRC_North_hit.C new file mode 100644 index 000000000..dedb0c3f1 --- /dev/null +++ b/src/plugins/monitoring/DIRC_online/DIRC_North_hit.C @@ -0,0 +1,73 @@ +// The following are special comments used by RootSpy to know +// which histograms to fetch for the macro. +// + +// hnamepath: /DIRC_online/Hit/NorthUpperBox/Hit_TimeOverThresholdVsChannel_LED +// hnamepath: /DIRC_online/Hit/NorthUpperBox/Hit_TimeVsChannel_LED +// hnamepath: /DIRC_online/Hit/NorthUpperBox/Hit_Time_NonLED +// hnamepath: /DIRC_online/Hit/NorthUpperBox/Hit_Time_LED +// hnamepath: /DIRC_online/Hit/NorthUpperBox/Hit_PixelOccupancy_LED +// hnamepath: /DIRC_online/Hit/NorthUpperBox/Hit_PixelOccupancy_NonLED + + +{ + TDirectory *dir = (TDirectory*)gDirectory->FindObjectAny("DIRC_online"); + if(dir) dir->cd(); + + TH2I* hTimeVsChannel_LED = (TH2I*)gDirectory->Get("Hit/NorthUpperBox/Hit_TimeVsChannel_LED"); + TH1I* hTS = (TH1I*)gDirectory->Get("Hit/NorthUpperBox/Hit_Time_NonLED"); + TH1I* hTS_LED = (TH1I*)gDirectory->Get("Hit/NorthUpperBox/Hit_Time_LED"); + + TH2I* hOcc_LED = (TH2I*)gDirectory->Get("Hit/NorthUpperBox/Hit_PixelOccupancy_LED"); + TH2I* hOcc = (TH2I*)gDirectory->Get("Hit/NorthUpperBox/Hit_PixelOccupancy_NonLED"); + + if(gPad == NULL){ + TCanvas *c1 = new TCanvas("c1","DIRC Hit Monitor",150,10,990,660); + c1->cd(0); + c1->Draw(); + c1->Update(); + } + + if(!gPad) return; + TCanvas* c1 = gPad->GetCanvas(); + c1->Divide(2,2); + + double tsize = 0.05; + gStyle->SetOptStat("emr"); + + if(hOcc_LED){ + c1->cd(1); + hOcc_LED->SetTitleSize(tsize,"xy"); + hOcc_LED->Draw("colz"); + } + + if(hTimeVsChannel_LED){ + c1->cd(2); + hTimeVsChannel_LED->SetTitleSize(tsize,"xy"); + hTimeVsChannel_LED->Draw("colz"); + } + + if(hOcc){ + c1->cd(3); + hOcc->SetTitleSize(tsize,"xy"); + hOcc->Draw("colz"); + } + + + if(hTS && hTS_LED){ + hTS->SetLineColor(kBlack); + hTS_LED->SetLineColor(kBlue); + c1->cd(4); + hTS_LED->SetTitleSize(tsize,"xy"); + hTS_LED->Draw(); + hTS->Scale(hTS_LED->GetMaximum()/hTS->GetMaximum()); + hTS->Draw("h same"); + + TLegend *leg = new TLegend(0.6, 0.6, 0.85, 0.8); + leg->AddEntry(hTS_LED,"LED trigger","l"); + leg->AddEntry(hTS,"Non-LED triggers","l"); + leg->Draw("same"); + } + + +} diff --git a/src/plugins/monitoring/DIRC_online/DIRC_digihit.C b/src/plugins/monitoring/DIRC_online/DIRC_digihit.C index 0019bef2d..88ab5d444 100644 --- a/src/plugins/monitoring/DIRC_online/DIRC_digihit.C +++ b/src/plugins/monitoring/DIRC_online/DIRC_digihit.C @@ -13,8 +13,6 @@ TDirectory *dir = (TDirectory*)gDirectory->FindObjectAny("DIRC_online"); if(dir) dir->cd(); - //TH2I* hOcc = (TH2I*)gDirectory->Get("DigiHit/DigiHit_NHitsVsBox"); - //TH2I* hTN = (TH2I*)gDirectory->Get("DigiHit/NorthUpperBox/TDCDigiHit_TimeVsChannel_NorthUpperBox"); TH1I* hDigiHit_Nhits_LED = (TH1I*)gDirectory->Get("DigiHit/DigiHit_NHits_LED"); TH1I* hDigiHit_Nhits = (TH1I*)gDirectory->Get("DigiHit/DigiHit_NHits_NonLED"); TH1I* hDigiHit_Time_LED = (TH1I*)gDirectory->Get("DigiHit/SouthLowerBox/TDCDigiHit_Time_LED"); @@ -36,25 +34,6 @@ double tsize = 0.05; gStyle->SetOptStat("emr"); -/* - if(hOcc){ - hOcc->SetFillColor(kBlue); - c1->cd(1); - hOcc->SetLabelSize(0.08,"x"); - hOcc->SetTitleSize(tsize,"xy"); - hOcc->GetXaxis()->SetBinLabel(1, "North"); - hOcc->GetXaxis()->SetBinLabel(2, "South"); - hOcc->Draw("colz"); - } - - if(hTN){ - hTN->SetFillColor(kBlue); - c1->cd(3); - hTN->SetTitleSize(tsize,"xy"); - hTN->Draw("colz"); - } -*/ - TLegend *leg = new TLegend(0.6, 0.6, 0.85, 0.8); leg->AddEntry(hTS_LED,"LED trigger","l"); leg->AddEntry(hTS,"Non-LED triggers","l"); diff --git a/src/plugins/monitoring/DIRC_online/DIRC_hit.C b/src/plugins/monitoring/DIRC_online/DIRC_hit.C index e775a89fc..4c2656c6c 100644 --- a/src/plugins/monitoring/DIRC_online/DIRC_hit.C +++ b/src/plugins/monitoring/DIRC_online/DIRC_hit.C @@ -2,7 +2,6 @@ // which histograms to fetch for the macro. // -// hnamepath: /DIRC_online/Hit/Hit_NHits_LED // hnamepath: /DIRC_online/Hit/SouthLowerBox/Hit_TimeOverThresholdVsChannel_LED // hnamepath: /DIRC_online/Hit/SouthLowerBox/Hit_TimeVsChannel_LED // hnamepath: /DIRC_online/Hit/SouthLowerBox/Hit_Time_NonLED @@ -15,17 +14,12 @@ TDirectory *dir = (TDirectory*)gDirectory->FindObjectAny("DIRC_online"); if(dir) dir->cd(); - TH2I* hOcc = (TH2I*)gDirectory->Get("Hit/Hit_NHits_LED"); - //TH2I* hTotN = (TH2I*)gDirectory->Get("Hit/NorthUpperBox/Hit_TimeOverThresholdVsChannel_NorthUpperBox"); - TH2I* hTotS = (TH2I*)gDirectory->Get("Hit/SouthLowerBox/Hit_TimeOverThresholdVsChannel_LED"); - //TH2I* hTN = (TH2I*)gDirectory->Get("Hit/NorthUpperBox/Hit_TimeVsChannel_NorthUpperBox"); - TH2I* hTimeVsChannelS_LED = (TH2I*)gDirectory->Get("Hit/SouthLowerBox/Hit_TimeVsChannel_LED"); + TH2I* hTimeVsChannel_LED = (TH2I*)gDirectory->Get("Hit/SouthLowerBox/Hit_TimeVsChannel_LED"); TH1I* hTS = (TH1I*)gDirectory->Get("Hit/SouthLowerBox/Hit_Time_NonLED"); TH1I* hTS_LED = (TH1I*)gDirectory->Get("Hit/SouthLowerBox/Hit_Time_LED"); - //TH2I* hOccN = (TH2I*)gDirectory->Get("Hit/NorthUpperBox/Hit_PixelOccupancy_NorthUpperBox"); - TH2I* hOccS_LED = (TH2I*)gDirectory->Get("Hit/SouthLowerBox/Hit_PixelOccupancy_LED"); - TH2I* hOccS = (TH2I*)gDirectory->Get("Hit/SouthLowerBox/Hit_PixelOccupancy_NonLED"); + TH2I* hOcc_LED = (TH2I*)gDirectory->Get("Hit/SouthLowerBox/Hit_PixelOccupancy_LED"); + TH2I* hOcc = (TH2I*)gDirectory->Get("Hit/SouthLowerBox/Hit_PixelOccupancy_NonLED"); if(gPad == NULL){ TCanvas *c1 = new TCanvas("c1","DIRC Hit Monitor",150,10,990,660); @@ -41,62 +35,22 @@ double tsize = 0.05; gStyle->SetOptStat("emr"); -/* - if(hOcc){ - c1->cd(1); - hOcc->SetLabelSize(0.08,"x"); - hOcc->SetTitleSize(tsize,"xy"); - hOcc->GetXaxis()->SetBinLabel(1, "North"); - hOcc->GetXaxis()->SetBinLabel(2, "South"); - hOcc->Draw("colz"); - } - - if(hTotN){ - c1->cd(2); - hTotN->SetTitleSize(tsize,"xy"); - hTotN->Draw("colz"); - } - - if(hOccN && hOccS){ - c1->cd(4); - TPad *p1 = new TPad("p1","p1",0.0,0.5,1.0,1.0); - p1->Draw(); - p1->cd(); - gStyle->SetOptStat(0); - if(hOccN && hOccN->Integral() > 10) hOccN->Draw("colz"); - - c1->cd(4); - TPad *p2 = new TPad("p2","p2",0.0,0.0,1.0,0.5); - p2->Draw(); - p2->cd(); - gStyle->SetOptStat(0); - if(hOccS && hOccS->Integral() > 10) hOccS->Draw("colz"); - } - - if(hTN){ - hTN->SetFillColor(kBlue); - c1->cd(5); - hTN->SetTitleSize(tsize,"xy"); - hTN->Draw("colz"); - } -*/ - - if(hOccS_LED){ + if(hOcc_LED){ c1->cd(1); - hOccS_LED->SetTitleSize(tsize,"xy"); - hOccS_LED->Draw("colz"); + hOcc_LED->SetTitleSize(tsize,"xy"); + hOcc_LED->Draw("colz"); } - if(hTimeVsChannelS_LED){ + if(hTimeVsChannel_LED){ c1->cd(2); - hTimeVsChannelS_LED->SetTitleSize(tsize,"xy"); - hTimeVsChannelS_LED->Draw("colz"); + hTimeVsChannel_LED->SetTitleSize(tsize,"xy"); + hTimeVsChannel_LED->Draw("colz"); } - if(hOccS){ + if(hOcc){ c1->cd(3); - hOccS->SetTitleSize(tsize,"xy"); - hOccS->Draw("colz"); + hOcc->SetTitleSize(tsize,"xy"); + hOcc->Draw("colz"); } diff --git a/src/plugins/monitoring/DIRC_online/JEventProcessor_DIRC_online.cc b/src/plugins/monitoring/DIRC_online/JEventProcessor_DIRC_online.cc index fbf8375e0..32f2efa40 100644 --- a/src/plugins/monitoring/DIRC_online/JEventProcessor_DIRC_online.cc +++ b/src/plugins/monitoring/DIRC_online/JEventProcessor_DIRC_online.cc @@ -138,7 +138,7 @@ jerror_t JEventProcessor_DIRC_online::init(void) { hHit_pixelOccupancy[i][j] = new TH2I("Hit_PixelOccupancy"+strN,"DIRCPmtHit pixel occupancy "+strT+"; pixel rows; pixel columns",Npixelcolumns,-0.5,-0.5+Npixelcolumns,Npixelrows,-0.5,-0.5+Npixelrows); hHit_TimeOverThreshold[i][j] = new TH1I("Hit_TimeOverThreshold"+strN,"DIRCPmtHit time-over-threshold "+strT+"; time-over-threshold (ns); hits",100,0.0,100.); hHit_TimeOverThresholdVsChannel[i][j] = new TH2I("Hit_TimeOverThresholdVsChannel"+strN,"DIRCPmtHit time-over-threshold vs channel "+strT+"; channel; time-over-threshold [ns]",Nchannels,-0.5,-0.5+Nchannels,100,0.0,100.); - hHit_tdcTime[i][j] = new TH1I("Hit_Time"+strN,"DIRCPmtHit time "+strT+";time [ns]; hits",6500,-500.0,6000.0); + hHit_tdcTime[i][j] = new TH1I("Hit_Time"+strN,"DIRCPmtHit time "+strT+";time [ns]; hits",500,0.0,500.0); hHit_tdcTimeVsEvent[i][j] = new TH2I("Hit_TimeVsEvent"+strN,"DIRCPmtHit time "+strT+"; event #; time [ns]; hits",1000,0,100e6,500,0.0,500.0); hHit_tdcTimeVsChannel[i][j] = new TH2I("Hit_TimeVsChannel"+strN,"DIRCPmtHit time vs. channel "+strT+"; channel;time [ns]",Nchannels,-0.5,-0.5+Nchannels,500,0.0,500.0); } @@ -338,7 +338,7 @@ jerror_t JEventProcessor_DIRC_online::evnt(JEventLoop *eventLoop, uint64_t event // Loop over calibrated hits to get mean for reference time double locRefTime = 0; int locNHits = 0; - double locFirstFiberTime = 128.1; // 205.5; used for no time offset + double locFirstFiberTime = 125.5; // 205.5; used for no time offset for (const auto& hit : hits) { int channel = (hit->ch < Nchannels) ? hit->ch : (hit->ch - Nchannels); int pmtrow = locDIRCGeometry->GetPmtRow(channel); From 568b0893895936f83c9e2a22ca729002025860bf Mon Sep 17 00:00:00 2001 From: Justin Stevens Date: Tue, 22 Oct 2019 12:16:44 -0400 Subject: [PATCH 12/19] remove thetaC corrections from this branch --- src/libraries/DIRC/DDIRCLut.cc | 11 ----------- src/libraries/DIRC/DDIRCLut.h | 5 ----- 2 files changed, 16 deletions(-) diff --git a/src/libraries/DIRC/DDIRCLut.cc b/src/libraries/DIRC/DDIRCLut.cc index 8b5a71ada..8dc8147dc 100644 --- a/src/libraries/DIRC/DDIRCLut.cc +++ b/src/libraries/DIRC/DDIRCLut.cc @@ -50,9 +50,6 @@ DDIRCLut::DDIRCLut() dCriticalAngle = asin(1.00028/1.47125); // n_quarzt = 1.47125; //(1.47125 <==> 390nm) dIndex = 1.473; - vector new_thetac(108); - dThetaC_offsets.push_back(new_thetac); dThetaC_offsets.push_back(new_thetac); - if(DIRC_DEBUG_HISTS) CreateDebugHistograms(); } @@ -67,12 +64,6 @@ bool DDIRCLut::brun(JEventLoop *loop) { loop->Get(locDIRCGeometry); dDIRCGeometry = locDIRCGeometry[0]; - // load constant tables - if (loop->GetCalib("/DIRC/North/thetac_offsets", dThetaC_offsets[0])) - jout << "Error loading /DIRC/North/thetac_offsets !" << endl; - if (loop->GetCalib("/DIRC/South/thetac_offsets", dThetaC_offsets[1])) - jout << "Error loading /DIRC/South/thetac_offsets !" << endl; - return true; } @@ -243,7 +234,6 @@ vector> DDIRCLut::CalcPhoton(const DDIRCPmtHit *locDIRCHit, return locDIRCPhotons; int pmt = channel/64; - double thetac_offset = dThetaC_offsets[box][pmt]; // use hit time to determine if reflected or not double hitTime = locDIRCHit->t - locFlightTime; @@ -303,7 +293,6 @@ vector> DDIRCLut::CalcPhoton(const DDIRCPmtHit *locDIRCHit, luttheta = dir.Angle(TVector3(-1,0,0)); if(luttheta > TMath::PiOver2()) luttheta = TMath::Pi()-luttheta; tangle = momInBar.Angle(dir); - tangle -= thetac_offset; //correction double bartime = lenz/cos(luttheta)/DIRC_LIGHT_V; double totalTime = bartime+evtime; diff --git a/src/libraries/DIRC/DDIRCLut.h b/src/libraries/DIRC/DDIRCLut.h index 5123614bb..c371c3f9c 100644 --- a/src/libraries/DIRC/DDIRCLut.h +++ b/src/libraries/DIRC/DDIRCLut.h @@ -25,9 +25,6 @@ using namespace jana; #include "TH2.h" #include "TF1.h" -// store constants so that they can be accessed by pixel number -typedef vector dirc_lut_constants_t; - class DDIRCLut: public JObject { public: @@ -37,8 +34,6 @@ class DDIRCLut: public JObject { DDIRCLut(); ~DDIRCLut(){}; - vector dThetaC_offsets; - bool brun(JEventLoop *loop); bool CreateDebugHistograms(); bool CalcLUT(TVector3 locProjPos, TVector3 locProjMom, const vector locDIRCHits, double locFlightTime, double locMass, shared_ptr& locDIRCMatchParams, const vector locDIRCBarHits, map, vector >& locDIRCTrackMatchParams) const; From 4ba12efdc94655db14fdf7ef51e395482cd35a6e Mon Sep 17 00:00:00 2001 From: "GlueX-Hall-D (Mark Ito)" Date: Tue, 22 Oct 2019 12:28:32 -0400 Subject: [PATCH 13/19] add truth information to dirc_tree --- .../dirc_tree/DEventProcessor_dirc_tree.cc | 30 ++++++++++++++++++- src/plugins/Analysis/dirc_tree/DrcEvent.h | 14 ++++++++- 2 files changed, 42 insertions(+), 2 deletions(-) diff --git a/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc b/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc index 1b5569504..87cdfcc2f 100644 --- a/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc +++ b/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc @@ -95,6 +95,10 @@ jerror_t DEventProcessor_dirc_tree::evnt(jana::JEventLoop* loop, uint64_t locEve loop->GetSingle(locDetectorMatches); DDetectorMatches locDetectorMatch = (DDetectorMatches)locDetectorMatches[0]; + // cheat and get truth info of track at bar + vector locDIRCBarHits; + loop->Get(locDIRCBarHits); + japp->RootWriteLock(); TClonesArray& cevt = *fcEvent; @@ -153,7 +157,7 @@ jerror_t DEventProcessor_dirc_tree::evnt(jana::JEventLoop* loop, uint64_t locEve DVector3 posInBar = locDIRCMatchParams->dExtrapolatedPos; DVector3 momInBar = locDIRCMatchParams->dExtrapolatedMom; double locExtrapolatedTime = locDIRCMatchParams->dExtrapolatedTime; - int locBar = dDIRCGeometry->GetBar(posInBar.Y()); + int locBar = locDIRCGeometry->GetBar(posInBar.Y()); fEvent = new DrcEvent(); fEvent->SetType(2); @@ -170,6 +174,30 @@ jerror_t DEventProcessor_dirc_tree::evnt(jana::JEventLoop* loop, uint64_t locEve fEvent->SetMissMass(locMissingP4.M2()); fEvent->SetChiSq(chisq); + if(locDIRCBarHits.size() > 0) { + + TVector3 bestMatchPos, bestMatchMom; + double bestFlightTime = 999.; + double bestMatchDist = 999.; + int bestBar = -1; + for(int i=0; i<(int)locDIRCBarHits.size(); i++) { + TVector3 locDIRCBarHitPos(locDIRCBarHits[i]->x, locDIRCBarHits[i]->y, locDIRCBarHits[i]->z); + TVector3 locDIRCBarHitMom(locDIRCBarHits[i]->px, locDIRCBarHits[i]->py, locDIRCBarHits[i]->pz); + if((posInBar - locDIRCBarHitPos).Mag() < bestMatchDist) { + bestMatchDist = (posInBar - locDIRCBarHitPos).Mag(); + bestMatchPos = locDIRCBarHitPos; + bestMatchMom = locDIRCBarHitMom; + bestFlightTime = locDIRCBarHits[i]->t; + bestBar = locDIRCBarHits[i]->bar; + } + } + + fEvent->SetMomentum_Truth(TVector3(bestMatchMom.X(), bestMatchMom.Y(), bestMatchMom.Z())); + fEvent->SetPosition_Truth(TVector3(bestMatchPos.X(), bestMatchPos.Y(), bestMatchPos.Z())); + fEvent->SetTime_Truth(bestFlightTime); + fEvent->SetId_Truth(bestBar); + } + DrcHit hit; for(const auto dhit : locDIRCPmtHits){ int ch=dhit->ch; diff --git a/src/plugins/Analysis/dirc_tree/DrcEvent.h b/src/plugins/Analysis/dirc_tree/DrcEvent.h index 5fd949c0a..13e0b69c3 100644 --- a/src/plugins/Analysis/dirc_tree/DrcEvent.h +++ b/src/plugins/Analysis/dirc_tree/DrcEvent.h @@ -40,6 +40,10 @@ class DrcEvent: public TObject { Int_t GetDcHits() const { return fDcHits; } TVector3 GetMomentum() const { return fMomentum; } TVector3 GetPosition() const { return fPosition; } + Int_t GetId_Truth() const { return fId_Truth; } + Double_t GetTime_Truth() const { return fTime_Truth; } + TVector3 GetMomentum_Truth() const { return fMomentum_Truth; } + TVector3 GetPosition_Truth() const { return fPosition_Truth; } Int_t GetHitSize() const { return fHitSize; } Double_t GetTest1() const { return fTest1; } Double_t GetTest2() const { return fTest2; } @@ -59,6 +63,10 @@ class DrcEvent: public TObject { void SetDcHits(Int_t val) { fDcHits = val; } void SetMomentum(TVector3 val){ fMomentum = val; } void SetPosition(TVector3 val){ fPosition = val; } + void SetId_Truth(Int_t val) { fId_Truth = val; } + void SetTime_Truth(Double_t val) { fTime_Truth = val; } + void SetMomentum_Truth(TVector3 val){ fMomentum_Truth = val; } + void SetPosition_Truth(TVector3 val){ fPosition_Truth = val; } void SetTest1(Double_t val) { fTest1 = val; } void SetTest2(Double_t val) { fTest2 = val; } @@ -74,12 +82,16 @@ class DrcEvent: public TObject { Double_t fChiSq; Double_t fTofTrackDist; Double_t fTofTrackDeltaT; - + Int_t fId_Truth; + Double_t fTime_Truth; + Int_t fHitSize; std::vector fHitArray; TVector3 fMomentum; TVector3 fPosition; + TVector3 fMomentum_Truth; + TVector3 fPosition_Truth; Double_t fTest1; Double_t fTest2; From edc56b6a5ad30a567d7dd25d8abb13ef42fcf849 Mon Sep 17 00:00:00 2001 From: "GlueX-Hall-D (Mark Ito)" Date: Tue, 22 Oct 2019 12:30:22 -0400 Subject: [PATCH 14/19] add TOF position residuals --- .../DCustomAction_dirc_reactions.cc | 18 +++++++++++++++++- .../DCustomAction_dirc_reactions.h | 2 ++ 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.cc b/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.cc index cfd7add08..cdf7caa4f 100644 --- a/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.cc +++ b/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.cc @@ -74,6 +74,8 @@ void DCustomAction_dirc_reactions::Initialize(JEventLoop* locEventLoop) CreateAndChangeTo_Directory(Form("Bar %d", locBar), Form("Bar %d", locBar)); double bar_y = dDIRCGeometry->GetBarY(locBar); + hDeltaThetaCVsChannelMap[locBar] = GetOrCreate_Histogram(Form("hDeltaThetaCVsChannel_%s_%d",locParticleName.data(),locBar), Form("cherenkov angle vs. channel; channel ID; %s Bar %d #Delta#theta_{C} [rad]", locParticleROOTName.data(),locBar), 6912, 0, 6912, 200,-0.2,0.2); + for(int locXbin=0; locXbin<40; locXbin++) { double xbin_min = -100.0 + locXbin*5.0; @@ -83,7 +85,8 @@ void DCustomAction_dirc_reactions::Initialize(JEventLoop* locEventLoop) TVector3 locBinVect((xbin_min+xbin_max)/2., bar_y, 586.5 - 65.); if(locBinVect.Theta()*180/TMath::Pi() > 12.) continue; - + + hDeltaTOF[locBar][locXbin] = GetOrCreate_Histogram(Form("hDeltaTOF_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; %s TOF #Delta x [cm]; %s TOF #Delta y [cm]; entries [#]", locBar,xbin_min,xbin_max,locParticleROOTName.data(),locParticleROOTName.data()), 40,-10,10,40,-10,10); hDiffMap[locBar][locXbin] = GetOrCreate_Histogram(Form("hDiff_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; %s t_{calc}-t_{measured} [ns]; entries [#]", locBar,xbin_min,xbin_max,locParticleROOTName.data()), 100,-20,20); hHitTimeMap[locBar][locXbin] = GetOrCreate_Histogram(Form("hHitTimeMap_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f]; %s t_{measured} [ns]; entries [#]", locBar,xbin_min,xbin_max,locParticleROOTName.data()), 100,0,150); hNphCMap[locBar][locXbin] = GetOrCreate_Histogram(Form("hNphC_%s_%d_%d",locParticleName.data(),locBar,locXbin), Form("Bar %d, xbin [%0.0f,%0.0f] # photons; %s # photons", locBar,xbin_min,xbin_max,locParticleROOTName.data()), 80, 0, 80); @@ -153,6 +156,16 @@ bool DCustomAction_dirc_reactions::Perform_Action(JEventLoop* locEventLoop, cons if(!foundTOF || locTOFHitMatchParams->dDeltaXToHit > 10.0 || locTOFHitMatchParams->dDeltaYToHit > 10.0) return true; + /* + // get FCAL for check on residual + shared_ptr locFCALShowerMatchParams; + bool foundFCAL = dParticleID->Get_BestFCALMatchParams(locTrackTimeBased, locDetectorMatches, locFCALShowerMatchParams); + const DFCALShower* locMatchedFCALShower; + if(foundFCAL) { + locMatchedFCALShower = locFCALShowerMatchParams->dFCALShower; + } + */ + Particle_t locPID = locTrackTimeBased->PID(); double locMass = ParticleMass(locPID); @@ -201,6 +214,7 @@ bool DCustomAction_dirc_reactions::Perform_Action(JEventLoop* locEventLoop, cons Lock_Action(); //ACQUIRE ROOT LOCK!! hExtrapolatedBarHitXY->Fill(posInBar.X(), posInBar.Y()); + if(DIRC_FILL_BAR_MAP) hDeltaTOF[locBar][locXbin]->Fill(locTOFHitMatchParams->dDeltaXToHit,locTOFHitMatchParams->dDeltaYToHit); Unlock_Action(); //RELEASE ROOT LOCK!! double locAngle = dDIRCLut->CalcAngle(locP, locMass); @@ -278,8 +292,10 @@ bool DCustomAction_dirc_reactions::Perform_Action(JEventLoop* locEventLoop, cons hThetaCVsP->Fill(momInBar.Mag(), locThetaC); hDeltaThetaCVsP->Fill(momInBar.Mag(), locThetaC-locExpectedThetaC); hDeltaThetaCVsChannel->Fill(locChannel, locThetaC-locExpectedThetaC); + if(DIRC_FILL_BAR_MAP) { hDeltaThetaCVsPMap[locBar][locXbin]->Fill(momInBar.Mag(), locThetaC-locExpectedThetaC); + hDeltaThetaCVsChannelMap[locBar]->Fill(locChannel, locThetaC-locExpectedThetaC); hThetaCVsPMap[locBar][locXbin]->Fill(momInBar.Mag(), locThetaC); } Unlock_Action(); //RELEASE ROOT LOCK!! diff --git a/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.h b/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.h index afc611515..0ad220a67 100644 --- a/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.h +++ b/src/plugins/Analysis/dirc_reactions/DCustomAction_dirc_reactions.h @@ -60,8 +60,10 @@ class DCustomAction_dirc_reactions : public DAnalysisAction TH1I *hDiff, *hNphC, *hThetaC, *hDeltaThetaC, *hLikelihood, *hLikelihoodDiff; TH2I *hThetaCVsP, *hDeltaThetaCVsP, *hLikelihoodDiffVsP, *hReactionLikelihoodDiffVsP; TH2I *hDeltaThetaCVsChannel; + TH2I *hDeltaThetaCVsChannelMap[48]; TH2I *hExtrapolatedBarHitXY, *hExtrapolatedBarHitXY_PreCut; + TH2I *hDeltaTOF[48][40]; TH1I *hDiffMap[48][40], *hHitTimeMap[48][40], *hHitTimeMapAll[48][40]; TH1I *hNphCMap[48][40], *hNphCMapSlot4[48][40], *hNphCMapSlot5[48][40]; TH3S *hPixelHitMap3D[48][40]; From df64191f31165583ce9e2d024f1f2ed1eb3e00d5 Mon Sep 17 00:00:00 2001 From: "GlueX-Hall-D (Mark Ito)" Date: Tue, 22 Oct 2019 12:31:28 -0400 Subject: [PATCH 15/19] add extrapolated X,Y position at DIRC to standard analysis TTree --- src/libraries/ANALYSIS/DEventWriterROOT.cc | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/libraries/ANALYSIS/DEventWriterROOT.cc b/src/libraries/ANALYSIS/DEventWriterROOT.cc index 9795bcbcd..9201befc6 100644 --- a/src/libraries/ANALYSIS/DEventWriterROOT.cc +++ b/src/libraries/ANALYSIS/DEventWriterROOT.cc @@ -672,6 +672,8 @@ void DEventWriterROOT::Create_Branches_ChargedHypotheses(DTreeBranchRegister& lo //DIRC: if(DIRC_OUTPUT) { locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "NumPhotons_DIRC"), locArraySizeString, dInitNumTrackArraySize); + locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "ExtrapolatedX_DIRC"), locArraySizeString, dInitNumTrackArraySize); + locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "ExtrapolatedY_DIRC"), locArraySizeString, dInitNumTrackArraySize); locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "ThetaC_DIRC"), locArraySizeString, dInitNumTrackArraySize); locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "Lele_DIRC"), locArraySizeString, dInitNumTrackArraySize); locBranchRegister.Register_FundamentalArray(Build_BranchName(locParticleBranchName, "Lpi_DIRC"), locArraySizeString, dInitNumTrackArraySize); @@ -1625,6 +1627,8 @@ void DEventWriterROOT::Fill_ChargedHypo(DTreeFillData* locTreeFillData, unsigned // DIRC if(DIRC_OUTPUT) { int locDIRCNumPhotons = 0; + double locDIRCExtrapolatedX = 999; + double locDIRCExtrapolatedY = 999; double locDIRCThetaC = 999.; double locDIRCLele = 999.; double locDIRCLpi = 999.; @@ -1632,6 +1636,8 @@ void DEventWriterROOT::Fill_ChargedHypo(DTreeFillData* locTreeFillData, unsigned double locDIRCLp = 999.; auto locDIRCMatchParams = locChargedTrackHypothesis->Get_DIRCMatchParams(); if(locDIRCMatchParams != NULL) { + locDIRCExtrapolatedX = locDIRCMatchParams->dExtrapolatedPos.X(); + locDIRCExtrapolatedY = locDIRCMatchParams->dExtrapolatedPos.Y(); locDIRCNumPhotons = locDIRCMatchParams->dNPhotons; locDIRCThetaC = locDIRCMatchParams->dThetaC; locDIRCLele = locDIRCMatchParams->dLikelihoodElectron; @@ -1640,6 +1646,8 @@ void DEventWriterROOT::Fill_ChargedHypo(DTreeFillData* locTreeFillData, unsigned locDIRCLp = locDIRCMatchParams->dLikelihoodProton; } locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "NumPhotons_DIRC"), locDIRCNumPhotons, locArrayIndex); + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "ExtrapolatedX_DIRC"), locDIRCExtrapolatedX, locArrayIndex); + locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "ExtrapolatedY_DIRC"), locDIRCExtrapolatedY, locArrayIndex); locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "ThetaC_DIRC"), locDIRCThetaC, locArrayIndex); locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "Lele_DIRC"), locDIRCLele, locArrayIndex); locTreeFillData->Fill_Array(Build_BranchName(locParticleBranchName, "Lpi_DIRC"), locDIRCLpi, locArrayIndex); From 9a363c88bdc179b23a1cfef1241c117cc1ce7964 Mon Sep 17 00:00:00 2001 From: "GlueX-Hall-D (Mark Ito)" Date: Tue, 22 Oct 2019 12:32:14 -0400 Subject: [PATCH 16/19] add ability to skip DIRC hits with flag --- src/libraries/DIRC/DDIRCPmtHit_factory.cc | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/libraries/DIRC/DDIRCPmtHit_factory.cc b/src/libraries/DIRC/DDIRCPmtHit_factory.cc index 0256bb19c..b9b04fa81 100644 --- a/src/libraries/DIRC/DDIRCPmtHit_factory.cc +++ b/src/libraries/DIRC/DDIRCPmtHit_factory.cc @@ -18,6 +18,9 @@ using namespace jana; //------------------ jerror_t DDIRCPmtHit_factory::init(void) { + DIRC_SKIP = false; + gPARMS->SetDefaultParameter("DIRC:SKIP",DIRC_SKIP); + // initialize calibration tables vector new_t0s(DIRC_MAX_CHANNELS); vector new_status(DIRC_MAX_CHANNELS); @@ -84,6 +87,9 @@ jerror_t DDIRCPmtHit_factory::evnt(JEventLoop *loop, uint64_t eventnumber) /// data in HDDM format. The HDDM event source will copy /// the precalibrated values directly into the _data vector. + if(DIRC_SKIP) + return NOERROR; + // check that SSP board timestamps match for all modules vector timestamps; loop->Get(timestamps); From 6feab403ace6afb09c162081b0f8e2edc8901a9d Mon Sep 17 00:00:00 2001 From: "GlueX-Hall-D (Mark Ito)" Date: Tue, 22 Oct 2019 12:34:11 -0400 Subject: [PATCH 17/19] add ability to skip DIRC hits with flag --- src/libraries/DIRC/DDIRCPmtHit_factory.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/libraries/DIRC/DDIRCPmtHit_factory.h b/src/libraries/DIRC/DDIRCPmtHit_factory.h index 433130632..329b54e16 100644 --- a/src/libraries/DIRC/DDIRCPmtHit_factory.h +++ b/src/libraries/DIRC/DDIRCPmtHit_factory.h @@ -44,6 +44,7 @@ class DDIRCPmtHit_factory:public jana::JFactory{ double t_base[2]; bool applyTimeOffset, applyTimewalk; + bool DIRC_SKIP; }; #endif // _DDIRCPmtHit_factory_ From 9563536c5b9029f23b47a46e70f4d3c16c34cedd Mon Sep 17 00:00:00 2001 From: "GlueX-Hall-D (Mark Ito)" Date: Tue, 22 Oct 2019 12:38:41 -0400 Subject: [PATCH 18/19] add possibility for DIRC bar rotation --- src/libraries/DIRC/DDIRCLut.cc | 44 ++++++++++++++++++++++++---------- src/libraries/DIRC/DDIRCLut.h | 8 +++---- 2 files changed, 34 insertions(+), 18 deletions(-) diff --git a/src/libraries/DIRC/DDIRCLut.cc b/src/libraries/DIRC/DDIRCLut.cc index 8b5a71ada..6dab62847 100644 --- a/src/libraries/DIRC/DDIRCLut.cc +++ b/src/libraries/DIRC/DDIRCLut.cc @@ -39,6 +39,10 @@ DDIRCLut::DDIRCLut() DIRC_SIGMA_THETAC = 0.01; gPARMS->SetDefaultParameter("DIRC:SIGMA_THETAC",DIRC_SIGMA_THETAC); + // Rotate tracks angle based on bar box survey data + DIRC_ROTATE_TRACK = false; + gPARMS->SetDefaultParameter("DIRC:ROTATE_TRACK",DIRC_ROTATE_TRACK); + // set PID for different passes in debuging histograms dFinalStatePIDs.push_back(Positron); dFinalStatePIDs.push_back(PiPlus); @@ -50,9 +54,6 @@ DDIRCLut::DDIRCLut() dCriticalAngle = asin(1.00028/1.47125); // n_quarzt = 1.47125; //(1.47125 <==> 390nm) dIndex = 1.473; - vector new_thetac(108); - dThetaC_offsets.push_back(new_thetac); dThetaC_offsets.push_back(new_thetac); - if(DIRC_DEBUG_HISTS) CreateDebugHistograms(); } @@ -67,11 +68,17 @@ bool DDIRCLut::brun(JEventLoop *loop) { loop->Get(locDIRCGeometry); dDIRCGeometry = locDIRCGeometry[0]; - // load constant tables - if (loop->GetCalib("/DIRC/North/thetac_offsets", dThetaC_offsets[0])) - jout << "Error loading /DIRC/North/thetac_offsets !" << endl; - if (loop->GetCalib("/DIRC/South/thetac_offsets", dThetaC_offsets[1])) - jout << "Error loading /DIRC/South/thetac_offsets !" << endl; + // rotation angles for bar boxes (reverse sign from survey) + for(int i=0; i<12; i++) { + dRotationX[i] = -0.2134*TMath::DegToRad(); // +pitch // -0.00219; + dRotationY[i] = 0.0963*TMath::DegToRad(); // -yaw // 0.00181 + dRotationZ[i] = -0.0355*TMath::DegToRad(); // -roll // -0.00062; + } + for(int i=12; i<24; i++) { + dRotationX[i] = -0.0791*TMath::DegToRad(); // +pitch // -0.00138; + dRotationY[i] = 0.1003*TMath::DegToRad(); // -yaw // 0.00175; + dRotationZ[i] = -0.0381*TMath::DegToRad(); // -roll // -0.00066; + } return true; } @@ -243,8 +250,14 @@ vector> DDIRCLut::CalcPhoton(const DDIRCPmtHit *locDIRCHit, return locDIRCPhotons; int pmt = channel/64; - double thetac_offset = dThetaC_offsets[box][pmt]; - + double rotX = 0, rotY = 0, rotZ=0; + int xbin = (int)(posInBar.X()/5.0) + 19; + if(bar>=0 && bar<=23 && xbin>=0 && xbin<=39) { + rotX = dRotationX[bar]; + rotY = dRotationY[bar]; + rotZ = dRotationZ[bar]; + } + // use hit time to determine if reflected or not double hitTime = locDIRCHit->t - locFlightTime; @@ -300,11 +313,16 @@ vector> DDIRCLut::CalcPhoton(const DDIRCPmtHit *locDIRCHit, if(r) dir.SetXYZ( -dir.X(), dir.Y(), dir.Z()); if(dir.Angle(fnY1) < dCriticalAngle || dir.Angle(fnZ1) < dCriticalAngle) continue; + TVector3 trackMom = momInBar; + if(DIRC_ROTATE_TRACK) { // rotate tracks to bar plane from survey data + trackMom.RotateX(rotX); + trackMom.RotateY(rotY); + trackMom.RotateZ(rotZ); + } + tangle = trackMom.Angle(dir); + luttheta = dir.Angle(TVector3(-1,0,0)); if(luttheta > TMath::PiOver2()) luttheta = TMath::Pi()-luttheta; - tangle = momInBar.Angle(dir); - tangle -= thetac_offset; //correction - double bartime = lenz/cos(luttheta)/DIRC_LIGHT_V; double totalTime = bartime+evtime; diff --git a/src/libraries/DIRC/DDIRCLut.h b/src/libraries/DIRC/DDIRCLut.h index 5123614bb..f8a2ddeca 100644 --- a/src/libraries/DIRC/DDIRCLut.h +++ b/src/libraries/DIRC/DDIRCLut.h @@ -25,9 +25,6 @@ using namespace jana; #include "TH2.h" #include "TF1.h" -// store constants so that they can be accessed by pixel number -typedef vector dirc_lut_constants_t; - class DDIRCLut: public JObject { public: @@ -37,8 +34,6 @@ class DDIRCLut: public JObject { DDIRCLut(); ~DDIRCLut(){}; - vector dThetaC_offsets; - bool brun(JEventLoop *loop); bool CreateDebugHistograms(); bool CalcLUT(TVector3 locProjPos, TVector3 locProjMom, const vector locDIRCHits, double locFlightTime, double locMass, shared_ptr& locDIRCMatchParams, const vector locDIRCBarHits, map, vector >& locDIRCTrackMatchParams) const; @@ -56,12 +51,15 @@ class DDIRCLut: public JObject { bool DIRC_DEBUG_HISTS; bool DIRC_TRUTH_BARHIT; bool DIRC_TRUTH_PIXELTIME; + bool DIRC_ROTATE_TRACK; double DIRC_CUT_TDIFFD; double DIRC_CUT_TDIFFR; double DIRC_SIGMA_THETAC; double DIRC_LIGHT_V; + double dRotationX[48], dRotationY[48], dRotationZ[48]; + int dMaxChannels; double dCriticalAngle, dIndex; From 8a20e897587fa292830eaebbc6477cc5cd08f1b9 Mon Sep 17 00:00:00 2001 From: Justin Stevens Date: Tue, 22 Oct 2019 12:57:12 -0400 Subject: [PATCH 19/19] get correct variable for DIRC geometry --- src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc b/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc index 8ce03bb7e..34a92de1f 100644 --- a/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc +++ b/src/plugins/Analysis/dirc_tree/DEventProcessor_dirc_tree.cc @@ -158,7 +158,7 @@ jerror_t DEventProcessor_dirc_tree::evnt(jana::JEventLoop* loop, uint64_t locEve DVector3 posInBar = locDIRCMatchParams->dExtrapolatedPos; DVector3 momInBar = locDIRCMatchParams->dExtrapolatedMom; double locExtrapolatedTime = locDIRCMatchParams->dExtrapolatedTime; - int locBar = locDIRCGeometry->GetBar(posInBar.Y()); + int locBar = locDIRCGeometryVec[0]->GetBar(posInBar.Y()); fEvent = new DrcEvent(); fEvent->SetType(2);