1 | #include "ELYSE/EventAction.hh" |
---|
2 | |
---|
3 | //Geant4 |
---|
4 | #include "G4Event.hh" |
---|
5 | #include "G4RunManager.hh" |
---|
6 | #include "G4EventManager.hh" |
---|
7 | #include "G4UImanager.hh" |
---|
8 | #include "G4TrajectoryContainer.hh" |
---|
9 | #include "G4VVisManager.hh" |
---|
10 | #include "G4ios.hh" |
---|
11 | #include "globals.hh" |
---|
12 | #include "G4ThreeVector.hh" |
---|
13 | #include "G4TransportationManager.hh" |
---|
14 | #include "G4Navigator.hh" |
---|
15 | #include "G4SDManager.hh" |
---|
16 | #include "G4DigiManager.hh" |
---|
17 | #include "G4UnitsTable.hh" |
---|
18 | #include "G4UIcmdWith3VectorAndUnit.hh" |
---|
19 | |
---|
20 | //std |
---|
21 | #include <set> |
---|
22 | #include <iomanip> |
---|
23 | #include <string> |
---|
24 | #include <vector> |
---|
25 | #include <stdlib.h> //abs |
---|
26 | |
---|
27 | // AIDA : |
---|
28 | #include <AIDA/ITree.h> |
---|
29 | #include <AIDA/IManagedObject.h> |
---|
30 | #include <AIDA/ITuple.h> |
---|
31 | |
---|
32 | //ELYSE |
---|
33 | #include "ELYSE/Cast.hh" |
---|
34 | #include "ELYSE/Analysis.hh" |
---|
35 | #include "ELYSE/Trajectory.hh" |
---|
36 | #include "ELYSE/RunAction.hh" |
---|
37 | #include "ELYSE/PrimaryGeneratorAction.hh" |
---|
38 | //#include "ELYSE/TrappVolHit.hh" |
---|
39 | #include "ELYSE/DetectorConstruction.hh" |
---|
40 | |
---|
41 | |
---|
42 | ELYSE::EventAction::EventAction(ELYSE::Analysis& aAnalysis, |
---|
43 | ELYSE::RunAction* myRun, |
---|
44 | ELYSE::DetectorConstruction* myDetector, |
---|
45 | ELYSE::PrimaryGeneratorAction* myGenerator) |
---|
46 | :fAnalysis(aAnalysis), |
---|
47 | runAction(myRun), |
---|
48 | generatorAction(myGenerator), |
---|
49 | detectorConstructor(myDetector), |
---|
50 | eventTuple(0) { |
---|
51 | |
---|
52 | //Get User Histo pointers |
---|
53 | AIDA::ITree* usrTree = aAnalysis.tree(); //get the tree |
---|
54 | if (!usrTree) { |
---|
55 | G4cout << "ELYSE::EventAction: cannot get Analysis Tree" << G4endl; |
---|
56 | return; |
---|
57 | } |
---|
58 | //see AIDA::IManagedObject a comment on dynamic_cast |
---|
59 | AIDA::IManagedObject* obj; //generic Object managed by a Tree |
---|
60 | |
---|
61 | obj = usrTree->find("Event"); |
---|
62 | if(!obj) { |
---|
63 | G4cout << "EventAction: WARNING: no tuple Event" << G4endl; |
---|
64 | usrTree->ls(); |
---|
65 | //exit(0); |
---|
66 | } else { |
---|
67 | eventTuple = CAST(obj,AIDA::ITuple); |
---|
68 | if (!eventTuple) { |
---|
69 | G4cout << "EventAction: FATAL: eventTuple not a Tuple" << G4endl; |
---|
70 | usrTree->ls(); |
---|
71 | exit(0); |
---|
72 | } |
---|
73 | } |
---|
74 | |
---|
75 | }//Ctor |
---|
76 | |
---|
77 | //------------------------------------------------------------------------------------------- |
---|
78 | |
---|
79 | ELYSE::EventAction::~EventAction(){ |
---|
80 | }//Dtor |
---|
81 | |
---|
82 | //------------------------------------------------------------------------------------------- |
---|
83 | |
---|
84 | void ELYSE::EventAction::BeginOfEventAction(const G4Event* evt){ |
---|
85 | |
---|
86 | G4cout << " (JEC) EventAction::Begin EventAction" << G4endl; |
---|
87 | if (!evt) evt->Print(); |
---|
88 | |
---|
89 | }//BeginOfEventAction |
---|
90 | |
---|
91 | //------------------------------------------------------------------------------------------- |
---|
92 | |
---|
93 | void ELYSE::EventAction::EndOfEventAction(const G4Event* evt) { |
---|
94 | G4cout << " (JEC) EventAction::End EventAction" << G4endl; |
---|
95 | |
---|
96 | //see Analysis.cxx to get the Tuple variables |
---|
97 | |
---|
98 | // -------------------- |
---|
99 | // Get Event Information |
---|
100 | // -------------------- |
---|
101 | //JEC FIXME: difference between event_id & vecRecNumber |
---|
102 | G4int event_id = evt->GetEventID(); |
---|
103 | |
---|
104 | G4ThreeVector vtx = generatorAction->GetVtx(); |
---|
105 | G4int vtxvol = EventFindStartingVolume(vtx); |
---|
106 | |
---|
107 | G4cout << ">>>>>>>>>>>>> New Event "; |
---|
108 | G4cout << " evt Id " << event_id |
---|
109 | << " start in volume " << vtxvol |
---|
110 | << G4endl; |
---|
111 | G4cout << "Vertex (" |
---|
112 | << vtx.x()/cm << " , " |
---|
113 | << vtx.y()/cm << " , " |
---|
114 | << vtx.z()/cm << " , " |
---|
115 | << ")" |
---|
116 | << G4endl; |
---|
117 | |
---|
118 | //JEC FIXME: introduce enumeration for the column shared by Analysis/EventAction & namespace protected |
---|
119 | |
---|
120 | if(!eventTuple) return; |
---|
121 | |
---|
122 | eventTuple->fill(0, event_id); //eventId |
---|
123 | eventTuple->fill(1, vtxvol); //vtxVol |
---|
124 | |
---|
125 | AIDA::ITuple* vtxPos = eventTuple->getTuple( 2 ); |
---|
126 | vtxPos->fill(0, vtx.x()/cm); //vtxPos |
---|
127 | vtxPos->fill(1, vtx.y()/cm); |
---|
128 | vtxPos->fill(2, vtx.z()/cm); |
---|
129 | vtxPos->addRow(); |
---|
130 | |
---|
131 | |
---|
132 | // -------------------- |
---|
133 | // Get Track Information |
---|
134 | // -------------------- |
---|
135 | |
---|
136 | //variables used to fill the tuple (track part) |
---|
137 | G4int pId, parent; |
---|
138 | G4float timeStart; |
---|
139 | G4double dx, dy, dz; |
---|
140 | G4double ETot, px, py, pz; |
---|
141 | G4int startVol, stopVol; |
---|
142 | |
---|
143 | //---------------- |
---|
144 | // The beam: 1 particle for the moment |
---|
145 | //---------------- |
---|
146 | AIDA::ITuple* track = eventTuple->getTuple( 4 ); |
---|
147 | |
---|
148 | G4int beampdg = generatorAction->GetBeamPDG(); |
---|
149 | pId = beampdg; //pId |
---|
150 | track->fill(0, pId); |
---|
151 | |
---|
152 | parent = 0; //parent (none) |
---|
153 | track->fill(1, parent); |
---|
154 | |
---|
155 | timeStart = 0; //creation time |
---|
156 | track->fill(2, timeStart); |
---|
157 | |
---|
158 | G4ThreeVector beamdir = generatorAction->GetBeamDir(); //direction |
---|
159 | AIDA::ITuple* direction = track->getTuple( 3 ); |
---|
160 | dx = beamdir.x(); |
---|
161 | dy = beamdir.y(); |
---|
162 | dz = beamdir.z(); |
---|
163 | |
---|
164 | direction->fill(0, dx); |
---|
165 | direction->fill(1, dy); |
---|
166 | direction->fill(2, dz); |
---|
167 | direction->addRow(); |
---|
168 | |
---|
169 | |
---|
170 | G4double beamenergy = generatorAction->GetBeamEnergy(); |
---|
171 | G4double mass = G4ParticleTable::GetParticleTable()->FindParticle(beampdg)->GetPDGMass(); |
---|
172 | |
---|
173 | track->fill(4, mass); |
---|
174 | |
---|
175 | ETot = beamenergy; |
---|
176 | G4double pTot = ETot; |
---|
177 | if ( (beampdg != 22) && //gamma & Optical Photon (To be verified) |
---|
178 | (std::abs(beampdg) != 12) && //(anti)nu-e |
---|
179 | (std::abs(beampdg) != 14) && //(anti)nu-mu |
---|
180 | (std::abs(beampdg) != 16) ) { //(anti)nu-tau |
---|
181 | pTot = ETot*ETot - mass*mass; |
---|
182 | if (pTot < 1e-6){ |
---|
183 | pTot = 0.; |
---|
184 | } else { |
---|
185 | pTot = std::sqrt(pTot); |
---|
186 | } |
---|
187 | } |
---|
188 | |
---|
189 | track->fill(5, pTot); |
---|
190 | track->fill(6, ETot); |
---|
191 | |
---|
192 | |
---|
193 | AIDA::ITuple* momentum = track->getTuple( 7 ); //momentum |
---|
194 | px = pTot * dx; |
---|
195 | py = pTot * dy; |
---|
196 | pz = pTot * dz; |
---|
197 | momentum->fill(0, px); |
---|
198 | momentum->fill(1, py); |
---|
199 | momentum->fill(2, pz); |
---|
200 | momentum->addRow(); |
---|
201 | |
---|
202 | AIDA::ITuple* startPos = track->getTuple( 8 ); //start position |
---|
203 | startPos->fill(0, vtx.x()/cm); |
---|
204 | startPos->fill(1, vtx.y()/cm); |
---|
205 | startPos->fill(2, vtx.z()/cm); |
---|
206 | startPos->addRow(); |
---|
207 | |
---|
208 | AIDA::ITuple* stopPos = track->getTuple( 9 ); //stop position |
---|
209 | stopPos->fill(0, vtx.x()/cm); |
---|
210 | stopPos->fill(1, vtx.y()/cm); |
---|
211 | stopPos->fill(2, vtx.z()/cm); |
---|
212 | stopPos->addRow(); |
---|
213 | |
---|
214 | startVol = stopVol = -1; |
---|
215 | track->fill(10, startVol); //Starting Volume |
---|
216 | track->fill(11, stopVol); //Stopping Volume |
---|
217 | |
---|
218 | //add the Beam track to tuple |
---|
219 | track->addRow(); |
---|
220 | |
---|
221 | G4cout << "----> Tk{Beam}: " |
---|
222 | << " pId " << pId |
---|
223 | << " parent " << parent |
---|
224 | << " creation time " << timeStart |
---|
225 | << " Volumes " << startVol << " " << stopVol << "\n" |
---|
226 | << " dx,dy,dz " << dx << " " << dy << " " << dz << "\n" |
---|
227 | << " ETot " << ETot |
---|
228 | << " Mass " << mass |
---|
229 | << " pTot " << pTot |
---|
230 | << " px,py,pz " << px << " " << py << " " << pz << "\n" |
---|
231 | << G4endl; |
---|
232 | |
---|
233 | // -------------------------- |
---|
234 | // Loop over Trajectories |
---|
235 | // -------------------------- |
---|
236 | |
---|
237 | G4TrajectoryContainer* trajectoryContainer = evt->GetTrajectoryContainer(); |
---|
238 | |
---|
239 | G4int n_trajectories = 0; |
---|
240 | if (trajectoryContainer) n_trajectories = trajectoryContainer->entries(); |
---|
241 | |
---|
242 | G4int trjIndex; //this is the trajectory Index |
---|
243 | G4int parentIndex; //this the parent trajectory index |
---|
244 | |
---|
245 | G4ThreeVector mom; //px,py,pz |
---|
246 | G4double mommag; //p = sqrt(px*px + ...) |
---|
247 | G4ThreeVector stop; //stopping point |
---|
248 | G4ThreeVector start; //starting point |
---|
249 | G4String stopVolumeName; //Volume name of the stopping point for retreive |
---|
250 | |
---|
251 | eventTuple->fill(5, n_trajectories); //nPart |
---|
252 | G4cout << "Final # of tracks : " << n_trajectories << G4endl; |
---|
253 | |
---|
254 | |
---|
255 | for (G4int i=0; i < n_trajectories; i++) { |
---|
256 | Trajectory* trj = (Trajectory*)((*(evt->GetTrajectoryContainer()))[i]); |
---|
257 | |
---|
258 | pId = trj->GetPDGEncoding(); //pId |
---|
259 | trjIndex = trj->GetTrackID(); |
---|
260 | |
---|
261 | if ( ! trj->GetSaveFlag() ) continue; //only consider trajectories which are flaged (see TrackingAction, SteppingAction) |
---|
262 | |
---|
263 | // initial point of the trajectory |
---|
264 | G4TrajectoryPoint* aa = (G4TrajectoryPoint*)trj->GetPoint(0) ; |
---|
265 | |
---|
266 | track->fill(0, pId); //pId |
---|
267 | |
---|
268 | parentIndex = trj->GetParentID(); |
---|
269 | track->fill(1,parentIndex); //parent |
---|
270 | |
---|
271 | timeStart = trj->GetGlobalTime()/ns; |
---|
272 | track->fill(2, timeStart); //creation time |
---|
273 | |
---|
274 | mom = trj->GetInitialMomentum(); |
---|
275 | mommag = mom.mag(); |
---|
276 | direction = track->getTuple( 3 ); //direction |
---|
277 | dx = mom.x()/mommag; |
---|
278 | dy = mom.y()/mommag; |
---|
279 | dz = mom.z()/mommag; |
---|
280 | direction->fill(0, dx); |
---|
281 | direction->fill(1, dy); |
---|
282 | direction->fill(2, dz); |
---|
283 | direction->addRow(); |
---|
284 | |
---|
285 | mass = trj->GetParticleDefinition()->GetPDGMass(); |
---|
286 | track->fill(4, mass); //mass |
---|
287 | |
---|
288 | pTot = mommag; |
---|
289 | track->fill(5, pTot); //p Total |
---|
290 | |
---|
291 | ETot = sqrt(mom.mag2() + mass*mass); |
---|
292 | track->fill(6, ETot); //Total energy |
---|
293 | |
---|
294 | momentum = track->getTuple( 7 ); //3-momentum |
---|
295 | px = mom.x(); |
---|
296 | py = mom.y(); |
---|
297 | pz = mom.z(); |
---|
298 | momentum->fill(0, px); |
---|
299 | momentum->fill(1, py); |
---|
300 | momentum->fill(2, pz); |
---|
301 | momentum->addRow(); |
---|
302 | |
---|
303 | start = aa->GetPosition(); |
---|
304 | startPos = track->getTuple( 8 ); //start position |
---|
305 | startPos->fill(0, start.x()/cm); |
---|
306 | startPos->fill(1, start.y()/cm); |
---|
307 | startPos->fill(2, start.z()/cm); |
---|
308 | startPos->addRow(); |
---|
309 | |
---|
310 | stop = trj->GetStoppingPoint(); |
---|
311 | stopPos = track->getTuple( 9 ); //stop position |
---|
312 | stopPos->fill(0, stop.x()/cm); |
---|
313 | stopPos->fill(1, stop.y()/cm); |
---|
314 | stopPos->fill(2, stop.z()/cm); |
---|
315 | stopPos->addRow(); |
---|
316 | |
---|
317 | startVol = EventFindStartingVolume(start); |
---|
318 | track->fill(10, startVol); //Starting Volume |
---|
319 | |
---|
320 | |
---|
321 | stopVolumeName = trj->GetStoppingVolume()->GetName(); |
---|
322 | stopVol = EventFindStoppingVolume(stopVolumeName); |
---|
323 | track->fill(11, stopVol); //Stopping Volume |
---|
324 | |
---|
325 | //add the new track to tuple |
---|
326 | track->addRow(); |
---|
327 | |
---|
328 | G4cout << "----> Tk{"<<i<<"}: " |
---|
329 | << " pId " << pId |
---|
330 | << " parent " << parent |
---|
331 | << " creation time " << timeStart |
---|
332 | << " Volumes " << startVol << " " << stopVol << "\n" |
---|
333 | << " Start Pos (" << start.x()/cm << "," << start.y() << "," << start.z() << ")\n" |
---|
334 | << " Stop Pos (" << stop.x()/cm << "," << stop.y() << "," << stop.z() << ")\n" |
---|
335 | << " dx,dy,dz " << dx << " " << dy << " " << dz << "\n" |
---|
336 | << " m " << mass |
---|
337 | << " ETot " << ETot |
---|
338 | << " pTot " << pTot |
---|
339 | << " px,py,pz " << px << " " << py << " " << pz << "\n" |
---|
340 | << G4endl; |
---|
341 | |
---|
342 | }//end of loop on trajectories |
---|
343 | |
---|
344 | |
---|
345 | |
---|
346 | // -------------------- |
---|
347 | // Get Hit Collection |
---|
348 | // -------------------- |
---|
349 | //JEC 22/2/07 not yet available |
---|
350 | |
---|
351 | // G4SDManager* SDman = G4SDManager::GetSDMpointer(); //JEC FIXME: use data member |
---|
352 | |
---|
353 | // // Get Hit collection of this event |
---|
354 | // G4HCofThisEvent* HCE = evt->GetHCofThisEvent(); |
---|
355 | // TrappVolHitsCollection* TrappVolHC = 0; |
---|
356 | |
---|
357 | // AIDA::ITuple* hit = eventTuple->getTuple(6); //hit |
---|
358 | // G4int nHits=0; |
---|
359 | |
---|
360 | // if (HCE) { |
---|
361 | // G4String name = "TrappingVolume"; //see TrappingVolume.cxx Ctor |
---|
362 | // G4int collectionID = SDman->GetCollectionID(name); |
---|
363 | // TrappVolHC = (TrappVolHitsCollection*)HCE->GetHC(collectionID); |
---|
364 | // } |
---|
365 | |
---|
366 | // if (TrappVolHC) { |
---|
367 | |
---|
368 | // nHits = TrappVolHC->entries(); |
---|
369 | // eventTuple->fill(5, nHits); //nHits |
---|
370 | |
---|
371 | // TrappVolHit* aHit; |
---|
372 | |
---|
373 | // for (G4int i=0; i<nHits ;i++) { |
---|
374 | // aHit = (*TrappVolHC)[i]; |
---|
375 | |
---|
376 | // //dump |
---|
377 | // aHit->Print(); |
---|
378 | |
---|
379 | // hit->fill(0,aHit->GetTime()); //time of hit creation |
---|
380 | // hit->fill(1,aHit->GetEtot()); //Total energy of the particle creating the hit |
---|
381 | // hit->fill(2,aHit->GetCoordinates().x()); //position of the hit |
---|
382 | // hit->fill(3,aHit->GetCoordinates().y()); |
---|
383 | // hit->fill(4,aHit->GetCoordinates().z()); |
---|
384 | // hit->fill(5,aHit->GetTrkId()); //Id of the particle |
---|
385 | // hit->fill(6,aHit->GetParentId()); //Id of the parent of the particle |
---|
386 | // hit->fill(7,aHit->GetPDGEncoding()); //PDG code of the particle |
---|
387 | // hit->fill(8,aHit->GetEdep()); //Energy deposit |
---|
388 | |
---|
389 | // hit->addRow(); |
---|
390 | // }//loop on hits |
---|
391 | // }//Hit container |
---|
392 | |
---|
393 | //Save the Event |
---|
394 | eventTuple->addRow(); |
---|
395 | |
---|
396 | }//EndOfEventAction |
---|
397 | |
---|
398 | //-------------------------------------------------------------------------------------------------- |
---|
399 | |
---|
400 | G4int ELYSE::EventAction::EventFindStartingVolume(G4ThreeVector vtx) { |
---|
401 | // Get volume of starting point (see GEANT4 FAQ) |
---|
402 | |
---|
403 | G4Navigator* tmpNavigator = |
---|
404 | G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); |
---|
405 | |
---|
406 | G4VPhysicalVolume* tmpVolume = tmpNavigator->LocateGlobalPointAndSetup(vtx); |
---|
407 | G4String vtxVolumeName = tmpVolume->GetName(); |
---|
408 | |
---|
409 | return FindVolume(vtxVolumeName); |
---|
410 | |
---|
411 | }//EventFindStartingVolume |
---|
412 | |
---|
413 | //-------------------------------------------------------------------------------------------------- |
---|
414 | |
---|
415 | G4int ELYSE::EventAction::EventFindStoppingVolume(G4String stopVolumeName) { |
---|
416 | return FindVolume(stopVolumeName); |
---|
417 | }//EventFindStoppingVolume |
---|
418 | |
---|
419 | |
---|
420 | //-------------------------------------------------------------------------------------------------- |
---|
421 | G4int ELYSE::EventAction::FindVolume(G4String volumeName) { |
---|
422 | G4int volumeCode = -1; |
---|
423 | |
---|
424 | G4cout << "(FindVolume): VolumeName: <" |
---|
425 | << volumeName |
---|
426 | << ">" |
---|
427 | << G4endl; |
---|
428 | |
---|
429 | if ( volumeName == "expHall" ) { |
---|
430 | volumeCode = 0; |
---|
431 | } else if ( volumeName == "SensiVol" ) { |
---|
432 | volumeCode = 1; |
---|
433 | } else if (volumeName == "TyvekSheet") { |
---|
434 | volumeCode = 2; |
---|
435 | } else if (volumeName == "waterTank") { |
---|
436 | volumeCode = 3; |
---|
437 | } |
---|
438 | |
---|
439 | return volumeCode; |
---|
440 | }//FindVolume |
---|