11 if(gPARMS->Exists(
"DIRC:TRUTH_BARHIT"))
19 locEventLoop->GetSingle(locParticleID);
26 vector<const DDIRCGeometry*> locDIRCGeometry;
27 locEventLoop->Get(locDIRCGeometry);
38 japp->RootWriteLock();
46 string locLikelihoodName =
"ln L(#pi) - ln L(K)";
48 locLikelihoodName =
"ln L(e) - ln L(#pi)";
50 locLikelihoodName =
"ln L(K) - ln L(p)";
52 hExtrapolatedBarHitXY_PreCut = GetOrCreate_Histogram<TH2I>(Form(
"hExtrapolatedBarHitXY_PreCut_%s",locParticleName.data()),
"; Bar Hit X (cm); Bar Hit Y (cm)", 200, -100, 100, 200, -100, 100);
53 hExtrapolatedBarHitXY = GetOrCreate_Histogram<TH2I>(Form(
"hExtrapolatedBarHitXY_%s",locParticleName.data()),
"; Bar Hit X (cm); Bar Hit Y (cm)", 200, -100, 100, 200, -100, 100);
55 hDiff = GetOrCreate_Histogram<TH1I>(Form(
"hDiff_%s",locParticleName.data()), Form(
"; %s t_{calc}-t_{measured} [ns]; entries [#]", locParticleROOTName.data()), 400,-20,20);
56 hNphC = GetOrCreate_Histogram<TH1I>(Form(
"hNphC_%s",locParticleName.data()), Form(
"# photons; %s # photons", locParticleROOTName.data()), 150, 0, 150);
57 hThetaC = GetOrCreate_Histogram<TH1I>(Form(
"hThetaC_%s",locParticleName.data()), Form(
"cherenkov angle; %s #theta_{C} [rad]", locParticleROOTName.data()), 250, 0.6, 1.0);
58 hDeltaThetaC = GetOrCreate_Histogram<TH1I>(Form(
"hDeltaThetaC_%s",locParticleName.data()), Form(
"cherenkov angle; %s #Delta#theta_{C} [rad]", locParticleROOTName.data()), 200,-0.2,0.2);
59 hLikelihood = GetOrCreate_Histogram<TH1I>(Form(
"hLikelihood_%s",locParticleName.data()), Form(
"; %s -lnL; entries [#]", locParticleROOTName.data()),1000,0.,1000.);
60 hLikelihoodDiff = GetOrCreate_Histogram<TH1I>(Form(
"hLikelihoodDiff_%s",locParticleName.data()), Form(
"; %s;entries [#]", locLikelihoodName.data()),100,-200.,200.);
62 hThetaCVsP = GetOrCreate_Histogram<TH2I>(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);
63 hDeltaThetaCVsP = GetOrCreate_Histogram<TH2I>(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);
64 hLikelihoodDiffVsP = GetOrCreate_Histogram<TH2I>(Form(
"hLikelihoodDiffVsP_%s",locParticleName.data()), Form(
"; p (GeV/c); %s", locLikelihoodName.data()), 120, 0.0, 12.0, 100, -200, 200);
65 hReactionLikelihoodDiffVsP = GetOrCreate_Histogram<TH2I>(Form(
"hReactionLikelihoodDiffVsP_%s",locParticleName.data()), Form(
"; p (GeV/c); %s", locLikelihoodName.data()), 120, 0.0, 12.0, 100, -200, 200);
70 for(
int locBar=0; locBar<24; locBar++) {
75 for(
int locXbin=0; locXbin<40; locXbin++) {
77 double xbin_min = -100.0 + locXbin*5.0;
78 double xbin_max = xbin_min + 5.0;
81 TVector3 locBinVect((xbin_min+xbin_max)/2., bar_y, 586.5 - 65.);
82 if(locBinVect.Theta()*180/TMath::Pi() > 12.)
85 hDiffMap[locBar][locXbin] = GetOrCreate_Histogram<TH1I>(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);
86 hHitTimeMap[locBar][locXbin] = GetOrCreate_Histogram<TH1I>(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);
87 hNphCMap[locBar][locXbin] = GetOrCreate_Histogram<TH1I>(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);
88 hNphCMapSlot4[locBar][locXbin] = GetOrCreate_Histogram<TH1I>(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);
89 hNphCMapSlot5[locBar][locXbin] = GetOrCreate_Histogram<TH1I>(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);
91 hDeltaThetaCVsPMap[locBar][locXbin] = GetOrCreate_Histogram<TH2I>(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);
92 hReactionLikelihoodDiffVsPMap[locBar][locXbin] = GetOrCreate_Histogram<TH2I>(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);
94 hPixelHitMap[locBar][locXbin] = GetOrCreate_Histogram<TH2S>(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);
95 hPixelHitMapReflected[locBar][locXbin] = GetOrCreate_Histogram<TH2S>(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);
97 hHitTimeMapAll[locBar][locXbin] = GetOrCreate_Histogram<TH1I>(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);
98 hPixelHitMapAll[locBar][locXbin] = GetOrCreate_Histogram<TH2S>(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);
99 hPixelHitMapAllReflected[locBar][locXbin] = GetOrCreate_Histogram<TH2S>(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);
103 gDirectory->cd(
"..");
117 locEventLoop->GetSingle(locDetectorMatches);
121 vector<const DDIRCTruthBarHit*> locDIRCBarHits;
122 locEventLoop->Get(locDIRCBarHits);
124 vector<const DDIRCPmtHit*> locDIRCPmtHits;
125 locEventLoop->Get(locDIRCPmtHits);
141 int locDCHits = locTrackTimeBased->
Ndof + 5;
142 double locTheta = locTrackTimeBased->
momentum().Theta()*180/TMath::Pi();
143 double locP = locParticle->lorentzMomentum().Vect().Mag();
144 if(locDCHits < 15 || locTheta < 1.0 || locTheta > 12.0 || locP > 12.0)
148 shared_ptr<const DTOFHitMatchParams> locTOFHitMatchParams;
150 if(!foundTOF || locTOFHitMatchParams->dDeltaXToHit > 10.0 || locTOFHitMatchParams->dDeltaYToHit > 10.0)
157 shared_ptr<const DDIRCMatchParams> locDIRCMatchParams;
159 vector<int> locUsedPixel;
163 DVector3 posInBar = locDIRCMatchParams->dExtrapolatedPos;
164 DVector3 momInBar = locDIRCMatchParams->dExtrapolatedMom;
165 double locExpectedThetaC = locDIRCMatchParams->dExpectedThetaC;
166 double locExtrapolatedTime = locDIRCMatchParams->dExtrapolatedTime;
170 int locXbin = (int)(posInBar.X()/5.0) + 19;
173 TVector3 locBarHitVect(posInBar.X(), posInBar.Y(), posInBar.Z() - 65.);
174 if(locBar < 0 || locBar >23 || locXbin < 0 || locXbin > 39 || locBarHitVect.Theta()*180/TMath::Pi() > 12.0)
185 map<shared_ptr<const DDIRCMatchParams>, vector<const DDIRCPmtHit*> > locDIRCTrackMatchParamsMap;
187 map<Particle_t, double> logLikelihoodSum;
189 int locNPhCSlot4 = 0;
190 int locNPhCSlot5 = 0;
193 for(uint loc_i=0; loc_i<locDIRCPmtHits.size(); loc_i++) {
194 vector<pair<double, double>> locDIRCPhotons =
dDIRCLut->
CalcPhoton(locDIRCPmtHits[loc_i], locExtrapolatedTime, posInBar, momInBar, locExpectedAngle, locAngle, locPID, logLikelihoodSum);
195 double locHitTime = locDIRCPmtHits[loc_i]->t - locExtrapolatedTime;
196 int locChannel = locDIRCPmtHits[loc_i]->ch;
197 if(locChannel >= 108*64) locChannel -= 108*64;
211 if(locHitTime > 0 && locHitTime < 100.0) {
212 if((pixel_row < 64 && pixel_col < 24) || (pixel_row < 32 && pixel_col > 23))
230 if(locDIRCPhotons.size() > 0) {
233 for(uint loc_j = 0; loc_j<locDIRCPhotons.size(); loc_j++) {
234 double locDeltaT = locDIRCPhotons[loc_j].first - locHitTime;
235 double locThetaC = locDIRCPhotons[loc_j].second;
238 hDiff->Fill(locDeltaT);
240 hDiffMap[locBar][locXbin]->Fill(locDeltaT);
246 if(fabs(locDeltaT) < 100.0) {
254 if(std::find(locUsedPixel.begin(), locUsedPixel.end(), locChannel) != locUsedPixel.end())
261 hPixelHitMap[locBar][locXbin]->Fill(pixel_row, pixel_col);
265 locUsedPixel.push_back(locChannel);
274 hNphC->Fill(locDIRCMatchParams->dNPhotons);
275 hThetaCVsP->Fill(momInBar.Mag(), locDIRCMatchParams->dThetaC);
277 hNphCMap[locBar][locXbin]->Fill(locDIRCMatchParams->dNPhotons);
284 hLikelihood->Fill(-1. * locDIRCMatchParams->dLikelihoodElectron);
285 hLikelihoodDiff->Fill(locDIRCMatchParams->dLikelihoodElectron - locDIRCMatchParams->dLikelihoodPion);
286 hLikelihoodDiffVsP->Fill(locP, locDIRCMatchParams->dLikelihoodElectron - locDIRCMatchParams->dLikelihoodPion);
292 hLikelihood->Fill(-1. * locDIRCMatchParams->dLikelihoodPion);
293 hLikelihoodDiff->Fill(locDIRCMatchParams->dLikelihoodPion - locDIRCMatchParams->dLikelihoodKaon);
294 hLikelihoodDiffVsP->Fill(locP, locDIRCMatchParams->dLikelihoodPion - locDIRCMatchParams->dLikelihoodKaon);
300 hLikelihood->Fill(-1. * locDIRCMatchParams->dLikelihoodKaon);
301 hLikelihoodDiff->Fill(locDIRCMatchParams->dLikelihoodPion - locDIRCMatchParams->dLikelihoodKaon);
302 hLikelihoodDiffVsP->Fill(locP, locDIRCMatchParams->dLikelihoodPion - locDIRCMatchParams->dLikelihoodKaon);
308 hLikelihood->Fill(-1. * locDIRCMatchParams->dLikelihoodProton);
309 hLikelihoodDiff->Fill(locDIRCMatchParams->dLikelihoodProton - locDIRCMatchParams->dLikelihoodKaon);
310 hLikelihoodDiffVsP->Fill(locP, locDIRCMatchParams->dLikelihoodKaon - locDIRCMatchParams->dLikelihoodProton);
const DAnalysisUtilities * dAnalysisUtilities
TDirectoryFile * ChangeTo_BaseDirectory(void)
vector< pair< double, double > > CalcPhoton(const DDIRCPmtHit *locDIRCHit, double locFlightTime, TVector3 posInBar, TVector3 momInBar, map< Particle_t, double > locExpectedAngle, double locAngle, Particle_t locPID, map< Particle_t, double > &logLikelihoodSum, int &nPhotonsThetaC, double &meanThetaC, double &meanDeltaT, bool &isGood) const
TH2I * hReactionLikelihoodDiffVsP
double GetBarY(int locBar) const
const DKinematicData * Get_FinalParticle_Measured(size_t locFinalParticleIndex) const
int GetPixelRow(int channel) const
static char * ParticleName_ROOT(Particle_t p)
double CalcAngle(TVector3 momInBar, double locMass) const
bool Get_DIRCTrackMatchParamsMap(map< shared_ptr< const DDIRCMatchParams >, vector< const DDIRCPmtHit * > > &locDIRCTrackMatchParamsMap)
bool Perform_Action(JEventLoop *locEventLoop, const DParticleCombo *locParticleCombo)
void Initialize(JEventLoop *locEventLoop)
TDirectoryFile * CreateAndChangeTo_Directory(TDirectoryFile *locBaseDirectory, string locDirName, string locDirTitle)
const DTrackTimeBased * Get_TrackTimeBased(void) const
TDirectoryFile * CreateAndChangeTo_ActionDirectory(void)
bool Get_DIRCMatchParams(const DTrackingData *locTrack, const DDetectorMatches *locDetectorMatches, shared_ptr< const DDIRCMatchParams > &locBestMatchParams) const
TH2I * hExtrapolatedBarHitXY
static char * ParticleType(Particle_t p)
TH2S * hPixelHitMapAllReflected[48][40]
int GetPixelColumn(int channel) const
const DChargedTrackHypothesis * Get_Hypothesis(Particle_t locPID) const
bool Get_UseKinFitResultsFlag(void) const
TH1I * hHitTimeMapAll[48][40]
TH1I * hNphCMapSlot5[48][40]
bool Get_BestTOFMatchParams(const DTrackingData *locTrack, const DDetectorMatches *locDetectorMatches, shared_ptr< const DTOFHitMatchParams > &locBestMatchParams) const
TH2S * hPixelHitMap[48][40]
TH2S * hPixelHitMapAll[48][40]
TH2I * hExtrapolatedBarHitXY_PreCut
const DDIRCGeometry * dDIRCGeometry
int Ndof
Number of degrees of freedom in the fit.
TH1I * hHitTimeMap[48][40]
static double ParticleMass(Particle_t p)
const JObject * Get_FinalParticle_SourceObject(size_t locFinalParticleIndex) const
TH1I * hNphCMapSlot4[48][40]
TH2S * hPixelHitMapReflected[48][40]
const DVector3 & momentum(void) const
const DDIRCLut * dDIRCLut
TH2I * hDeltaThetaCVsPMap[48][40]
TH2I * hReactionLikelihoodDiffVsPMap[48][40]
map< Particle_t, double > CalcExpectedAngles(TVector3 momInBar) const
const DParticleComboStep * Get_ParticleComboStep(size_t locStepIndex) const
const DKinematicData * Get_FinalParticle(size_t locFinalParticleIndex) const
deque< Particle_t > dFinalStatePIDs
TH2I * hLikelihoodDiffVsP
Particle_t PID(void) const
int dParticleComboStepIndex
const DParticleID * dParticleID
int GetBar(float y) const