| 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: CCaloSD.cc 27 // File: CCaloSD.cc
28 // Description: Stores hits of calorimetric ty 28 // Description: Stores hits of calorimetric type in appropriate container
29 ////////////////////////////////////////////// 29 ///////////////////////////////////////////////////////////////////////////////
30 30
31 #include "CCaloSD.hh" 31 #include "CCaloSD.hh"
32 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh"
33 #include "G4VProcess.hh" 33 #include "G4VProcess.hh"
34 #include "G4SDManager.hh" 34 #include "G4SDManager.hh"
35 #include "G4VTouchable.hh" 35 #include "G4VTouchable.hh"
36 #include "CCalVOrganization.hh" 36 #include "CCalVOrganization.hh"
37 #include "CCalSDList.hh" 37 #include "CCalSDList.hh"
38 38
39 //#define debug 39 //#define debug
40 //#define ddebug 40 //#define ddebug
41 41
42 CCaloSD::CCaloSD(G4String name, CCalVOrganizat 42 CCaloSD::CCaloSD(G4String name, CCalVOrganization* numberingScheme):
43 G4VSensitiveDetector(name), HCID(-1), SDname 43 G4VSensitiveDetector(name), HCID(-1), SDname(name), theHC(0),
44 CurrentHit(0), theTrack(0), CurrentPV(0), Pr 44 CurrentHit(0), theTrack(0), CurrentPV(0), PreviousPV(0), UnitID(0),
45 PreviousUnitID(0), PreStepPoint(0), PostStep 45 PreviousUnitID(0), PreStepPoint(0), PostStepPoint(0),
46 theDescription(numberingScheme) { 46 theDescription(numberingScheme) {
47 47
48 collectionName.insert(name); 48 collectionName.insert(name);
49 49
50 G4cout << "********************************* 50 G4cout << "*******************************************************" << G4endl;
51 G4cout << "* 51 G4cout << "* *" << G4endl;
52 G4cout << "* Constructing a CCaloSD with na 52 G4cout << "* Constructing a CCaloSD with name " << name << G4endl;
53 G4cout << "* 53 G4cout << "* *" << G4endl;
54 G4cout << "********************************* 54 G4cout << "*******************************************************" << G4endl;
55 55
56 CCalSDList::getInstance()->addCalo(name); 56 CCalSDList::getInstance()->addCalo(name);
57 } 57 }
58 58
59 CCaloSD::~CCaloSD() { 59 CCaloSD::~CCaloSD() {
60 delete theDescription; 60 delete theDescription;
61 } 61 }
62 62
63 void CCaloSD::Initialize(G4HCofThisEvent*HCE) 63 void CCaloSD::Initialize(G4HCofThisEvent*HCE) {
64 64
65 #ifdef debug 65 #ifdef debug
66 G4cout << "CCaloSD : Initialize called for " 66 G4cout << "CCaloSD : Initialize called for " << SDname << G4endl;
67 #endif 67 #endif
68 //This initialization is performed at the be 68 //This initialization is performed at the beginning of an event
69 //------------------------------------------ 69 //------------------------------------------------------------
70 70
71 theHC = new CCalG4HitCollection(SDname, coll 71 theHC = new CCalG4HitCollection(SDname, collectionName[0]);
72 if (HCID<0) { 72 if (HCID<0) {
73 HCID = G4SDManager::GetSDMpointer()->GetCo 73 HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
74 } 74 }
75 HCE->AddHitsCollection( HCID, theHC ); 75 HCE->AddHitsCollection( HCID, theHC );
76 76
77 TSID = -2; 77 TSID = -2;
78 PrimID = -2; 78 PrimID = -2;
79 /////PrimaryID = -99; <--- initialized by S 79 /////PrimaryID = -99; <--- initialized by StackingAction.
80 } 80 }
81 81
82 G4bool CCaloSD::ProcessHits(G4Step*aStep, G4To 82 G4bool CCaloSD::ProcessHits(G4Step*aStep, G4TouchableHistory*) {
83 83
84 if(aStep->GetTotalEnergyDeposit() > CLHEP::e 84 if(aStep->GetTotalEnergyDeposit() > CLHEP::eV) {
85 getStepInfo(aStep); 85 getStepInfo(aStep);
86 if (hitExists() == false && EdepositEM+Ede 86 if (hitExists() == false && EdepositEM+EdepositEHAD>0.f) {
87 createNewHit(); 87 createNewHit();
88 } 88 }
89 } 89 }
90 return true; 90 return true;
91 } 91 }
92 92
93 void CCaloSD::getStepInfo(const G4Step* aStep) 93 void CCaloSD::getStepInfo(const G4Step* aStep) {
94 94
95 PreStepPoint = aStep->GetPreStepPoint(); 95 PreStepPoint = aStep->GetPreStepPoint();
96 PostStepPoint= aStep->GetPostStepPoint(); 96 PostStepPoint= aStep->GetPostStepPoint();
97 HitPoint = PreStepPoint->GetPosition(); 97 HitPoint = PreStepPoint->GetPosition();
98 98
99 theTrack = aStep->GetTrack(); 99 theTrack = aStep->GetTrack();
100 CurrentPV= PreStepPoint->GetPhysicalVolume() 100 CurrentPV= PreStepPoint->GetPhysicalVolume();
101 101
102 G4String pname = theTrack->GetDefinition()-> 102 G4String pname = theTrack->GetDefinition()->GetParticleName();
103 G4double de = aStep->GetTotalEnergyDeposit() 103 G4double de = aStep->GetTotalEnergyDeposit();
104 //G4cout << "##### Process new step dE= " << 104 //G4cout << "##### Process new step dE= " << de
105 // << " " << pname << " inside " << 105 // << " " << pname << " inside " << GetName() << G4endl;
106 106
107 const G4VTouchable* touch = aStep->GetPreSte 107 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
108 G4double weight = 1.; 108 G4double weight = 1.;
109 if (touch->GetVolume(0)->GetName() == "Cryst 109 if (touch->GetVolume(0)->GetName() == "CrystalMatrixCrystal") {
110 weight = curve_LY(PreStepPoint); 110 weight = curve_LY(PreStepPoint);
111 } 111 }
112 112
113 if (pname == "e-" || pname == "e+" || pname 113 if (pname == "e-" || pname == "e+" || pname == "gamma" ){
114 EdepositEM = weight*de; 114 EdepositEM = weight*de;
115 EdepositEHAD = 0.f; 115 EdepositEHAD = 0.f;
116 } else { 116 } else {
117 EdepositEM = 0.f; 117 EdepositEM = 0.f;
118 EdepositEHAD = weight*de; 118 EdepositEHAD = weight*de;
119 } 119 }
120 120
121 TSlice = (PostStepPoint) ? PostStepPoint->Ge 121 TSlice = (PostStepPoint) ? PostStepPoint->GetGlobalTime()/nanosecond : 0.;
122 //G4cout << " W= " << weight << " T= " < 122 //G4cout << " W= " << weight << " T= " << TSlice << G4endl;
123 G4int it = (G4int)TSlice; 123 G4int it = (G4int)TSlice;
124 TSliceID = std::min(it, 10000); 124 TSliceID = std::min(it, 10000);
125 //G4cout << " tID= " << TSliceID << G4endl; 125 //G4cout << " tID= " << TSliceID << G4endl;
126 126
127 UnitID = (theDescription) ? theDescription-> 127 UnitID = (theDescription) ? theDescription->GetUnitID(aStep) : 0;
128 } 128 }
129 129
130 G4bool CCaloSD::hitExists() { 130 G4bool CCaloSD::hitExists() {
131 131
132 if ( CurrentPV==PreviousPV && PrimaryID == P 132 if ( CurrentPV==PreviousPV && PrimaryID == PrimID && TSliceID == TSID &&
133 UnitID==PreviousUnitID) { 133 UnitID==PreviousUnitID) {
134 updateHit(); 134 updateHit();
135 return true; 135 return true;
136 } 136 }
137 137
138 if (PrimaryID != PrimID) { ResetForNewPrimar 138 if (PrimaryID != PrimID) { ResetForNewPrimary(); }
139 139
140 // look in HC whether a hit with the same pr 140 // look in HC whether a hit with the same primID,UnitID,TSliceID
141 // exists already: 141 // exists already:
142 142
143 G4bool found = false; 143 G4bool found = false;
144 for (std::size_t j=0; j<theHC->entries(); ++ 144 for (std::size_t j=0; j<theHC->entries(); ++j) {
145 145
146 CCalG4Hit* aPreviousHit = (*theHC)[j]; 146 CCalG4Hit* aPreviousHit = (*theHC)[j];
147 if (aPreviousHit->getTrackID() == Primary 147 if (aPreviousHit->getTrackID() == PrimaryID &&
148 aPreviousHit->getTimeSliceID() == TSli 148 aPreviousHit->getTimeSliceID() == TSliceID &&
149 aPreviousHit->getUnitID()== UnitID 149 aPreviousHit->getUnitID()== UnitID ) {
150 CurrentHit = aPreviousHit; 150 CurrentHit = aPreviousHit;
151 found = true; 151 found = true;
152 break; 152 break;
153 } 153 }
154 } 154 }
155 155
156 if (found) { 156 if (found) {
157 updateHit(); 157 updateHit();
158 return true; 158 return true;
159 } else { 159 } else {
160 return false; 160 return false;
161 } 161 }
162 } 162 }
163 163
164 void CCaloSD::ResetForNewPrimary() { 164 void CCaloSD::ResetForNewPrimary() {
165 165
166 EntrancePoint = SetToLocal(HitPoint); 166 EntrancePoint = SetToLocal(HitPoint);
167 IncidentEnergy = PreStepPoint->GetKineticEne 167 IncidentEnergy = PreStepPoint->GetKineticEnergy();
168 168
169 } 169 }
170 170
171 void CCaloSD::StoreHit(CCalG4Hit* hit){ 171 void CCaloSD::StoreHit(CCalG4Hit* hit){
172 172
173 if (PrimID<0) return; 173 if (PrimID<0) return;
174 if (hit == 0 ) { 174 if (hit == 0 ) {
175 G4cout << "CCaloSD: hit to be stored is NU 175 G4cout << "CCaloSD: hit to be stored is NULL !!" <<G4endl;
176 return; 176 return;
177 } 177 }
178 178
179 theHC->insert( hit ); 179 theHC->insert( hit );
180 } 180 }
181 181
182 void CCaloSD::createNewHit() { 182 void CCaloSD::createNewHit() {
183 183
184 #ifdef debug 184 #ifdef debug
185 G4int currentCopyNo = -999; 185 G4int currentCopyNo = -999;
186 G4int motherCopyNo = -999; 186 G4int motherCopyNo = -999;
187 G4TouchableHistory* theTouchable 187 G4TouchableHistory* theTouchable
188 = (G4TouchableHistory*)( theTrack->GetTouc 188 = (G4TouchableHistory*)( theTrack->GetTouchable() );
189 if ( theTouchable ) { 189 if ( theTouchable ) {
190 currentCopyNo = theTouchable->GetReplicaNu 190 currentCopyNo = theTouchable->GetReplicaNumber( 0 );
191 if ( theTouchable->GetHistoryDepth() > 0 ) 191 if ( theTouchable->GetHistoryDepth() > 0 ) {
192 motherCopyNo = theTouchable->GetReplicaN 192 motherCopyNo = theTouchable->GetReplicaNumber( 1 );
193 } 193 }
194 } 194 }
195 G4cout << "CCaloSD createNewHit for" 195 G4cout << "CCaloSD createNewHit for"
196 << " PV " << CurrentPV->GetName() 196 << " PV " << CurrentPV->GetName()
197 << " PVid = " << currentCopyNo 197 << " PVid = " << currentCopyNo
198 << " MVid = " << motherCopyNo 198 << " MVid = " << motherCopyNo
199 << " Unit " << UnitID <<G4endl; 199 << " Unit " << UnitID <<G4endl;
200 G4cout << " primary " << PrimaryID 200 G4cout << " primary " << PrimaryID
201 << " time slice " << TSliceID 201 << " time slice " << TSliceID
202 << " For Track " << theTrack->GetTra 202 << " For Track " << theTrack->GetTrackID()
203 << " which is a " << theTrack->GetDe 203 << " which is a " << theTrack->GetDefinition()->GetParticleName();
204 204
205 if (theTrack->GetTrackID()==1) { 205 if (theTrack->GetTrackID()==1) {
206 G4cout << " of energy " << theTrack->G 206 G4cout << " of energy " << theTrack->GetTotalEnergy();
207 } else { 207 } else {
208 G4cout << " daughter of part. " << theTrac 208 G4cout << " daughter of part. " << theTrack->GetParentID();
209 } 209 }
210 210
211 G4cout << " and created by " ; 211 G4cout << " and created by " ;
212 if (theTrack->GetCreatorProcess()!=NULL) 212 if (theTrack->GetCreatorProcess()!=NULL)
213 G4cout << theTrack->GetCreatorProcess()->G 213 G4cout << theTrack->GetCreatorProcess()->GetProcessName() ;
214 else 214 else
215 G4cout << "NO process"; 215 G4cout << "NO process";
216 G4cout << G4endl; 216 G4cout << G4endl;
217 #endif 217 #endif
218 218
219 CurrentHit = new CCalG4Hit; 219 CurrentHit = new CCalG4Hit;
220 CurrentHit->setTrackID(PrimaryID); 220 CurrentHit->setTrackID(PrimaryID);
221 CurrentHit->setTimeSlice(TSlice); 221 CurrentHit->setTimeSlice(TSlice);
222 CurrentHit->setUnitID(UnitID); 222 CurrentHit->setUnitID(UnitID);
223 CurrentHit->setEntry(EntrancePoint); 223 CurrentHit->setEntry(EntrancePoint);
224 CurrentHit->setIncidentEnergy(IncidentEnergy 224 CurrentHit->setIncidentEnergy(IncidentEnergy);
225 updateHit(); 225 updateHit();
226 226
227 StoreHit(CurrentHit); 227 StoreHit(CurrentHit);
228 228
229 } 229 }
230 230
231 void CCaloSD::updateHit() { 231 void CCaloSD::updateHit() {
232 if (EdepositEM+EdepositEHAD != 0) { 232 if (EdepositEM+EdepositEHAD != 0) {
233 CurrentHit->addEnergyDeposit(EdepositEM,Ed 233 CurrentHit->addEnergyDeposit(EdepositEM,EdepositEHAD);
234 #ifdef debug 234 #ifdef debug
235 G4cout << "Energy deposit in Unit " << Uni 235 G4cout << "Energy deposit in Unit " << UnitID << " em " << EdepositEM/MeV
236 << " hadronic " << EdepositEHAD/MeV < 236 << " hadronic " << EdepositEHAD/MeV << " MeV" << G4endl;
237 #endif 237 #endif
238 } 238 }
239 239
240 // buffer for next steps: 240 // buffer for next steps:
241 TSID = TSliceID; 241 TSID = TSliceID;
242 PrimID = PrimaryID; 242 PrimID = PrimaryID;
243 PreviousPV = CurrentPV; 243 PreviousPV = CurrentPV;
244 PreviousUnitID = UnitID; 244 PreviousUnitID = UnitID;
245 } 245 }
246 246
247 G4ThreeVector CCaloSD::SetToLocal(const G4Thre 247 G4ThreeVector CCaloSD::SetToLocal(const G4ThreeVector& global) const{
248 248
249 G4ThreeVector localPoint; 249 G4ThreeVector localPoint;
250 const G4VTouchable* touch= PreStepPoint->G 250 const G4VTouchable* touch= PreStepPoint->GetTouchable();
251 localPoint=touch->GetHistory()->GetTopTransf 251 localPoint=touch->GetHistory()->GetTopTransform().TransformPoint(global);
252 252
253 return localPoint; 253 return localPoint;
254 } 254 }
255 255
256 void CCaloSD::EndOfEvent(G4HCofThisEvent*) { 256 void CCaloSD::EndOfEvent(G4HCofThisEvent*) {
257 summarize(); 257 summarize();
258 } 258 }
259 259
260 void CCaloSD::summarize() { 260 void CCaloSD::summarize() {
261 } 261 }
262 262
263 void CCaloSD::clear() { 263 void CCaloSD::clear() {
264 } 264 }
265 265
266 void CCaloSD::DrawAll() { 266 void CCaloSD::DrawAll() {
267 } 267 }
268 268
269 void CCaloSD::PrintAll() { 269 void CCaloSD::PrintAll() {
270 G4cout << "CCaloSD: Collection " << theHC->G 270 G4cout << "CCaloSD: Collection " << theHC->GetName() << G4endl;
271 theHC->PrintAllHits(); 271 theHC->PrintAllHits();
272 } 272 }
273 273
274 void CCaloSD::SetOrganization(CCalVOrganizatio 274 void CCaloSD::SetOrganization(CCalVOrganization* org){
275 275
276 if (theDescription!=0) 276 if (theDescription!=0)
277 delete theDescription; 277 delete theDescription;
278 theDescription = org; 278 theDescription = org;
279 } 279 }
280 280
281 G4double CCaloSD::curve_LY(const G4StepPoint* 281 G4double CCaloSD::curve_LY(const G4StepPoint* stepPoint) {
282 282
283 G4double weight = 1.; 283 G4double weight = 1.;
284 G4ThreeVector localPoint = SetToLocal(stepPo 284 G4ThreeVector localPoint = SetToLocal(stepPoint->GetPosition());
285 const G4double crlength = 230.; 285 const G4double crlength = 230.;
286 G4double dapd = 0.5 * crlength - localPoint. 286 G4double dapd = 0.5 * crlength - localPoint.z();
287 if (dapd >= -0.1 || dapd <= crlength+0.1) { 287 if (dapd >= -0.1 || dapd <= crlength+0.1) {
288 if (dapd <= 100.) 288 if (dapd <= 100.)
289 weight = 1.05 - dapd * 0.0005; 289 weight = 1.05 - dapd * 0.0005;
290 } else { 290 } else {
291 G4cout << "CCaloSD, light coll curve : wro 291 G4cout << "CCaloSD, light coll curve : wrong distance to APD " << dapd
292 << " crlength = " << crlength 292 << " crlength = " << crlength
293 << " z of localPoint = " << localPo 293 << " z of localPoint = " << localPoint.z()
294 << " take weight = " << weight << G 294 << " take weight = " << weight << G4endl;
295 } 295 }
296 #ifdef ddebug 296 #ifdef ddebug
297 G4cout << "CCaloSD, light coll curve : " << 297 G4cout << "CCaloSD, light coll curve : " << dapd
298 << " crlength = " << crlength 298 << " crlength = " << crlength
299 << " z of localPoint = " << localPoin 299 << " z of localPoint = " << localPoint.z()
300 << " take weight = " << weight << G4e 300 << " take weight = " << weight << G4endl;
301 #endif 301 #endif
302 return weight; 302 return weight;
303 } 303 }
304 304