Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
eventgenerator/HepMC/HepMCEx01/src/ExN04StackingAction.cc
Go to the documentation of this file.
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 /// \file eventgenerator/HepMC/HepMCEx01/src/ExN04StackingAction.cc
27 /// \brief Implementation of the ExN04StackingAction class
28 //
29 // $Id: ExN04StackingAction.cc 77801 2013-11-28 13:33:20Z gcosmo $
30 //
31 
32 #include "G4SDManager.hh"
33 #include "G4RunManager.hh"
34 #include "G4Event.hh"
35 #include "G4HCofThisEvent.hh"
36 #include "G4Track.hh"
37 #include "G4TrackStatus.hh"
38 #include "G4ParticleDefinition.hh"
39 #include "G4ParticleTypes.hh"
40 #include "ExN04StackingActionMessenger.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "ExN04StackingAction.hh"
43 
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46  : trkHits(0), muonHits(0), stage(0)
47 {
48  angRoI = 30.0*deg;
49  reqMuon = 2;
50  reqIso = 10;
51  theMessenger = new ExN04StackingActionMessenger(this);
52 }
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 { delete theMessenger; }
57 
58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 {
62  G4ClassificationOfNewTrack classification = fWaiting;
63  switch(stage)
64  {
65  case 0: // Stage 0 : Primary muons only
66  if(aTrack->GetParentID()==0)
67  {
68  G4ParticleDefinition * particleType = aTrack->GetDefinition();
69  if((particleType==G4MuonPlus::MuonPlusDefinition())
70  ||(particleType==G4MuonMinus::MuonMinusDefinition()))
71  { classification = fUrgent; }
72  }
73  break;
74 
75  case 1: // Stage 1 : Charged primaries only
76  // Suspended tracks will be sent to the waiting stack
77  if(aTrack->GetParentID()!=0) { break; }
78  if(aTrack->GetTrackStatus()==fSuspend) { break; }
79  if(aTrack->GetDefinition()->GetPDGCharge()==0.) { break; }
80  classification = fUrgent;
81  break;
82 
83  default: // Stage 2 : Accept all primaries
84  // Accept all secondaries in RoI
85  // Kill secondaries outside RoI
86  if(aTrack->GetParentID()==0)
87  {
88  classification = fUrgent;
89  break;
90  }
91  if((angRoI<0.)||InsideRoI(aTrack,angRoI))
92  {
93  classification = fUrgent;
94  break;
95  }
96  classification = fKill;
97  }
98  return classification;
99 }
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102 G4bool ExN04StackingAction::InsideRoI(const G4Track * aTrack,G4double ang)
103 {
104  if(!muonHits)
105  { muonHits = (ExN04MuonHitsCollection*)GetCollection("muonCollection"); }
106  if(!muonHits)
107  { G4cerr << "muonCollection NOT FOUND" << G4endl;
108  return true; }
109 
110  G4int nhits = muonHits->entries();
111 
112  const G4ThreeVector trPos = aTrack->GetPosition();
113  for(G4int i=0;i<nhits;i++)
114  {
115  G4ThreeVector muHitPos = (*muonHits)[i]->GetPos();
116  G4double angl = muHitPos.angle(trPos);
117  if(angl<ang) { return true; }
118  }
119 
120  return false;
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
124 G4VHitsCollection* ExN04StackingAction::GetCollection(G4String colName)
125 {
128  int colID = SDMan->GetCollectionID(colName);
129  if(colID>=0)
130  {
131  const G4Event* currentEvent = runMan->GetCurrentEvent();
132  G4HCofThisEvent* HCE = currentEvent->GetHCofThisEvent();
133  return HCE->GetHC(colID);
134  }
135  return 0;
136 }
137 
138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
140 {
141  stage++;
142  G4int nhits;
143  if(stage==1)
144  {
145  // Stage 0->1 : check if at least "reqMuon" hits on muon chamber
146  // otherwise abort current event
147  if(!muonHits)
148  { muonHits = (ExN04MuonHitsCollection*)GetCollection("muonCollection"); }
149  if(!muonHits)
150  { G4cerr << "muonCollection NOT FOUND" << G4endl;
151  return; }
152  nhits = muonHits->entries();
153  G4cout << "Stage 0->1 : " << nhits << " hits found in the muon chamber."
154  << G4endl;
155  if(nhits<reqMuon)
156  {
157  stackManager->clear();
158  G4cout << "++++++++ event aborted" << G4endl;
159  return;
160  }
162  return;
163  }
164 
165  else if(stage==2)
166  {
167  // Stage 1->2 : check the isolation of muon tracks
168  // at least "reqIsoMuon" isolated muons
169  // otherwise abort current event.
170  // Isolation requires "reqIso" or less hits
171  // (including own hits) in the RoI region
172  // in the tracker layers.
173  nhits = muonHits->entries();
174  if(!trkHits)
175  { trkHits =
176  (ExN04TrackerHitsCollection*)GetCollection("trackerCollection"); }
177  if(!trkHits)
178  { G4cerr << "trackerCollection NOT FOUND" << G4endl;
179  return; }
180  G4int nTrkhits = trkHits->entries();
181  G4int isoMuon = 0;
182  for(G4int j=0;j<nhits;j++)
183  {
184  G4ThreeVector hitPos = (*muonHits)[j]->GetPos();
185  G4int nhitIn = 0;
186  for(G4int jj=0;(jj<nTrkhits)&&(nhitIn<=reqIso);jj++)
187  {
188  G4ThreeVector trkhitPos = (*trkHits)[jj]->GetPos();
189  if(trkhitPos.angle(hitPos)<angRoI) nhitIn++;
190  }
191  if(nhitIn<=reqIso) isoMuon++;
192  }
193  G4cout << "Stage 1->2 : " << isoMuon << " isolated muon found." << G4endl;
194  if(isoMuon<reqIsoMuon)
195  {
196  stackManager->clear();
197  G4cout << "++++++++ event aborted" << G4endl;
198  return;
199  }
201  return;
202  }
203 
204  else
205  {
206  // Other stage change : just re-classify
208  }
209 }
210 
211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
213 {
214  stage = 0;
215  trkHits = 0;
216  muonHits = 0;
217 }
G4ParticleDefinition * GetDefinition() const
G4int GetParentID() const
static G4MuonPlus * MuonPlusDefinition()
Definition: G4MuonPlus.cc:94
virtual G4ClassificationOfNewTrack ClassifyNewTrack(const G4Track *aTrack)
G4int GetCollectionID(G4String colName)
Definition: G4SDManager.cc:131
double angle(const Hep3Vector &) const
const G4ThreeVector & GetPosition() const
G4TrackStatus GetTrackStatus() const
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:74
static G4SDManager * GetSDMpointer()
Definition: G4SDManager.cc:40
static G4MuonMinus * MuonMinusDefinition()
Definition: G4MuonMinus.cc:95
#define G4endl
Definition: G4ios.hh:61
G4HCofThisEvent * GetHCofThisEvent() const
Definition: G4Event.hh:174
G4StackManager * stackManager
double G4double
Definition: G4Types.hh:76
const G4Event * GetCurrentEvent() const
G4double GetPDGCharge() const
G4GLOB_DLL std::ostream G4cerr