source: trunk/examples/advanced/composite_calorimeter/src/CCaloSD.cc@ 1313

Last change on this file since 1313 was 807, checked in by garnier, 17 years ago

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