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