Skip to content

Commit

Permalink
Merge pull request #100 from JeffersonLab/bcal_timing_fixes_rtj
Browse files Browse the repository at this point in the history 
* change bcal shower timing to use energy-weighted times,
  • Loading branch information
@sdobbs
sdobbs authored Apr 4, 2019
2 parents 02f6a22 + 6bc99ee commit 16b2bc9
Show file tree
Hide file tree
Showing 3 changed files with 123 additions and 4 deletions.
10 changes: 9 additions & 1 deletion src/GlueXHitBCALcell.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,10 @@ int GlueXHitBCALcell::operator==(const GlueXHitBCALcell &right) const
if (hits[ih].E_GeV != right.hits[ih].E_GeV ||
hits[ih].zlocal_cm != right.hits[ih].zlocal_cm ||
hits[ih].t_ns != right.hits[ih].t_ns ||
hits[ih].Eup_GeV != right.hits[ih].Eup_GeV ||
hits[ih].Edown_GeV != right.hits[ih].Edown_GeV ||
hits[ih].tup_ns != right.hits[ih].tup_ns ||
hits[ih].tdown_ns != right.hits[ih].tdown_ns ||
hits[ih].incidentId_ != right.hits[ih].incidentId_ )
{
return 0;
Expand Down Expand Up @@ -82,9 +86,13 @@ void GlueXHitBCALcell::Print() const
std::vector<hitinfo_t>::const_iterator hiter;
for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
G4cout << " E = " << hiter->E_GeV << " GeV" << G4endl
<< " zlocal = " << hiter->t_ns << " cm" << G4endl
<< " t = " << hiter->t_ns << " ns" << G4endl
<< " zlocal = " << hiter->t_ns << " cm" << G4endl
<< " incidentId = " << hiter->incidentId_ << G4endl
<< " E(upstream end) = " << hiter->Eup_GeV << " GeV" << G4endl
<< " E(downstream end) = " << hiter->Edown_GeV << " GeV" << G4endl
<< " t(upstream end) = " << hiter->tup_ns << " ns" << G4endl
<< " t(downstream end) = " << hiter->tdown_ns << " ns" << G4endl
<< G4endl;
}
}
Expand Down
4 changes: 4 additions & 0 deletions src/GlueXHitBCALcell.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,10 @@ class GlueXHitBCALcell : public G4VHit
G4double t_ns; // pulse leading-edge time (ns)
G4double zlocal_cm; // z coordinate of the hit in local refsys
G4double incidentId_; // id of particle that generated this shower
G4double Eup_GeV; // upstream end energy deposition (GeV)
G4double Edown_GeV; // downstream end energy deposition (GeV)
G4double tup_ns; // upstream end pulse leading-edge time (ns)
G4double tdown_ns; // downstream end pulse leading-edge time (ns)
};
std::vector<hitinfo_t> hits;

Expand Down
113 changes: 110 additions & 3 deletions src/GlueXSensitiveDetectorBCAL.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,8 @@

#include <JANA/JApplication.h>

#define USE_ENERGY_WEIGHTED_TIMES 1

// Cutoff on the total number of allowed hits
int GlueXSensitiveDetectorBCAL::MAX_HITS = 100;

Expand Down Expand Up @@ -234,7 +236,7 @@ G4bool GlueXSensitiveDetectorBCAL::ProcessHits(G4Step* step,
cell = (*fCellsMap)[key];
}

// Add the hit to the hits vector, maintaining strict time ordering
// Add the hit to the bcal truth hits list, maintaining strict time ordering

int merge_hit = 0;
std::vector<GlueXHitBCALcell::hitinfo_t>::iterator hiter;
Expand All @@ -249,17 +251,28 @@ G4bool GlueXSensitiveDetectorBCAL::ProcessHits(G4Step* step,
}
if (merge_hit) {
// Use the time from the earlier hit but add the energy deposition
hiter->E_GeV += dEsum/GeV;
#if USE_ENERGY_WEIGHTED_TIMES
hiter->t_ns = (hiter->t_ns * hiter->E_GeV +
t/ns * dEsum/GeV) /
(hiter->E_GeV + dEsum/GeV);
hiter->zlocal_cm = (hiter->zlocal_cm * hiter->E_GeV +
xlocal[2]/cm * dEsum/GeV) /
(hiter->E_GeV + dEsum/GeV);
#else
if (hiter->t_ns*ns > t) {
hiter->t_ns = t/ns;
hiter->zlocal_cm = xlocal[2]/cm;
hiter->incidentId_ = trackinfo->GetBCALincidentID();
}
#endif
// correction factor makes shower yields match hdgeant
hiter->E_GeV += dEsum/GeV * 1.015;
}
else if ((int)cell->hits.size() < MAX_HITS) {
// create new hit
hiter = cell->hits.insert(hiter, GlueXHitBCALcell::hitinfo_t());
hiter->E_GeV = dEsum/GeV;
// correction factor makes shower yields match hdgeant
hiter->E_GeV = dEsum/GeV * 1.015;
hiter->t_ns = t/ns;
hiter->zlocal_cm = xlocal[2]/cm;
hiter->incidentId_ = trackinfo->GetBCALincidentID();
Expand All @@ -269,6 +282,90 @@ G4bool GlueXSensitiveDetectorBCAL::ProcessHits(G4Step* step,
<< "max hit count " << MAX_HITS << " exceeded, truncating!"
<< G4endl;
}

// Add the hit to the upstream sipm hits list, with strict time ordering

double udist = MODULE_FULL_LENGTH/2 + xlocal[2];
double dEup = dEsum * exp(-udist/ATTENUATION_LENGTH);
double tup = t + udist / C_EFFECTIVE;
merge_hit = 0;
for (hiter = cell->hits.begin(); hiter != cell->hits.end(); ++hiter) {
if (fabs(hiter->tup_ns*ns - tup) < TWO_HIT_TIME_RESOL) {
merge_hit = 1;
break;
}
else if (hiter->tup_ns*ns > tup) {
break;
}
}
if (merge_hit) {
// Use the time from the earlier hit but add the energy deposition
#if USE_ENERGY_WEIGHTED_TIMES
hiter->tup_ns = (hiter->tup_ns * hiter->Eup_GeV +
tup/ns * dEup/GeV) /
(hiter->Eup_GeV + dEup/GeV);
#else
if (hiter->tup_ns*ns > tup) {
hiter->tup_ns = tup/ns;
}
#endif
// correction factor makes shower yields match hdgeant
hiter->Eup_GeV += dEup/GeV * 1.015;
}
else if ((int)cell->hits.size() < MAX_HITS) {
// create new hit
hiter = cell->hits.insert(hiter, GlueXHitBCALcell::hitinfo_t());
// correction factor makes shower yields match hdgeant
hiter->Eup_GeV = dEup/GeV * 1.015;
hiter->tup_ns = tup/ns;
}
else {
G4cerr << "GlueXSensitiveDetectorBCAL::ProcessHits error: "
<< "max hit count " << MAX_HITS << " exceeded, truncating!"
<< G4endl;
}

// Add the hit to the downstream sipm hits list, with strict time ordering

double ddist = MODULE_FULL_LENGTH/2 - xlocal[2];
double dEdown = dEsum * exp(-ddist/ATTENUATION_LENGTH);
double tdown = t + ddist / C_EFFECTIVE;
merge_hit = 0;
for (hiter = cell->hits.begin(); hiter != cell->hits.end(); ++hiter) {
if (fabs(hiter->tdown_ns*ns - tdown) < TWO_HIT_TIME_RESOL) {
merge_hit = 1;
break;
}
else if (hiter->tdown_ns*ns > tdown) {
break;
}
}
if (merge_hit) {
// Use the time from the earlier hit but add the energy deposition
#if USE_ENERGY_WEIGHTED_TIMES
hiter->tdown_ns = (hiter->tdown_ns * hiter->Edown_GeV +
tdown/ns * dEdown/GeV) /
(hiter->Edown_GeV + dEdown/GeV);
#else
if (hiter->tdown_ns*ns > tdown) {
hiter->tdown_ns = tdown/ns;
}
#endif
// correction factor makes shower yields match hdgeant
hiter->Edown_GeV += dEdown/GeV * 1.015;
}
else if ((int)cell->hits.size() < MAX_HITS) {
// create new hit
hiter = cell->hits.insert(hiter, GlueXHitBCALcell::hitinfo_t());
// correction factor makes shower yields match hdgeant
hiter->Edown_GeV = dEdown/GeV * 1.015;
hiter->tdown_ns = tdown/ns;
}
else {
G4cerr << "GlueXSensitiveDetectorBCAL::ProcessHits error: "
<< "max hit count " << MAX_HITS << " exceeded, truncating!"
<< G4endl;
}
}
return true;
}
Expand Down Expand Up @@ -340,6 +437,16 @@ void GlueXSensitiveDetectorBCAL::EndOfEvent(G4HCofThisEvent*)
thit(0).setT(hits[ih].t_ns);
thit(0).setZLocal(hits[ih].zlocal_cm);
thit(0).setIncident_id(hits[ih].incidentId_);
if (hits[ih].Eup_GeV >= THRESH_MEV/1e3) {
hddm_s::BcalSiPMUpHitList uhit = cell(0).addBcalSiPMUpHits(1);
uhit(0).setE(hits[ih].Eup_GeV);
uhit(0).setT(hits[ih].tup_ns);
}
if (hits[ih].Edown_GeV >= THRESH_MEV/1e3) {
hddm_s::BcalSiPMDownHitList dhit = cell(0).addBcalSiPMDownHits(1);
dhit(0).setE(hits[ih].Edown_GeV);
dhit(0).setT(hits[ih].tdown_ns);
}
}
}
}
Expand Down

0 comments on commit 16b2bc9

Please sign in to comment.