source: trunk/examples/extended/analysis/A01/src/A01EventAction.cc@ 812

Last change on this file since 812 was 807, checked in by garnier, 17 years ago

update
File size: 9.3 KB
RevLine 
[807]1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id: A01EventAction.cc,v 1.10 2007/05/17 09:55:14 duns Exp $
27// --------------------------------------------------------------
28//
29
30#include "A01EventAction.hh"
31#include "A01EventActionMessenger.hh"
32#ifdef G4ANALYSIS_USE
33#include "A01AnalysisManager.hh"
34#endif // G4ANALYSIS_USE
35
36#include "G4Event.hh"
37#include "G4EventManager.hh"
38#include "G4HCofThisEvent.hh"
39#include "G4VHitsCollection.hh"
40#include "G4TrajectoryContainer.hh"
41#include "G4Trajectory.hh"
42#include "G4VVisManager.hh"
43#include "G4SDManager.hh"
44#include "G4UImanager.hh"
45#include "G4ios.hh"
46
47#include "A01HodoscopeHit.hh"
48#include "A01DriftChamberHit.hh"
49#include "A01EmCalorimeterHit.hh"
50
51
52//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
53#include "A01HadCalorimeterHit.hh"
54//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
55
56A01EventAction::A01EventAction()
57{
58 G4String colName;
59 G4SDManager* SDman = G4SDManager::GetSDMpointer();
60 HHC1ID = SDman->GetCollectionID(colName="hodoscope1/hodoscopeColl");
61 HHC2ID = SDman->GetCollectionID(colName="hodoscope2/hodoscopeColl");
62 DHC1ID = SDman->GetCollectionID(colName="chamber1/driftChamberColl");
63 DHC2ID = SDman->GetCollectionID(colName="chamber2/driftChamberColl");
64 ECHCID = SDman->GetCollectionID(colName="EMcalorimeter/EMcalorimeterColl");
65//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
66 HCHCID = SDman->GetCollectionID(colName="HadCalorimeter/HadCalorimeterColl");
67//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
68 verboseLevel = 1;
69 messenger = new A01EventActionMessenger(this);
70
71#ifdef G4ANALYSIS_USE
72 plotter = 0;
73 tuple = 0;
74 dc1Hits = dc2Hits = 0;
75 dc1XY = dc2XY = evstof = 0;
76
77 // Do some analysis
78
79 A01AnalysisManager* analysisManager = A01AnalysisManager::getInstance();
80 IHistogramFactory* hFactory = analysisManager->getHistogramFactory();
81
82 if (hFactory)
83 {
84 // Create some histograms
85 dc1Hits = hFactory->createHistogram1D("Drift Chamber 1 # Hits",50,0,50);
86 dc2Hits = hFactory->createHistogram1D("Drift Chamber 2 # Hits",50,0,50);
87
88 // Create some clouds (Scatter Plots)
89 dc1XY = hFactory->createCloud2D("Drift Chamber 1 X vs Y");
90 dc2XY = hFactory->createCloud2D("Drift Chamber 2 X vs Y");
91 evstof = hFactory->createCloud2D("EDep vs Time-of-flight");
92
93 plotter = analysisManager->getPlotter();
94 if (plotter)
95 {
96 plotter->createRegions(3,2);
97 plotter->region(0)->plot(*dc1Hits);
98 plotter->region(1)->plot(*dc2Hits);
99 plotter->region(2)->plot(*dc1XY);
100 plotter->region(3)->plot(*dc2XY);
101 plotter->region(4)->plot(*evstof);
102 plotter->show();
103 }
104 }
105
106 // Create a Tuple
107
108 ITupleFactory* tFactory = analysisManager->getTupleFactory();
109 if (tFactory)
110 {
111 tuple = tFactory->create("MyTuple","MyTuple","int dc1Hits, dc2Hits, double ECEnergy, HCEnergy, time1, time2","");
112 }
113#endif // G4ANALYSIS_USE
114}
115
116A01EventAction::~A01EventAction()
117{
118#ifdef G4ANALYSIS_USE
119 A01AnalysisManager::dispose();
120#endif // G4ANALYSIS_USE
121 delete messenger;
122}
123
124void A01EventAction::BeginOfEventAction(const G4Event*)
125{
126}
127
128void A01EventAction::EndOfEventAction(const G4Event* evt)
129{
130 G4HCofThisEvent * HCE = evt->GetHCofThisEvent();
131 A01HodoscopeHitsCollection* HHC1 = 0;
132 A01HodoscopeHitsCollection* HHC2 = 0;
133 A01DriftChamberHitsCollection* DHC1 = 0;
134 A01DriftChamberHitsCollection* DHC2 = 0;
135 A01EmCalorimeterHitsCollection* ECHC = 0;
136 A01HadCalorimeterHitsCollection* HCHC = 0;
137 if(HCE)
138 {
139 HHC1 = (A01HodoscopeHitsCollection*)(HCE->GetHC(HHC1ID));
140 HHC2 = (A01HodoscopeHitsCollection*)(HCE->GetHC(HHC2ID));
141 DHC1 = (A01DriftChamberHitsCollection*)(HCE->GetHC(DHC1ID));
142 DHC2 = (A01DriftChamberHitsCollection*)(HCE->GetHC(DHC2ID));
143 ECHC = (A01EmCalorimeterHitsCollection*)(HCE->GetHC(ECHCID));
144 HCHC = (A01HadCalorimeterHitsCollection*)(HCE->GetHC(HCHCID));
145 }
146
147#ifdef G4ANALYSIS_USE
148 // Fill some histograms
149
150 if (DHC1 && dc1Hits)
151 {
152 int n_hit = DHC1->entries();
153 dc1Hits->fill(n_hit);
154 for(int i1=0;i1<n_hit;i1++)
155 {
156 A01DriftChamberHit* aHit = (*DHC1)[i1];
157 G4ThreeVector localPos = aHit->GetLocalPos();
158 if (dc1XY) dc1XY->fill(localPos.x(), localPos.y());
159 }
160 }
161 if (DHC2 && dc2Hits)
162 {
163 int n_hit = DHC2->entries();
164 dc2Hits->fill(n_hit);
165 for(int i1=0;i1<n_hit;i1++)
166 {
167 A01DriftChamberHit* aHit = (*DHC2)[i1];
168 G4ThreeVector localPos = aHit->GetLocalPos();
169 if (dc2XY) dc2XY->fill(localPos.x(), localPos.y());
170 }
171 }
172
173 // Fill the tuple
174
175 if (tuple)
176 {
177 if (DHC1) tuple->fill(0,DHC1->entries());
178 if (DHC2) tuple->fill(1,DHC2->entries());
179 if(ECHC)
180 {
181 int iHit = 0;
182 double totalE = 0.;
183 for(int i1=0;i1<80;i1++)
184 {
185 A01EmCalorimeterHit* aHit = (*ECHC)[i1];
186 double eDep = aHit->GetEdep();
187 if(eDep>0.)
188 {
189 iHit++;
190 totalE += eDep;
191 }
192 }
193 tuple->fill(2,totalE);
194
195 if (HHC1 && HHC2 && HHC1->entries()==1 && HHC2->entries()==1)
196 {
197 double tof = (*HHC2)[0]->GetTime() - (*HHC1)[0]->GetTime();
198 if (evstof) evstof->fill(tof,totalE);
199 }
200 }
201 if(HCHC)
202 {
203 int iHit = 0;
204 double totalE = 0.;
205 for(int i1=0;i1<20;i1++)
206 {
207 A01HadCalorimeterHit* aHit = (*HCHC)[i1];
208 double eDep = aHit->GetEdep();
209 if(eDep>0.)
210 {
211 iHit++;
212 totalE += eDep;
213 }
214 }
215 tuple->fill(3,totalE);
216 }
217 if (HHC1 && HHC1->entries()==1) tuple->fill(4,(*HHC1)[0]->GetTime());
218 if (HHC2 && HHC2->entries()==1) tuple->fill(5,(*HHC2)[0]->GetTime());
219 tuple->addRow();
220 }
221 if (plotter) plotter->refresh();
222#endif // G4ANALYSIS_USE
223
224
225 // Diagnostics
226
227 if (verboseLevel==0 || evt->GetEventID() % verboseLevel != 0) return;
228
229 G4PrimaryParticle* primary = evt->GetPrimaryVertex(0)->GetPrimary(0);
230 G4cout << G4endl
231 << ">>> Event " << evt->GetEventID() << " >>> Simulation truth : "
232 << primary->GetG4code()->GetParticleName()
233 << " " << primary->GetMomentum() << G4endl;
234
235 if(HHC1)
236 {
237 int n_hit = HHC1->entries();
238 G4cout << "Hodoscope 1 has " << n_hit << " hits." << G4endl;
239 for(int i1=0;i1<n_hit;i1++)
240 {
241 A01HodoscopeHit* aHit = (*HHC1)[i1];
242 aHit->Print();
243 }
244 }
245 if(HHC2)
246 {
247 int n_hit = HHC2->entries();
248 G4cout << "Hodoscope 2 has " << n_hit << " hits." << G4endl;
249 for(int i1=0;i1<n_hit;i1++)
250 {
251 A01HodoscopeHit* aHit = (*HHC2)[i1];
252 aHit->Print();
253 }
254 }
255 if(DHC1)
256 {
257 int n_hit = DHC1->entries();
258 G4cout << "Drift Chamber 1 has " << n_hit << " hits." << G4endl;
259 for(int i2=0;i2<5;i2++)
260 {
261 for(int i1=0;i1<n_hit;i1++)
262 {
263 A01DriftChamberHit* aHit = (*DHC1)[i1];
264 if(aHit->GetLayerID()==i2) aHit->Print();
265 }
266 }
267 }
268 if(DHC2)
269 {
270 int n_hit = DHC2->entries();
271 G4cout << "Drift Chamber 2 has " << n_hit << " hits." << G4endl;
272 for(int i2=0;i2<5;i2++)
273 {
274 for(int i1=0;i1<n_hit;i1++)
275 {
276 A01DriftChamberHit* aHit = (*DHC2)[i1];
277 if(aHit->GetLayerID()==i2) aHit->Print();
278 }
279 }
280 }
281 if(ECHC)
282 {
283 int iHit = 0;
284 double totalE = 0.;
285 for(int i1=0;i1<80;i1++)
286 {
287 A01EmCalorimeterHit* aHit = (*ECHC)[i1];
288 double eDep = aHit->GetEdep();
289 if(eDep>0.)
290 {
291 iHit++;
292 totalE += eDep;
293 }
294 }
295 G4cout << "EM Calorimeter has " << iHit << " hits. Total Edep is "
296 << totalE/MeV << " (MeV)" << G4endl;
297 }
298 if(HCHC)
299 {
300 int iHit = 0;
301 double totalE = 0.;
302 for(int i1=0;i1<20;i1++)
303 {
304 A01HadCalorimeterHit* aHit = (*HCHC)[i1];
305 double eDep = aHit->GetEdep();
306 if(eDep>0.)
307 {
308 iHit++;
309 totalE += eDep;
310 }
311 }
312 G4cout << "Hadron Calorimeter has " << iHit << " hits. Total Edep is "
313 << totalE/MeV << " (MeV)" << G4endl;
314 }
315}
316
317
Note: See TracBrowser for help on using the repository browser.