source: ELYSE/HEAD/source/EventAction.cxx@ 626

Last change on this file since 626 was 294, checked in by campagne, 19 years ago

Unuseful Sensitive Detector for the moment (JEC)
File size: 13.5 KB
Line 
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
42ELYSE::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
79ELYSE::EventAction::~EventAction(){
80}//Dtor
81
82//-------------------------------------------------------------------------------------------
83
84void 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
93void 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
400G4int 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
415G4int ELYSE::EventAction::EventFindStoppingVolume(G4String stopVolumeName) {
416 return FindVolume(stopVolumeName);
417}//EventFindStoppingVolume
418
419
420//--------------------------------------------------------------------------------------------------
421G4int 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
Note: See TracBrowser for help on using the repository browser.