source: trunk/examples/advanced/gammaray_telescope/src/GammaRayTelEventAction.cc @ 1282

Last change on this file since 1282 was 1230, checked in by garnier, 14 years ago

update to geant4.9.3
File size: 8.5 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//
27// $Id: GammaRayTelEventAction.cc,v 1.19 2006/06/29 15:56:39 gunter Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29// ------------------------------------------------------------
30//      GEANT 4 class implementation file
31//      CERN Geneva Switzerland
32//
33//
34//      ------------ GammaRayTelEventAction  ------
35//           by  R.Giannitrapani, F.Longo & G.Santin (13 nov 2000)
36//
37// - inclusion of Digits by F.Longo & R.Giannitrapani (24 oct 2001)
38//
39// - Modification of analysis management by G.Santin (18 Nov 2001)
40//
41// ************************************************************
42
43#include "GammaRayTelEventAction.hh"
44#include "GammaRayTelTrackerHit.hh"
45#include "GammaRayTelAnticoincidenceHit.hh"
46#include "GammaRayTelCalorimeterHit.hh"
47
48#ifdef G4ANALYSIS_USE
49#include "GammaRayTelAnalysis.hh"
50#endif
51
52#include "G4Event.hh"
53#include "G4EventManager.hh"
54#include "G4HCofThisEvent.hh"
55#include "G4VHitsCollection.hh"
56#include "G4TrajectoryContainer.hh"
57#include "G4Trajectory.hh"
58#include "G4VVisManager.hh"
59#include "G4SDManager.hh"
60#include "G4UImanager.hh"
61#include "G4ios.hh"
62#include "G4UnitsTable.hh"
63#include "Randomize.hh"
64
65#include "GammaRayTelDigi.hh"
66#include "GammaRayTelDigitizer.hh"
67#include "G4DigiManager.hh"
68
69// This file is a global variable in which we store energy deposition per hit
70// and other relevant information
71
72#ifdef G4STORE_DATA
73extern std::ofstream outFile;
74#endif
75
76//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77
78GammaRayTelEventAction::GammaRayTelEventAction()
79  :trackerCollID(-1),calorimeterCollID(-1),               
80  anticoincidenceCollID(-1), drawFlag("all")
81{ 
82  G4DigiManager * fDM = G4DigiManager::GetDMpointer();
83  GammaRayTelDigitizer * myDM = new GammaRayTelDigitizer( "GammaRayTelDigitizer" );
84  fDM->AddNewModule(myDM);
85}
86
87//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88
89GammaRayTelEventAction::~GammaRayTelEventAction()
90{
91}
92
93//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94
95void GammaRayTelEventAction::BeginOfEventAction(const G4Event* evt)
96{
97
98  G4int evtNb = evt->GetEventID();
99  G4cout << "Event: " << evtNb << G4endl;
100  G4SDManager * SDman = G4SDManager::GetSDMpointer(); 
101
102  if (trackerCollID==-1) {
103    trackerCollID = SDman->GetCollectionID("TrackerCollection");
104  }
105  if(anticoincidenceCollID==-1) {
106    anticoincidenceCollID =
107      SDman->GetCollectionID("AnticoincidenceCollection");
108  }
109  if(calorimeterCollID==-1) {
110    calorimeterCollID =
111      SDman->GetCollectionID("CalorimeterCollection");
112  }
113}
114
115//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116
117void GammaRayTelEventAction::EndOfEventAction(const G4Event* evt)
118{
119  G4int event_id = evt->GetEventID();
120
121  G4TrajectoryContainer * trajectoryContainer = evt->GetTrajectoryContainer();
122  G4int n_trajectories = 0;
123  if (trajectoryContainer) n_trajectories = trajectoryContainer->entries();
124 
125 
126  G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
127  GammaRayTelTrackerHitsCollection* THC = 0;
128  GammaRayTelCalorimeterHitsCollection* CHC = 0;
129  GammaRayTelAnticoincidenceHitsCollection* AHC = 0;
130
131
132  G4DigiManager * fDM = G4DigiManager::GetDMpointer();
133
134  if (HCE)
135    {
136      THC = (GammaRayTelTrackerHitsCollection*)(HCE->GetHC(trackerCollID));
137      CHC = (GammaRayTelCalorimeterHitsCollection*)
138        (HCE->GetHC(calorimeterCollID));
139      AHC = (GammaRayTelAnticoincidenceHitsCollection*)
140        (HCE->GetHC(anticoincidenceCollID));
141     
142      if (THC)
143        {
144          int n_hit = THC->entries();
145          G4cout << "Number of tracker hits in this event =  " << n_hit << G4endl;
146          G4double ESil=0;
147          G4int NStrip, NPlane, IsX;
148         
149          // This is a cycle on all the tracker hits of this event
150         
151          for (int i=0;i<n_hit;i++) 
152            {
153              // Here we put the hit data in a an ASCII file for
154              // later analysis
155              ESil = (*THC)[i]->GetEdepSil();
156              NStrip = (*THC)[i]->GetNStrip();
157              NPlane = (*THC)[i]->GetNSilPlane();
158              IsX = (*THC)[i]->GetPlaneType();
159             
160#ifdef G4STORE_DATA
161              outFile << std::setw(7) << event_id << " " << 
162                ESil/keV << " " << NStrip << 
163                " " << NPlane << " " << IsX << " " <<
164                (*THC)[i]->GetPos().x()/mm <<" "<<
165                (*THC)[i]->GetPos().y()/mm <<" "<<
166                (*THC)[i]->GetPos().z()/mm <<" "<<
167                G4endl;
168#else     
169              G4cout << std::setw(7) << event_id << " " << 
170                ESil/keV << " " << NStrip << 
171                " " << NPlane << " " << IsX << " " <<
172                (*THC)[i]->GetPos().x()/mm <<" "<<
173                (*THC)[i]->GetPos().y()/mm <<" "<<
174                (*THC)[i]->GetPos().z()/mm <<" "<<
175                G4endl;
176#endif
177             
178#ifdef G4ANALYSIS_USE
179             
180              // Here we fill the histograms of the Analysis manager
181              GammaRayTelAnalysis* analysis = GammaRayTelAnalysis::getInstance();
182             
183              if(IsX) 
184                {
185                  if (analysis->GetHisto2DMode()=="position")
186                    analysis->InsertPositionXZ((*THC)[i]->GetPos().x()/mm,(*THC)[i]->GetPos().z()/mm);
187                  else
188                    analysis->InsertPositionXZ(NStrip, NPlane);               
189                  if (NPlane == 0) analysis->InsertEnergy(ESil/keV);
190                  analysis->InsertHits(NPlane);
191                } 
192              else 
193                {
194                  if (analysis->GetHisto2DMode()=="position")
195                    analysis->InsertPositionYZ((*THC)[i]->GetPos().y()/mm,(*THC)[i]->GetPos().z()/mm); 
196                  else 
197                    analysis->InsertPositionYZ(NStrip, NPlane);               
198                  if (NPlane == 0) analysis->InsertEnergy(ESil/keV);
199                  analysis->InsertHits(NPlane);
200                }
201             
202#ifdef G4ANALYSIS_USE
203              analysis->setNtuple( ESil/keV, NPlane, (*THC)[i]->GetPos().x()/mm,
204                                   (*THC)[i]->GetPos().y()/mm,
205                                   (*THC)[i]->GetPos().z()/mm);
206#endif
207             
208#endif
209         
210            }
211          // Here we call the analysis manager function for visualization
212#ifdef G4ANALYSIS_USE
213          GammaRayTelAnalysis* analysis = GammaRayTelAnalysis::getInstance();
214          analysis->EndOfEvent(n_hit);
215#endif
216         
217        }
218     
219      GammaRayTelDigitizer * myDM = 
220        (GammaRayTelDigitizer*)fDM->FindDigitizerModule( "GammaRayTelDigitizer" );
221      myDM->Digitize();
222     
223      G4int myDigiCollID = fDM->GetDigiCollectionID("DigitsCollection");
224     
225      // G4cout << "digi collecion" << myDigiCollID << G4endl;
226     
227      GammaRayTelDigitsCollection * DC = (GammaRayTelDigitsCollection*)fDM->GetDigiCollection( myDigiCollID );
228     
229      if(DC) {
230        //    G4cout << "Total Digits " << DC->entries() << G4endl;
231        G4int n_digi =  DC->entries();
232        G4int NStrip, NPlane, IsX;
233        for (G4int i=0;i<n_digi;i++) {
234          // Here we put the digi data in a an ASCII file for
235          // later analysis
236          NStrip = (*DC)[i]->GetStripNumber();
237          NPlane = (*DC)[i]->GetPlaneNumber();
238          IsX = (*DC)[i]->GetPlaneType();
239         
240          //      outFile << std::setw(7) << event_id << " " << NStrip <<
241          //    " " << NPlane << " " << IsX << " " << G4endl;   
242        }
243      }
244    }
245 
246  if(G4VVisManager::GetConcreteInstance()) {
247    for(G4int i=0; i<n_trajectories; i++) { 
248      G4Trajectory* trj = (G4Trajectory *)((*(evt->GetTrajectoryContainer()))[i]);
249      if (drawFlag == "all") trj->DrawTrajectory(50);
250      else if ((drawFlag == "charged")&&(trj->GetCharge() != 0.))
251        trj->DrawTrajectory(50); 
252    }
253  }
254}
255
256
257
258
259
260
261
262
263
264
265
266
267
Note: See TracBrowser for help on using the repository browser.