| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: RE05StackingAction.cc 77724 2013-11-27 17:11:16Z gcosmo $
26 // 27 //
27 /// \file RE05/src/RE05StackingAction.cc 28 /// \file RE05/src/RE05StackingAction.cc
28 /// \brief Implementation of the RE05StackingA 29 /// \brief Implementation of the RE05StackingAction class
29 // 30 //
30 31
31 #include "RE05StackingAction.hh" 32 #include "RE05StackingAction.hh"
32 << 33 #include "G4SDManager.hh"
33 #include "RE05StackingActionMessenger.hh" << 34 #include "G4RunManager.hh"
34 <<
35 #include "G4Event.hh" 35 #include "G4Event.hh"
36 #include "G4HCofThisEvent.hh" 36 #include "G4HCofThisEvent.hh"
>> 37 #include "G4Track.hh"
>> 38 #include "G4TrackStatus.hh"
37 #include "G4ParticleDefinition.hh" 39 #include "G4ParticleDefinition.hh"
38 #include "G4ParticleTypes.hh" 40 #include "G4ParticleTypes.hh"
39 #include "G4RunManager.hh" << 41 #include "RE05StackingActionMessenger.hh"
40 #include "G4SDManager.hh" <<
41 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh"
42 #include "G4Track.hh" <<
43 #include "G4TrackStatus.hh" <<
44 #include "G4ios.hh" 43 #include "G4ios.hh"
45 44
46 //....oooOO0OOooo........oooOO0OOooo........oo <<
47 <<
48 RE05StackingAction::RE05StackingAction() 45 RE05StackingAction::RE05StackingAction()
49 : G4UserStackingAction(), << 46 : trkHits(0), muonHits(0), stage(0)
50 fTrkHits(0), << 47 {
51 fMuonHits(0), << 48 angRoI = 30.0*deg;
52 fMessenger(0), << 49 reqMuon = 2;
53 fStage(0), << 50 reqIso = 10;
54 fReqMuon(2), << 51 theMessenger = new RE05StackingActionMessenger(this);
55 fReqIso(10), <<
56 fAngRoI(30.0 * deg) <<
57 { <<
58 fMessenger = new RE05StackingActionMessenger <<
59 } 52 }
60 53
61 //....oooOO0OOooo........oooOO0OOooo........oo <<
62 <<
63 RE05StackingAction::~RE05StackingAction() 54 RE05StackingAction::~RE05StackingAction()
64 { << 55 { delete theMessenger; }
65 delete fMessenger; <<
66 } <<
67 <<
68 //....oooOO0OOooo........oooOO0OOooo........oo <<
69 56
70 G4ClassificationOfNewTrack RE05StackingAction: << 57 G4ClassificationOfNewTrack
>> 58 RE05StackingAction::ClassifyNewTrack(const G4Track * aTrack)
71 { 59 {
72 G4ClassificationOfNewTrack classification = 60 G4ClassificationOfNewTrack classification = fWaiting;
73 switch (fStage) { << 61 switch(stage)
74 case 0: // Stage 0 : Primary muons only << 62 {
75 if (aTrack->GetParentID() == 0) { << 63 case 0: // Stage 0 : Primary muons only
76 G4ParticleDefinition* particleType = a << 64 if(aTrack->GetParentID()==0)
77 if ((particleType == G4MuonPlus::MuonP << 65 {
78 || (particleType == G4MuonMinus::M << 66 G4ParticleDefinition * particleType = aTrack->GetDefinition();
79 { << 67 if((particleType==G4MuonPlus::MuonPlusDefinition())
80 classification = fUrgent; << 68 ||(particleType==G4MuonMinus::MuonMinusDefinition()))
81 } << 69 { classification = fUrgent; }
82 } << 70 }
>> 71 break;
>> 72
>> 73 case 1: // Stage 1 : Charged primaries only
>> 74 // Suspended tracks will be sent to the waiting stack
>> 75 if(aTrack->GetParentID()!=0) { break; }
>> 76 if(aTrack->GetTrackStatus()==fSuspend) { break; }
>> 77 if(aTrack->GetDefinition()->GetPDGCharge()==0.) { break; }
>> 78 classification = fUrgent;
>> 79 break;
>> 80
>> 81 default: // Stage 2 : Accept all primaries
>> 82 // Accept all secondaries in RoI
>> 83 // Kill secondaries outside RoI
>> 84 if(aTrack->GetParentID()==0)
>> 85 {
>> 86 classification = fUrgent;
83 break; 87 break;
84 << 88 }
85 case 1: // Stage 1 : Charged primaries on << 89 if((angRoI<0.)||InsideRoI(aTrack,angRoI))
86 // Suspended tracks wil << 90 {
87 if (aTrack->GetParentID() != 0) { <<
88 break; <<
89 } <<
90 if (aTrack->GetTrackStatus() == fSuspend <<
91 break; <<
92 } <<
93 if (aTrack->GetDefinition()->GetPDGCharg <<
94 break; <<
95 } <<
96 classification = fUrgent; 91 classification = fUrgent;
97 break; 92 break;
98 << 93 }
99 default: // Stage 2 : Accept all primarie << 94 classification = fKill;
100 // Accept all secondar <<
101 // Kill secondaries ou <<
102 if (aTrack->GetParentID() == 0) { <<
103 classification = fUrgent; <<
104 break; <<
105 } <<
106 if ((fAngRoI < 0.) || InsideRoI(aTrack, <<
107 classification = fUrgent; <<
108 break; <<
109 } <<
110 classification = fKill; <<
111 } 95 }
112 return classification; 96 return classification;
113 } 97 }
114 98
115 //....oooOO0OOooo........oooOO0OOooo........oo << 99 G4bool RE05StackingAction::InsideRoI(const G4Track * aTrack,G4double ang)
116 <<
117 G4bool RE05StackingAction::InsideRoI(const G4T <<
118 { 100 {
119 if (!fMuonHits) { << 101 if(!muonHits)
120 fMuonHits = (RE05MuonHitsCollection*)GetCo << 102 { muonHits = (RE05MuonHitsCollection*)GetCollection("muonCollection"); }
121 } << 103 if(!muonHits)
122 if (!fMuonHits) { << 104 { G4cerr << "muonCollection NOT FOUND" << G4endl;
123 G4cerr << "muonCollection NOT FOUND" << G4 << 105 return true; }
124 return true; <<
125 } <<
126 106
127 G4int nhits = fMuonHits->entries(); << 107 G4int nhits = muonHits->entries();
128 108
129 const G4ThreeVector trPos = aTrack->GetPosit 109 const G4ThreeVector trPos = aTrack->GetPosition();
130 for (G4int i = 0; i < nhits; i++) { << 110 for(G4int i=0;i<nhits;i++)
131 G4ThreeVector muHitPos = (*fMuonHits)[i]-> << 111 {
>> 112 G4ThreeVector muHitPos = (*muonHits)[i]->GetPos();
132 G4double angl = muHitPos.angle(trPos); 113 G4double angl = muHitPos.angle(trPos);
133 if (angl < ang) { << 114 if(angl<ang) { return true; }
134 return true; <<
135 } <<
136 } 115 }
137 116
138 return false; 117 return false;
139 } 118 }
140 119
141 //....oooOO0OOooo........oooOO0OOooo........oo <<
142 <<
143 G4VHitsCollection* RE05StackingAction::GetColl 120 G4VHitsCollection* RE05StackingAction::GetCollection(G4String colName)
144 { 121 {
145 G4SDManager* SDMan = G4SDManager::GetSDMpoin 122 G4SDManager* SDMan = G4SDManager::GetSDMpointer();
146 G4RunManager* runMan = G4RunManager::GetRunM 123 G4RunManager* runMan = G4RunManager::GetRunManager();
147 int colID = SDMan->GetCollectionID(colName); 124 int colID = SDMan->GetCollectionID(colName);
148 if (colID >= 0) { << 125 if(colID>=0)
>> 126 {
149 const G4Event* currentEvent = runMan->GetC 127 const G4Event* currentEvent = runMan->GetCurrentEvent();
150 G4HCofThisEvent* HCE = currentEvent->GetHC 128 G4HCofThisEvent* HCE = currentEvent->GetHCofThisEvent();
151 return HCE->GetHC(colID); 129 return HCE->GetHC(colID);
152 } 130 }
153 return 0; 131 return 0;
154 } 132 }
155 133
156 //....oooOO0OOooo........oooOO0OOooo........oo <<
157 <<
158 void RE05StackingAction::NewStage() 134 void RE05StackingAction::NewStage()
159 { 135 {
160 fStage++; << 136 stage++;
161 G4int nhits; 137 G4int nhits;
162 if (fStage == 1) { << 138 if(stage==1)
163 // Stage 0->1 : check if at least "fReqMuo << 139 {
164 // otherwise abort current ev << 140 // Stage 0->1 : check if at least "reqMuon" hits on muon chamber
165 if (!fMuonHits) { << 141 // otherwise abort current event
166 fMuonHits = (RE05MuonHitsCollection*)Get << 142 if(!muonHits)
167 } << 143 { muonHits = (RE05MuonHitsCollection*)GetCollection("muonCollection"); }
168 if (!fMuonHits) { << 144 if(!muonHits)
169 G4cerr << "muonCollection NOT FOUND" << << 145 { G4cerr << "muonCollection NOT FOUND" << G4endl;
170 return; << 146 return; }
171 } << 147 nhits = muonHits->entries();
172 nhits = fMuonHits->entries(); << 148 //// G4cout << "Stage 0->1 : " << nhits << " hits found in the muon chamber."
173 //// G4cout << "Stage 0->1 : " << nhits << 149 //// << G4endl;
174 //// << G4endl; << 150 if(nhits<reqMuon)
175 if (nhits < fReqMuon) { << 151 {
176 stackManager->clear(); 152 stackManager->clear();
177 //// G4cout << "++++++++ event abor << 153 //// G4cout << "++++++++ event aborted" << G4endl;
178 return; 154 return;
179 } 155 }
180 stackManager->ReClassify(); 156 stackManager->ReClassify();
181 return; 157 return;
182 } 158 }
183 159
184 else if (fStage == 2) { << 160 else if(stage==2)
185 // Stage 1->2 : check the isolation of muo << 161 {
186 // at least "fReqIsoMuon" iso << 162 // Stage 1->2 : check the isolation of muon tracks
187 // otherwise abort current ev << 163 // at least "reqIsoMuon" isolated muons
188 // Isolation requires "fReqIs << 164 // otherwise abort current event.
189 // (including own hits) in th << 165 // Isolation requires "reqIso" or less hits
190 // in the tracker layers. << 166 // (including own hits) in the RoI region
191 nhits = fMuonHits->entries(); << 167 // in the tracker layers.
192 if (!fTrkHits) { << 168 nhits = muonHits->entries();
193 fTrkHits = (RE05TrackerHitsCollection*)G << 169 if(!trkHits)
194 } << 170 { trkHits = (RE05TrackerHitsCollection*)GetCollection("trackerCollection"); }
195 if (!fTrkHits) { << 171 if(!trkHits)
196 G4cerr << "trackerCollection NOT FOUND" << 172 { G4cerr << "trackerCollection NOT FOUND" << G4endl;
197 return; << 173 return; }
198 } << 174 G4int nTrkhits = trkHits->entries();
199 G4int nTrkhits = fTrkHits->entries(); <<
200 G4int isoMuon = 0; 175 G4int isoMuon = 0;
201 for (G4int j = 0; j < nhits; j++) { << 176 for(G4int j=0;j<nhits;j++)
202 G4ThreeVector hitPos = (*fMuonHits)[j]-> << 177 {
>> 178 G4ThreeVector hitPos = (*muonHits)[j]->GetPos();
203 G4int nhitIn = 0; 179 G4int nhitIn = 0;
204 for (G4int jj = 0; (jj < nTrkhits) && (n << 180 for(G4int jj=0;(jj<nTrkhits)&&(nhitIn<=reqIso);jj++)
205 G4ThreeVector trkhitPos = (*fTrkHits)[ << 181 {
206 if (trkhitPos.angle(hitPos) < fAngRoI) << 182 G4ThreeVector trkhitPos = (*trkHits)[jj]->GetPos();
207 } << 183 if(trkhitPos.angle(hitPos)<angRoI) nhitIn++;
208 if (nhitIn <= fReqIso) isoMuon++; << 184 }
209 } << 185 if(nhitIn<=reqIso) isoMuon++;
210 //// G4cout << "Stage 1->2 : " << isoMu << 186 }
211 if (isoMuon < fReqIsoMuon) { << 187 //// G4cout << "Stage 1->2 : " << isoMuon << " isolated muon found." << G4endl;
>> 188 if(isoMuon<reqIsoMuon)
>> 189 {
212 stackManager->clear(); 190 stackManager->clear();
213 //// G4cout << "++++++++ event abor << 191 //// G4cout << "++++++++ event aborted" << G4endl;
214 return; 192 return;
215 } 193 }
216 stackManager->ReClassify(); 194 stackManager->ReClassify();
217 return; 195 return;
218 } 196 }
219 197
220 else { << 198 else
221 // Other fStage change : just re-classify << 199 {
>> 200 // Other stage change : just re-classify
222 stackManager->ReClassify(); 201 stackManager->ReClassify();
223 } 202 }
224 } 203 }
225 << 204
226 //....oooOO0OOooo........oooOO0OOooo........oo <<
227 <<
228 void RE05StackingAction::PrepareNewEvent() 205 void RE05StackingAction::PrepareNewEvent()
229 { << 206 {
230 fStage = 0; << 207 stage = 0;
231 fTrkHits = 0; << 208 trkHits = 0;
232 fMuonHits = 0; << 209 muonHits = 0;
233 } 210 }
234 211
235 //....oooOO0OOooo........oooOO0OOooo........oo << 212
236 213