| Geant4 Cross Reference |
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 //
28 // ------------------------------------------------------------
29 // GEANT 4 class implementation
30 //
31 // ---------------- GFlashShowerModel ----------------
32 //
33 // Authors: E.Barberio & Joanna Weng - 9.11.2004
34 // ------------------------------------------------------------
35
36 #include "G4Electron.hh"
37 #include "G4Positron.hh"
38 #include "G4NeutrinoE.hh"
39 #include "G4NeutrinoMu.hh"
40 #include "G4NeutrinoTau.hh"
41 #include "G4AntiNeutrinoE.hh"
42 #include "G4AntiNeutrinoMu.hh"
43 #include "G4AntiNeutrinoTau.hh"
44 #include "G4PionZero.hh"
45 #include "G4VProcess.hh"
46 #include "G4ios.hh"
47 #include "G4LogicalVolume.hh"
48 #include "geomdefs.hh"
49
50 #include "GFlashShowerModel.hh"
51 #include "GFlashHomoShowerParameterisation.hh"
52 #include "GFlashSamplingShowerParameterisation.hh"
53 #include "GFlashEnergySpot.hh"
54
55 GFlashShowerModel::GFlashShowerModel(G4String modelName, G4Envelope* envelope)
56 : G4VFastSimulationModel(modelName, envelope), PBound(0), Parameterisation(0), HMaker(0)
57 {
58 FlagParamType = 0;
59 FlagParticleContainment = 1;
60 StepInX0 = 0.1;
61 EnergyStop = 0.0;
62 Messenger = new GFlashShowerModelMessenger(this);
63 }
64
65 GFlashShowerModel::GFlashShowerModel(G4String modelName)
66 : G4VFastSimulationModel(modelName), PBound(0), Parameterisation(0), HMaker(0)
67 {
68 FlagParamType = 1;
69 FlagParticleContainment = 1;
70 StepInX0 = 0.1;
71 EnergyStop = 0.0;
72 Messenger = new GFlashShowerModelMessenger(this);
73 }
74
75 GFlashShowerModel::~GFlashShowerModel()
76 {
77 delete Messenger;
78 }
79
80 G4bool GFlashShowerModel::IsApplicable(const G4ParticleDefinition& particleType)
81 {
82 return &particleType == G4Electron::ElectronDefinition()
83 || &particleType == G4Positron::PositronDefinition();
84 }
85
86 /**********************************************************************/
87 /* Checks whether conditions of fast parameterisation are fullfilled */
88 /**********************************************************************/
89
90 G4bool GFlashShowerModel::ModelTrigger(const G4FastTrack& fastTrack)
91
92 {
93 G4bool select = false;
94 if (FlagParamType != 0) {
95 G4double ParticleEnergy = fastTrack.GetPrimaryTrack()->GetKineticEnergy();
96 G4ParticleDefinition& ParticleType = *(fastTrack.GetPrimaryTrack()->GetDefinition());
97 if (ParticleEnergy > PBound->GetMinEneToParametrise(ParticleType)
98 && ParticleEnergy < PBound->GetMaxEneToParametrise(ParticleType))
99 {
100 // check conditions depending on particle flavour
101 // performance to be optimized @@@@@@@
102 Parameterisation->GenerateLongitudinalProfile(ParticleEnergy);
103 select = CheckParticleDefAndContainment(fastTrack);
104 if (select) EnergyStop = PBound->GetEneToKill(ParticleType);
105 }
106 }
107
108 return select;
109 }
110
111 G4bool GFlashShowerModel::CheckParticleDefAndContainment(const G4FastTrack& fastTrack)
112 {
113 G4bool filter = false;
114 G4ParticleDefinition* ParticleType = fastTrack.GetPrimaryTrack()->GetDefinition();
115
116 if (ParticleType == G4Electron::ElectronDefinition()
117 || ParticleType == G4Positron::PositronDefinition())
118 {
119 filter = true;
120 if (FlagParticleContainment == 1) {
121 filter = CheckContainment(fastTrack);
122 }
123 }
124 return filter;
125 }
126
127 G4bool GFlashShowerModel::CheckContainment(const G4FastTrack& fastTrack)
128 {
129 G4bool filter = false;
130 // track informations
131 G4ThreeVector DirectionShower = fastTrack.GetPrimaryTrackLocalDirection();
132 G4ThreeVector InitialPositionShower = fastTrack.GetPrimaryTrackLocalPosition();
133
134 G4ThreeVector OrthoShower, CrossShower;
135 // Returns orthogonal vector
136 OrthoShower = DirectionShower.orthogonal();
137 // Shower in direction perpendicular to OrthoShower and DirectionShower
138 CrossShower = DirectionShower.cross(OrthoShower);
139
140 G4double R = Parameterisation->GetAveR90();
141 G4double Z = Parameterisation->GetAveT90();
142 G4int CosPhi[4] = {1, 0, -1, 0};
143 G4int SinPhi[4] = {0, 1, 0, -1};
144
145 G4ThreeVector Position;
146 G4int NlateralInside = 0;
147 // pointer to solid we're in
148 G4VSolid* SolidCalo = fastTrack.GetEnvelopeSolid();
149 for (int i = 0; i < 4; i++) {
150 // polar coordinates
151 Position = InitialPositionShower + Z * DirectionShower + R * CosPhi[i] * OrthoShower
152 + R * SinPhi[i] * CrossShower;
153
154 if (SolidCalo->Inside(Position) != kOutside) NlateralInside++;
155 }
156
157 // choose to parameterise or flag when all inetc...
158 if (NlateralInside == 4) filter = true;
159 // std::cout << " points = " <<NlateralInside << std::endl;
160 return filter;
161 }
162
163 void GFlashShowerModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep)
164 {
165 // parametrise electrons
166 if (fastTrack.GetPrimaryTrack()->GetDefinition() == G4Electron::ElectronDefinition()
167 || fastTrack.GetPrimaryTrack()->GetDefinition() == G4Positron::PositronDefinition())
168 ElectronDoIt(fastTrack, fastStep);
169 }
170
171 void GFlashShowerModel::ElectronDoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep)
172 {
173 // std::cout<<"--- ElectronDoit --- "<<std::endl;
174
175 fastStep.KillPrimaryTrack();
176 fastStep.ProposePrimaryTrackPathLength(0.0);
177 fastStep.ProposeTotalEnergyDeposited(fastTrack.GetPrimaryTrack()->GetKineticEnergy());
178
179 //-----------------------------
180 // Get track parameters
181 //-----------------------------
182 // E,vect{p} and t,vec(x)
183 G4double Energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy();
184
185 // axis of the shower, in global reference frame:
186 G4ThreeVector DirectionShower = fastTrack.GetPrimaryTrack()->GetMomentumDirection();
187 G4ThreeVector OrthoShower, CrossShower;
188 OrthoShower = DirectionShower.orthogonal();
189 CrossShower = DirectionShower.cross(OrthoShower);
190
191 //--------------------------------
192 /// Generate longitudinal profile
193 //--------------------------------
194 Parameterisation->GenerateLongitudinalProfile(Energy);
195 // performance iteration @@@@@@@
196
197 /// Initialisation of long. loop variables
198 G4VSolid* SolidCalo = fastTrack.GetEnvelopeSolid();
199 G4ThreeVector pos = fastTrack.GetPrimaryTrackLocalPosition();
200 G4ThreeVector dir = fastTrack.GetPrimaryTrackLocalDirection();
201 G4double Bound = SolidCalo->DistanceToOut(pos, dir);
202
203 G4double Dz = 0.00;
204 G4double ZEndStep = 0.00;
205
206 G4double EnergyNow = Energy;
207 G4double EneIntegral = 0.00;
208 G4double LastEneIntegral = 0.00;
209 G4double DEne = 0.00;
210
211 G4double NspIntegral = 0.00;
212 G4double LastNspIntegral = 0.00;
213 G4double DNsp = 0.00;
214
215 // starting point of the shower:
216 G4ThreeVector PositionShower = fastTrack.GetPrimaryTrack()->GetPosition();
217 G4ThreeVector NewPositionShower = PositionShower;
218 G4double StepLenght = 0.00;
219
220 //--------------------------
221 /// Begin Longitudinal Loop
222 //-------------------------
223
224 do {
225 // determine step size=min(1Xo,next boundary)
226 G4double stepLength = StepInX0 * Parameterisation->GetX0();
227 if (Bound < stepLength) {
228 Dz = Bound;
229 Bound = 0.00;
230 }
231 else {
232 Dz = stepLength;
233 Bound = Bound - Dz;
234 }
235 ZEndStep = ZEndStep + Dz;
236
237 // Determine Energy Release in Step
238 if (EnergyNow > EnergyStop) {
239 LastEneIntegral = EneIntegral;
240 EneIntegral = Parameterisation->IntegrateEneLongitudinal(ZEndStep);
241 DEne = std::min(EnergyNow, (EneIntegral - LastEneIntegral) * Energy);
242 LastNspIntegral = NspIntegral;
243 NspIntegral = Parameterisation->IntegrateNspLongitudinal(ZEndStep);
244 DNsp =
245 std::max(1., std::floor((NspIntegral - LastNspIntegral) * Parameterisation->GetNspot()));
246 }
247 // end of the shower
248 else {
249 DEne = EnergyNow;
250 DNsp = std::max(1., std::floor((1. - NspIntegral) * Parameterisation->GetNspot()));
251 }
252 EnergyNow = EnergyNow - DEne;
253
254 // Apply sampling fluctuation - only in sampling calorimeters
255 //
256 GFlashSamplingShowerParameterisation* sp =
257 dynamic_cast<GFlashSamplingShowerParameterisation*>(Parameterisation);
258 if (sp) {
259 G4double DEneSampling = sp->ApplySampling(DEne, Energy);
260 DEne = DEneSampling;
261 }
262
263 // move particle in the middle of the step
264 StepLenght = StepLenght + Dz / 2.00;
265 NewPositionShower = NewPositionShower + StepLenght * DirectionShower;
266 StepLenght = Dz / 2.00;
267
268 // generate spots & hits:
269 for (G4int i = 0; i < DNsp; ++i) {
270 GFlashEnergySpot Spot;
271
272 // Spot energy: the same for all spots
273 Spot.SetEnergy(DEne / DNsp);
274 G4double PhiSpot = Parameterisation->GeneratePhi(); // phi of spot
275 G4double RSpot = Parameterisation // radius of spot
276 ->GenerateRadius(i, Energy, ZEndStep - Dz / 2.);
277
278 // check reference-> may be need to introduce rot matrix @@@
279 // Position: equally spaced in z
280
281 G4ThreeVector SpotPosition =
282 NewPositionShower + Dz / DNsp * DirectionShower * (i + 1 / 2. - DNsp / 2.)
283 + RSpot * std::cos(PhiSpot) * OrthoShower + RSpot * std::sin(PhiSpot) * CrossShower;
284 Spot.SetPosition(SpotPosition);
285
286 // Generate Hits of this spot
287 HMaker->make(&Spot, &fastTrack);
288 }
289 } while (EnergyNow > 0.0 && Bound > 0.0);
290
291 //---------------
292 /// End Loop
293 //---------------
294 }
295
296 /*
297
298 void
299 GFlashShowerModel::GammaDoIt(const G4FastTrack& fastTrack,
300 G4FastStep& fastStep)
301 {
302
303 if( fastTrack.GetPrimaryTrack()->GetKineticEnergy() > EnergyStop )
304 return;
305
306 //deposita in uno spot unico l'energia
307 //con andamento exp decrescente.
308
309 // Kill the particle to be parametrised
310 fastStep.KillPrimaryTrack();
311 fastStep.SetPrimaryTrackPathLength(0.0);
312 fastStep.SetTotalEnergyDeposited(fastTrack.GetPrimaryTrack()
313 ->GetKineticEnergy());
314 // other settings????
315 feSpotList.clear();
316
317 //-----------------------------
318 // Get track parameters
319 //-----------------------------
320
321 // E,vect{p} and t,vec(x)
322 G4double Energy =
323 fastTrack.GetPrimaryTrack()->GetKineticEnergy();
324 // axis of the shower, in global reference frame:
325 G4ThreeVector DirectionShower =
326 fastTrack.GetPrimaryTrack()->GetMomentumDirection();
327 // starting point of the shower:
328 G4ThreeVector PositionShower =
329 fastTrack.GetPrimaryTrack()->GetPosition();
330
331 //G4double DEneSampling = Parameterisation->ApplySampling(Energy,Energy);
332 //if(DEneSampling <= 0.00) DEneSampling=Energy;
333
334 if(Energy > 0.0)
335 {
336 G4double dist = Parameterisation->GenerateExponential(Energy);
337
338 GFlashEnergySpot Spot;
339 Spot.SetEnergy( Energy );
340 G4ThreeVector SpotPosition = PositionShower + dist*DirectionShower;
341 Spot.SetPosition(SpotPosition);
342
343 // Record the Spot:
344 feSpotList.push_back(Spot);
345
346 //Generate Hits of this spot
347 HMaker->make(Spot);
348 }
349 }
350
351 */
352