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