ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Par01EMShowerModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Par01EMShowerModel.cc
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 //
28 //
29 //
30 //
31 #include "Par01EMShowerModel.hh"
32 #include "Par01EnergySpot.hh"
33 
34 #include "Randomize.hh"
35 
36 #include "G4Electron.hh"
37 #include "G4Positron.hh"
38 #include "G4Gamma.hh"
39 #include "G4TransportationManager.hh"
40 #include "G4VSensitiveDetector.hh"
41 #include "G4TouchableHandle.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4NistManager.hh"
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 
48 Par01EMShowerModel::Par01EMShowerModel(G4String modelName, G4Region* envelope)
49 : G4VFastSimulationModel(modelName, envelope)
50 {
51  fFakeStep = new G4Step();
52  fFakePreStepPoint = fFakeStep->GetPreStepPoint();
53  fFakePostStepPoint = fFakeStep->GetPostStepPoint();
54  fTouchableHandle = new G4TouchableHistory();
55  fpNavigator = new G4Navigator();
56  fNaviSetup = false;
57  fCsI = 0;
58 }
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
62 Par01EMShowerModel::Par01EMShowerModel(G4String modelName)
63 : G4VFastSimulationModel(modelName)
64 {
65  fFakeStep = new G4Step();
66  fFakePreStepPoint = fFakeStep->GetPreStepPoint();
67  fFakePostStepPoint = fFakeStep->GetPostStepPoint();
68  fTouchableHandle = new G4TouchableHistory();
69  fpNavigator = new G4Navigator();
70  fNaviSetup = false;
71  fCsI = 0;
72 }
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75 
76 Par01EMShowerModel::~Par01EMShowerModel()
77 {
78  delete fFakeStep;
79  delete fpNavigator;
80 }
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83 
84 G4bool Par01EMShowerModel::IsApplicable(const G4ParticleDefinition& particleType)
85 {
86  return
87  &particleType == G4Electron::ElectronDefinition() ||
88  &particleType == G4Positron::PositronDefinition() ||
89  &particleType == G4Gamma::GammaDefinition();
90 }
91 
92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
93 
94 G4bool Par01EMShowerModel::ModelTrigger(const G4FastTrack& fastTrack)
95 {
96  // Applies the parameterisation above 100 MeV:
97  return fastTrack.GetPrimaryTrack()->GetKineticEnergy() > 100*MeV;
98 }
99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101 
102 void Par01EMShowerModel::DoIt(const G4FastTrack& fastTrack,
103  G4FastStep& fastStep)
104 {
105  // Kill the parameterised particle:
106  fastStep.KillPrimaryTrack();
107  fastStep.ProposePrimaryTrackPathLength(0.0);
108  fastStep.ProposeTotalEnergyDeposited(fastTrack.GetPrimaryTrack()->GetKineticEnergy());
109 
110  // split into "energy spots" energy according to the shower shape:
111  Explode(fastTrack);
112 
113  // and put those energy spots into the crystals:
114  BuildDetectorResponse();
115 
116 }
117 
118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
119 
120 void Par01EMShowerModel::Explode(const G4FastTrack& fastTrack)
121 {
122  //-----------------------------------------------------
123  //
124  //-----------------------------------------------------
125 
126  // Reduced quantities:
127  // -- critical energy in CsI:
128  G4double Ec = 800*MeV/(54. + 1.2); // 54 = mean Z of CsI
129  G4double Energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy();
130  G4double y = Energy/Ec;
131 
132  // compute value of parameter "a" of longitudinal profile, b assumed = 0.5
133  G4double a, tmax, b(0.5), C;
134  if (fastTrack.GetPrimaryTrack()->GetDefinition() == G4Gamma::GammaDefinition()) C = 0.5;
135  else C = -0.5;
136  tmax = 1.0 * (std::log(y) + C);
137  a = 1.0 + b*tmax;
138 
139  // t : reduced quantity = z/X0:
140  G4double t, bt;
141  if ( fCsI == 0 ) fCsI = G4NistManager::Instance()->FindOrBuildMaterial("G4_CESIUM_IODIDE");
142  G4double X0 = fCsI->GetRadlen();
143  // Moliere radius:
144  G4double Es = 21*MeV;
145  G4double Rm = X0*Es/Ec;
146 
147  // axis of the shower, in global reference frame:
148  G4ThreeVector xShower, yShower, zShower;
149  zShower = fastTrack.GetPrimaryTrack()->GetMomentumDirection();
150  xShower = zShower.orthogonal();
151  yShower = zShower.cross(xShower);
152  // starting point of the shower:
153  G4ThreeVector sShower = fastTrack.GetPrimaryTrack()->GetPosition();
154 
155  // We shoot 100 spots of energy:
156  G4int nSpots = 100;
157  G4double deposit = Energy/double(nSpots);
158  Par01EnergySpot eSpot;
159  eSpot.SetEnergy(deposit);
160  G4ThreeVector ePoint;
161  G4double z, r, phi;
162 
163  feSpotList.clear();
164  for (int i = 0; i < nSpots; i++)
165  {
166  // Longitudinal profile:
167  // -- shoot z according to Gamma distribution:
168  bt = G4RandGamma::shoot(a,1.0);
169  t = bt/b;
170  z = t*X0;
171 
172  // transverse profile:
173  // we set 90% of energy in one Rm,
174  // the rest between 1 and 3.5 Rm:
175  G4double xr = G4UniformRand();
176  if (xr < 0.9) r = xr/0.9*Rm;
177  else r = ((xr - 0.9)/0.1*2.5 + 1.0)*Rm;
178  phi = G4UniformRand()*twopi;
179 
180  // build the position:
181  ePoint = sShower +
182  z*zShower +
183  r*std::cos(phi)*xShower + r*std::sin(phi)*yShower;
184 
185  // and the energy spot:
186  eSpot.SetPosition(ePoint);
187 
188  // Records the eSpot:
189  feSpotList.push_back(eSpot);
190  }
191 }
192 
193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
194 
195 void Par01EMShowerModel::BuildDetectorResponse()
196 {
197  // Does the assignation of the energy spots to the sensitive volumes:
198  for (size_t i = 0; i < feSpotList.size(); i++)
199  {
200  // Draw the energy spot:
201  // feSpotList[i].Draw();
202  // feSpotList[i].Print();
203 
204  // "converts" the energy spot into the fake
205  // G4Step to pass to sensitive detector:
206  AssignSpotAndCallHit(feSpotList[i]);
207  }
208 }
209 
210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
211 
212 void Par01EMShowerModel::AssignSpotAndCallHit(const Par01EnergySpot &eSpot)
213 {
214  //
215  // "converts" the energy spot into the fake
216  // G4Step to pass to sensitive detector:
217  //
218  FillFakeStep(eSpot);
219 
220  //
221  // call sensitive part: taken/adapted from the stepping:
222  // Send G4Step information to Hit/Dig if the volume is sensitive
223  //
224  G4VPhysicalVolume* pCurrentVolume =
225  fFakeStep->GetPreStepPoint()->GetPhysicalVolume();
226  G4VSensitiveDetector* pSensitive;
227 
228  if( pCurrentVolume != 0 )
229  {
230  pSensitive = pCurrentVolume->GetLogicalVolume()->
231  GetSensitiveDetector();
232  if( pSensitive != 0 )
233  {
234  pSensitive->Hit(fFakeStep);
235  }
236  }
237 }
238 
239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
240 
241 void Par01EMShowerModel::FillFakeStep(const Par01EnergySpot &eSpot)
242 {
243  //-----------------------------------------------------------
244  // find in which volume the spot is.
245  //-----------------------------------------------------------
246  if (!fNaviSetup)
247  {
248  fpNavigator->
249  SetWorldVolume(G4TransportationManager::GetTransportationManager()->
250  GetNavigatorForTracking()->GetWorldVolume());
251  fpNavigator->
252  LocateGlobalPointAndUpdateTouchableHandle(eSpot.GetPosition(),
253  G4ThreeVector(0.,0.,0.),
254  fTouchableHandle,
255  false);
256  fNaviSetup = true;
257  }
258  else
259  {
260  fpNavigator->
261  LocateGlobalPointAndUpdateTouchableHandle(eSpot.GetPosition(),
262  G4ThreeVector(0.,0.,0.),
263  fTouchableHandle);
264  }
265  //--------------------------------------
266  // Fills attribute of the G4Step needed
267  // by our sensitive detector:
268  //-------------------------------------
269  // set touchable volume at PreStepPoint:
270  fFakePreStepPoint->SetTouchableHandle(fTouchableHandle);
271  // set total energy deposit:
272  fFakeStep->SetTotalEnergyDeposit(eSpot.GetEnergy());
273 }