Context Navigation

EventAction.cxx @ 294

Last change on this file since 294 was 294, checked in by campagne, 17 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
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 = ETotETot - massmass;
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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: