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