diff --git a/src/GlueXHitBCALcell.cc b/src/GlueXHitBCALcell.cc index 732fc8c..5a7c2d8 100644 --- a/src/GlueXHitBCALcell.cc +++ b/src/GlueXHitBCALcell.cc @@ -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; @@ -82,9 +86,13 @@ void GlueXHitBCALcell::Print() const std::vector::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; } } diff --git a/src/GlueXHitBCALcell.hh b/src/GlueXHitBCALcell.hh index e614931..4549836 100644 --- a/src/GlueXHitBCALcell.hh +++ b/src/GlueXHitBCALcell.hh @@ -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 hits; diff --git a/src/GlueXSensitiveDetectorBCAL.cc b/src/GlueXSensitiveDetectorBCAL.cc index 32260d9..ebb2ffb 100644 --- a/src/GlueXSensitiveDetectorBCAL.cc +++ b/src/GlueXSensitiveDetectorBCAL.cc @@ -20,6 +20,8 @@ #include +#define USE_ENERGY_WEIGHTED_TIMES 1 + // Cutoff on the total number of allowed hits int GlueXSensitiveDetectorBCAL::MAX_HITS = 100; @@ -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::iterator hiter; @@ -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(); @@ -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; } @@ -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); + } } } }