1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | #include "GlueXSensitiveDetectorPSC.hh" |
8 | #include "GlueXDetectorConstruction.hh" |
9 | #include "GlueXPrimaryGeneratorAction.hh" |
10 | #include "GlueXUserEventInformation.hh" |
11 | #include "GlueXUserTrackInformation.hh" |
12 | #include "HddmOutput.hh" |
13 | |
14 | #include "G4VPhysicalVolume.hh" |
15 | #include "G4PVPlacement.hh" |
16 | #include "G4EventManager.hh" |
17 | #include "G4HCofThisEvent.hh" |
18 | #include "G4Step.hh" |
19 | #include "G4SDManager.hh" |
20 | #include "G4ios.hh" |
21 | |
22 | #include <JANA/JApplication.h> |
23 | |
24 | |
25 | int GlueXSensitiveDetectorPSC::MAX_HITS = 100; |
26 | |
27 | |
28 | double GlueXSensitiveDetectorPSC::TWO_HIT_TIME_RESOL = 25*ns; |
29 | |
30 | |
31 | double GlueXSensitiveDetectorPSC::THRESH_MEV = 0.010; |
32 | |
33 | |
34 | int GlueXSensitiveDetectorPSC::NUM_MODULES_PER_ARM = 8; |
35 | |
36 | int GlueXSensitiveDetectorPSC::instanceCount = 0; |
37 | G4Mutex GlueXSensitiveDetectorPSC::fMutex = G4MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } }; |
38 | |
39 | GlueXSensitiveDetectorPSC::GlueXSensitiveDetectorPSC(const G4String& name) |
40 | : G4VSensitiveDetector(name), |
41 | fCounterHitsMap(0), fPointsMap(0) |
42 | { |
43 | collectionName.insert("PSCPaddleHitsCollection"); |
44 | collectionName.insert("PSCPointsCollection"); |
45 | |
46 | |
47 | |
48 | |
49 | |
50 | |
51 | |
52 | G4AutoLock barrier(&fMutex); |
53 | if (instanceCount++ == 0) { |
54 | int runno = HddmOutput::getRunNo(); |
55 | extern jana::JApplication *japp; |
56 | if (japp == 0) { |
57 | G4cerr(*G4cerr_p) << "Error in GlueXSensitiveDetector constructor - " |
58 | << "jana global DApplication object not set, " |
59 | << "cannot continue." << G4endlstd::endl; |
60 | exit(-1); |
61 | } |
62 | jana::JCalibration *jcalib = japp->GetJCalibration(runno); |
| Value stored to 'jcalib' during its initialization is never read |
63 | if (japp == 0) { |
64 | jcalib = 0; |
65 | G4cout(*G4cout_p) << "PSC: ALL parameters loaded from ccdb" << G4endlstd::endl; |
66 | } |
67 | } |
68 | } |
69 | |
70 | GlueXSensitiveDetectorPSC::GlueXSensitiveDetectorPSC( |
71 | const GlueXSensitiveDetectorPSC &src) |
72 | : G4VSensitiveDetector(src), |
73 | fCounterHitsMap(src.fCounterHitsMap), fPointsMap(src.fPointsMap) |
74 | { |
75 | G4AutoLock barrier(&fMutex); |
76 | ++instanceCount; |
77 | } |
78 | |
79 | GlueXSensitiveDetectorPSC &GlueXSensitiveDetectorPSC::operator=(const |
80 | GlueXSensitiveDetectorPSC &src) |
81 | { |
82 | G4AutoLock barrier(&fMutex); |
83 | *(G4VSensitiveDetector*)this = src; |
84 | fCounterHitsMap = src.fCounterHitsMap; |
85 | fPointsMap = src.fPointsMap; |
86 | return *this; |
87 | } |
88 | |
89 | GlueXSensitiveDetectorPSC::~GlueXSensitiveDetectorPSC() |
90 | { |
91 | G4AutoLock barrier(&fMutex); |
92 | --instanceCount; |
93 | } |
94 | |
95 | void GlueXSensitiveDetectorPSC::Initialize(G4HCofThisEvent* hce) |
96 | { |
97 | fCounterHitsMap = new |
98 | GlueXHitsMapPSCpaddle(SensitiveDetectorName, collectionName[0]); |
99 | fPointsMap = new |
100 | GlueXHitsMapPSCpoint(SensitiveDetectorName, collectionName[1]); |
101 | G4SDManager *sdm = G4SDManager::GetSDMpointer(); |
102 | hce->AddHitsCollection(sdm->GetCollectionID(collectionName[0]), fCounterHitsMap); |
103 | hce->AddHitsCollection(sdm->GetCollectionID(collectionName[1]), fPointsMap); |
104 | } |
105 | |
106 | G4bool GlueXSensitiveDetectorPSC::ProcessHits(G4Step* step, |
107 | G4TouchableHistory* ROhist) |
108 | { |
109 | double dEsum = step->GetTotalEnergyDeposit(); |
110 | if (dEsum == 0) |
111 | return false; |
112 | |
113 | const G4ThreeVector &pin = step->GetPreStepPoint()->GetMomentum(); |
114 | const G4ThreeVector &xin = step->GetPreStepPoint()->GetPosition(); |
115 | const G4ThreeVector &xout = step->GetPostStepPoint()->GetPosition(); |
116 | double Ein = step->GetPreStepPoint()->GetTotalEnergy(); |
117 | double tin = step->GetPreStepPoint()->GetGlobalTime(); |
118 | double tout = step->GetPostStepPoint()->GetGlobalTime(); |
119 | G4ThreeVector x = (xin + xout) / 2; |
120 | G4ThreeVector dx = xout - xin; |
121 | double t = (tin + tout) / 2; |
122 | |
123 | const G4VTouchable* touch = step->GetPreStepPoint()->GetTouchable(); |
124 | const G4AffineTransform &local_from_global = touch->GetHistory() |
125 | ->GetTopTransform(); |
126 | G4ThreeVector xlocal = local_from_global.TransformPoint(x); |
127 | |
128 | |
129 | |
130 | |
131 | |
132 | |
133 | |
134 | if (tout > 1.0*s) |
135 | t = tin; |
136 | |
137 | double dr = dx.mag(); |
138 | double dEdx = (dr > 1e-3*cm)? dEsum/dr : 0; |
139 | |
140 | |
141 | |
142 | |
143 | int module = GetIdent("module", touch); |
144 | int arm = (module - 1) / NUM_MODULES_PER_ARM; |
145 | module = (module - 1) % NUM_MODULES_PER_ARM + 1; |
146 | G4Track *track = step->GetTrack(); |
147 | G4int trackID = track->GetTrackID(); |
148 | int pdgtype = track->GetDynamicParticle()->GetPDGcode(); |
149 | int g3type = GlueXPrimaryGeneratorAction::ConvertPdgToGeant3(pdgtype); |
150 | GlueXUserTrackInformation *trackinfo = (GlueXUserTrackInformation*) |
151 | track->GetUserInformation(); |
152 | int itrack = trackinfo->GetGlueXTrackID(); |
153 | if (trackinfo->GetGlueXHistory() == 0 && itrack > 0) { |
154 | G4int key = fPointsMap->entries(); |
155 | GlueXHitPSCpoint* lastPoint = (*fPointsMap)[key - 1]; |
156 | if (lastPoint == 0 || lastPoint->track_ != trackID || |
157 | fabs(lastPoint->t_ns - t/ns) > 0.2 || |
158 | fabs(lastPoint->x_cm - x[0]/cm) > 5.0 || |
159 | fabs(lastPoint->y_cm - x[1]/cm) > 5.0 || |
160 | fabs(lastPoint->z_cm - x[2]/cm) > 5.0) |
161 | { |
162 | GlueXHitPSCpoint newPoint; |
163 | newPoint.arm_ = arm; |
164 | newPoint.module_ = module; |
165 | newPoint.ptype_G3 = g3type; |
166 | newPoint.track_ = trackID; |
167 | newPoint.trackID_ = itrack; |
168 | newPoint.primary_ = (track->GetParentID() == 0); |
169 | newPoint.t_ns = t/ns; |
170 | newPoint.x_cm = x[0]/cm; |
171 | newPoint.y_cm = x[1]/cm; |
172 | newPoint.z_cm = x[2]/cm; |
173 | newPoint.px_GeV = pin[0]/GeV; |
174 | newPoint.py_GeV = pin[1]/GeV; |
175 | newPoint.pz_GeV = pin[2]/GeV; |
176 | newPoint.E_GeV = Ein/GeV; |
177 | newPoint.dEdx_GeV_cm = dEdx/(GeV/cm); |
178 | fPointsMap->add(key, newPoint); |
179 | } |
180 | } |
181 | |
182 | |
183 | |
184 | if (dEsum > 0) { |
185 | int key = GlueXHitPSCpaddle::GetKey(arm, module); |
186 | GlueXHitPSCpaddle *paddle = (*fCounterHitsMap)[key]; |
187 | if (paddle == 0) { |
188 | GlueXHitPSCpaddle newpaddle(arm, module); |
189 | fCounterHitsMap->add(key, newpaddle); |
190 | paddle = (*fCounterHitsMap)[key]; |
191 | } |
192 | |
193 | |
194 | |
195 | int merge_hit = 0; |
196 | std::vector<GlueXHitPSCpaddle::hitinfo_t>::iterator hiter; |
197 | for (hiter = paddle->hits.begin(); hiter != paddle->hits.end(); ++hiter) { |
198 | if (fabs(hiter->t_ns*ns - t) < TWO_HIT_TIME_RESOL) { |
199 | merge_hit = 1; |
200 | break; |
201 | } |
202 | else if (hiter->t_ns*ns > t) { |
203 | break; |
204 | } |
205 | } |
206 | if (merge_hit) { |
207 | |
208 | hiter->t_ns = (hiter->t_ns * hiter->dE_GeV + t/ns * dEsum/GeV) / |
209 | (hiter->dE_GeV + dEsum/GeV); |
210 | hiter->dE_GeV += dEsum/GeV; |
211 | } |
212 | else if ((int)paddle->hits.size() < MAX_HITS) { |
213 | |
214 | hiter = paddle->hits.insert(hiter, GlueXHitPSCpaddle::hitinfo_t()); |
215 | hiter->dE_GeV = dEsum/GeV; |
216 | hiter->t_ns = t/ns; |
217 | hiter->itrack_ = itrack; |
218 | hiter->ptype_G3 = g3type; |
219 | } |
220 | else { |
221 | G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorPSC::ProcessHits error: " |
222 | << "max hit count " << MAX_HITS << " exceeded, truncating!" |
223 | << G4endlstd::endl; |
224 | } |
225 | } |
226 | return true; |
227 | } |
228 | |
229 | void GlueXSensitiveDetectorPSC::EndOfEvent(G4HCofThisEvent*) |
230 | { |
231 | std::map<int,GlueXHitPSCpaddle*> *paddles = fCounterHitsMap->GetMap(); |
232 | std::map<int,GlueXHitPSCpoint*> *points = fPointsMap->GetMap(); |
233 | if (paddles->size() == 0 && points->size() == 0) |
234 | return; |
235 | std::map<int,GlueXHitPSCpaddle*>::iterator siter; |
236 | std::map<int,GlueXHitPSCpoint*>::iterator piter; |
237 | |
238 | if (verboseLevel > 1) { |
239 | G4cout(*G4cout_p) << G4endlstd::endl |
240 | << "--------> Hits Collection: in this event there are " |
241 | << paddles->size() << " paddles with hits in the PSC: " |
242 | << G4endlstd::endl; |
243 | for (siter = paddles->begin(); siter != paddles->end(); ++siter) |
244 | siter->second->Print(); |
245 | |
246 | G4cout(*G4cout_p) << G4endlstd::endl |
247 | << "--------> Hits Collection: in this event there are " |
248 | << points->size() << " truth points in the PSC: " |
249 | << G4endlstd::endl; |
250 | for (piter = points->begin(); piter != points->end(); ++piter) |
251 | piter->second->Print(); |
252 | } |
253 | |
254 | |
255 | |
256 | G4EventManager* mgr = G4EventManager::GetEventManager(); |
257 | G4VUserEventInformation* info = mgr->GetUserInformation(); |
258 | hddm_s::HDDM *record = ((GlueXUserEventInformation*)info)->getOutputRecord(); |
259 | if (record == 0) { |
260 | G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorPSC::EndOfEvent error - " |
261 | << "hits seen but no output hddm record to save them into, " |
262 | << "cannot continue!" << G4endlstd::endl; |
263 | exit(1); |
264 | } |
265 | |
266 | if (record->getPhysicsEvents().size() == 0) |
267 | record->addPhysicsEvents(); |
268 | if (record->getHitViews().size() == 0) |
269 | record->getPhysicsEvent().addHitViews(); |
270 | hddm_s::HitView &hitview = record->getPhysicsEvent().getHitView(); |
271 | if (hitview.getPairSpectrometerCoarses().size() == 0) |
272 | hitview.addPairSpectrometerCoarses(); |
273 | hddm_s::PairSpectrometerCoarse &psc = hitview.getPairSpectrometerCoarse(); |
274 | |
275 | |
276 | for (siter = paddles->begin(); siter != paddles->end(); ++siter) { |
277 | std::vector<GlueXHitPSCpaddle::hitinfo_t> &hits = siter->second->hits; |
278 | |
279 | for (unsigned int ih=0; ih < hits.size(); ++ih) { |
280 | if (hits[ih].dE_GeV*1e3 <= THRESH_MEV) { |
281 | hits.erase(hits.begin() + ih); |
282 | --ih; |
283 | } |
284 | } |
285 | if (hits.size() > 0) { |
286 | hddm_s::PscPaddleList paddle = psc.addPscPaddles(1); |
287 | paddle(0).setArm(siter->second->arm_); |
288 | paddle(0).setModule(siter->second->module_); |
289 | for (int ih=0; ih < (int)hits.size(); ++ih) { |
290 | hddm_s::PscTruthHitList thit = paddle(0).addPscTruthHits(1); |
291 | thit(0).setDE(hits[ih].dE_GeV); |
292 | thit(0).setT(hits[ih].t_ns); |
293 | thit(0).setItrack(hits[ih].itrack_); |
294 | thit(0).setPtype(hits[ih].ptype_G3); |
295 | } |
296 | } |
297 | } |
298 | |
299 | |
300 | for (piter = points->begin(); piter != points->end(); ++piter) { |
301 | hddm_s::PscTruthPointList point = psc.addPscTruthPoints(1); |
302 | point(0).setE(piter->second->E_GeV); |
303 | point(0).setDEdx(piter->second->dEdx_GeV_cm); |
304 | point(0).setPrimary(piter->second->primary_); |
305 | point(0).setPtype(piter->second->ptype_G3); |
306 | point(0).setPx(piter->second->px_GeV); |
307 | point(0).setPy(piter->second->py_GeV); |
308 | point(0).setPz(piter->second->pz_GeV); |
309 | point(0).setX(piter->second->x_cm); |
310 | point(0).setY(piter->second->y_cm); |
311 | point(0).setZ(piter->second->z_cm); |
312 | point(0).setT(piter->second->t_ns); |
313 | point(0).setArm(piter->second->arm_); |
314 | point(0).setModule(piter->second->module_); |
315 | point(0).setTrack(piter->second->track_); |
316 | hddm_s::TrackIDList tid = point(0).addTrackIDs(); |
317 | tid(0).setItrack(piter->second->trackID_); |
318 | } |
319 | } |
320 | |
321 | int GlueXSensitiveDetectorPSC::GetIdent(std::string div, |
322 | const G4VTouchable *touch) |
323 | { |
324 | const HddsG4Builder* bldr = GlueXDetectorConstruction::GetBuilder(); |
325 | std::map<std::string, std::vector<int> >::const_iterator iter; |
326 | std::map<std::string, std::vector<int> > *identifiers; |
327 | int max_depth = touch->GetHistoryDepth(); |
328 | for (int depth = 0; depth < max_depth; ++depth) { |
329 | G4VPhysicalVolume *pvol = touch->GetVolume(depth); |
330 | G4LogicalVolume *lvol = pvol->GetLogicalVolume(); |
331 | int volId = fVolumeTable[lvol]; |
332 | if (volId == 0) { |
333 | volId = bldr->getVolumeId(lvol); |
334 | fVolumeTable[lvol] = volId; |
335 | } |
336 | identifiers = &Refsys::fIdentifierTable[volId]; |
337 | if ((iter = identifiers->find(div)) != identifiers->end()) { |
338 | int copyNum = touch->GetCopyNumber(depth); |
339 | copyNum += (dynamic_cast<G4PVPlacement*>(pvol))? -1 : 0; |
340 | return iter->second[copyNum]; |
341 | } |
342 | } |
343 | return -1; |
344 | } |