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

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

update
File size: 9.3 KB
Line 
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.