| Geant4 Cross Reference |
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 /// \file RE05/src/RE05StackingAction.cc
28 /// \brief Implementation of the RE05StackingAction class
29 //
30
31 #include "RE05StackingAction.hh"
32
33 #include "RE05StackingActionMessenger.hh"
34
35 #include "G4Event.hh"
36 #include "G4HCofThisEvent.hh"
37 #include "G4ParticleDefinition.hh"
38 #include "G4ParticleTypes.hh"
39 #include "G4RunManager.hh"
40 #include "G4SDManager.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4Track.hh"
43 #include "G4TrackStatus.hh"
44 #include "G4ios.hh"
45
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47
48 RE05StackingAction::RE05StackingAction()
49 : G4UserStackingAction(),
50 fTrkHits(0),
51 fMuonHits(0),
52 fMessenger(0),
53 fStage(0),
54 fReqMuon(2),
55 fReqIso(10),
56 fAngRoI(30.0 * deg)
57 {
58 fMessenger = new RE05StackingActionMessenger(this);
59 }
60
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62
63 RE05StackingAction::~RE05StackingAction()
64 {
65 delete fMessenger;
66 }
67
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69
70 G4ClassificationOfNewTrack RE05StackingAction::ClassifyNewTrack(const G4Track* aTrack)
71 {
72 G4ClassificationOfNewTrack classification = fWaiting;
73 switch (fStage) {
74 case 0: // Stage 0 : Primary muons only
75 if (aTrack->GetParentID() == 0) {
76 G4ParticleDefinition* particleType = aTrack->GetDefinition();
77 if ((particleType == G4MuonPlus::MuonPlusDefinition())
78 || (particleType == G4MuonMinus::MuonMinusDefinition()))
79 {
80 classification = fUrgent;
81 }
82 }
83 break;
84
85 case 1: // Stage 1 : Charged primaries only
86 // Suspended tracks will be sent to the waiting stack
87 if (aTrack->GetParentID() != 0) {
88 break;
89 }
90 if (aTrack->GetTrackStatus() == fSuspend) {
91 break;
92 }
93 if (aTrack->GetDefinition()->GetPDGCharge() == 0.) {
94 break;
95 }
96 classification = fUrgent;
97 break;
98
99 default: // Stage 2 : Accept all primaries
100 // Accept all secondaries in RoI
101 // Kill secondaries outside RoI
102 if (aTrack->GetParentID() == 0) {
103 classification = fUrgent;
104 break;
105 }
106 if ((fAngRoI < 0.) || InsideRoI(aTrack, fAngRoI)) {
107 classification = fUrgent;
108 break;
109 }
110 classification = fKill;
111 }
112 return classification;
113 }
114
115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
116
117 G4bool RE05StackingAction::InsideRoI(const G4Track* aTrack, G4double ang)
118 {
119 if (!fMuonHits) {
120 fMuonHits = (RE05MuonHitsCollection*)GetCollection("muonCollection");
121 }
122 if (!fMuonHits) {
123 G4cerr << "muonCollection NOT FOUND" << G4endl;
124 return true;
125 }
126
127 G4int nhits = fMuonHits->entries();
128
129 const G4ThreeVector trPos = aTrack->GetPosition();
130 for (G4int i = 0; i < nhits; i++) {
131 G4ThreeVector muHitPos = (*fMuonHits)[i]->GetPos();
132 G4double angl = muHitPos.angle(trPos);
133 if (angl < ang) {
134 return true;
135 }
136 }
137
138 return false;
139 }
140
141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
142
143 G4VHitsCollection* RE05StackingAction::GetCollection(G4String colName)
144 {
145 G4SDManager* SDMan = G4SDManager::GetSDMpointer();
146 G4RunManager* runMan = G4RunManager::GetRunManager();
147 int colID = SDMan->GetCollectionID(colName);
148 if (colID >= 0) {
149 const G4Event* currentEvent = runMan->GetCurrentEvent();
150 G4HCofThisEvent* HCE = currentEvent->GetHCofThisEvent();
151 return HCE->GetHC(colID);
152 }
153 return 0;
154 }
155
156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
157
158 void RE05StackingAction::NewStage()
159 {
160 fStage++;
161 G4int nhits;
162 if (fStage == 1) {
163 // Stage 0->1 : check if at least "fReqMuon" hits on muon chamber
164 // otherwise abort current event
165 if (!fMuonHits) {
166 fMuonHits = (RE05MuonHitsCollection*)GetCollection("muonCollection");
167 }
168 if (!fMuonHits) {
169 G4cerr << "muonCollection NOT FOUND" << G4endl;
170 return;
171 }
172 nhits = fMuonHits->entries();
173 //// G4cout << "Stage 0->1 : " << nhits << " hits found in the muon chamber."
174 //// << G4endl;
175 if (nhits < fReqMuon) {
176 stackManager->clear();
177 //// G4cout << "++++++++ event aborted" << G4endl;
178 return;
179 }
180 stackManager->ReClassify();
181 return;
182 }
183
184 else if (fStage == 2) {
185 // Stage 1->2 : check the isolation of muon tracks
186 // at least "fReqIsoMuon" isolated muons
187 // otherwise abort current event.
188 // Isolation requires "fReqIso" or less hits
189 // (including own hits) in the RoI region
190 // in the tracker layers.
191 nhits = fMuonHits->entries();
192 if (!fTrkHits) {
193 fTrkHits = (RE05TrackerHitsCollection*)GetCollection("trackerCollection");
194 }
195 if (!fTrkHits) {
196 G4cerr << "trackerCollection NOT FOUND" << G4endl;
197 return;
198 }
199 G4int nTrkhits = fTrkHits->entries();
200 G4int isoMuon = 0;
201 for (G4int j = 0; j < nhits; j++) {
202 G4ThreeVector hitPos = (*fMuonHits)[j]->GetPos();
203 G4int nhitIn = 0;
204 for (G4int jj = 0; (jj < nTrkhits) && (nhitIn <= fReqIso); jj++) {
205 G4ThreeVector trkhitPos = (*fTrkHits)[jj]->GetPos();
206 if (trkhitPos.angle(hitPos) < fAngRoI) nhitIn++;
207 }
208 if (nhitIn <= fReqIso) isoMuon++;
209 }
210 //// G4cout << "Stage 1->2 : " << isoMuon << " isolated muon found." << G4endl;
211 if (isoMuon < fReqIsoMuon) {
212 stackManager->clear();
213 //// G4cout << "++++++++ event aborted" << G4endl;
214 return;
215 }
216 stackManager->ReClassify();
217 return;
218 }
219
220 else {
221 // Other fStage change : just re-classify
222 stackManager->ReClassify();
223 }
224 }
225
226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
227
228 void RE05StackingAction::PrepareNewEvent()
229 {
230 fStage = 0;
231 fTrkHits = 0;
232 fMuonHits = 0;
233 }
234
235 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
236