source: trunk/source/parameterisations/gflash/src/GFlashShowerModel.cc@ 848

Last change on this file since 848 was 823, checked in by garnier, 17 years ago

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