ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Channeling.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Channeling.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 //
26 
27 #include "G4Channeling.hh"
28 
29 #include "Randomize.hh"
30 
31 #include "G4ChannelingTrackData.hh"
32 #include "G4TouchableHistory.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4LambdacPlus.hh"
35 
36 G4Channeling::G4Channeling():
37 G4VDiscreteProcess("channeling"),
38 fChannelingID(-1),
39 fTimeStepMin(0.),
40 fTimeStepMax(0.),
41 fTransverseVariationMax(2.E-2 * CLHEP::angstrom),
42 k010(G4ThreeVector(0.,1.,0.)){
43  fChannelingID = G4PhysicsModelCatalog::GetIndex("channeling");
44  if(fChannelingID == -1){
45  fChannelingID = G4PhysicsModelCatalog::Register("channeling");
46  }
47  fSpin = G4ThreeVector(0.,0.,0.);
48 }
49 
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 
52 G4Channeling::~G4Channeling(){;}
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55 
56 G4ChannelingTrackData* G4Channeling::GetTrackData(const G4Track& aTrack){
57  G4ChannelingTrackData* trackdata =
58  (G4ChannelingTrackData*)(aTrack.GetAuxiliaryTrackInformation(fChannelingID));
59  if(trackdata == nullptr){
60  trackdata = new G4ChannelingTrackData();
61  aTrack.SetAuxiliaryTrackInformation(fChannelingID,trackdata);
62  }
63  return trackdata;
64 }
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67 
68 void G4Channeling::GetEF(const G4Track& aTrack,
69  G4ThreeVector& pos,
70  G4ThreeVector& out){
71  out = G4ThreeVector((GetMatData(aTrack)->GetEFX()->GetEC(pos)),
72  (GetMatData(aTrack)->GetEFY()->GetEC(pos)),
73  0.);
74 }
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 
78 void G4Channeling::PosToLattice(G4StepPoint* step,G4ThreeVector& pos){
79  G4TouchableHistory* theTouchable = (G4TouchableHistory*)(step->GetTouchable());
80 
81  pos -= theTouchable->GetTranslation();
82  pos = ((*theTouchable->GetRotation()).inverse())(pos);
83 }
84 
85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86 
87 G4bool G4Channeling::UpdateParameters(const G4Track& aTrack){
88 
89  G4LogicalCrystalVolume* aLCV = (G4LogicalCrystalVolume*)(aTrack.GetVolume()->GetLogicalVolume());
90 
91  G4StepPoint* postStepPoint = aTrack.GetStep()->GetPostStepPoint();
92  G4StepPoint* preStepPoint = aTrack.GetStep()->GetPreStepPoint();
93 
94  G4ThreeVector posPost = postStepPoint->GetPosition();
95  aLCV->RotateToLattice(posPost);
96  G4ThreeVector posPre = preStepPoint->GetPosition();
97  aLCV->RotateToLattice(posPre);
98 
99  G4double integrationLimit = std::fabs(posPost.z() - posPre.z());
100 
101  if(integrationLimit > 0.){
102  //----------------------------------------
103  // Check if the crystal is bent
104  //----------------------------------------
105  G4bool isBent = GetMatData(aTrack)->IsBent();
106 
107  //----------------------------------------
108  // Get the momentum in the world reference
109  // frame and rotate to the solid reference frame
110  //----------------------------------------
111  G4TouchableHistory* theTouchable = (G4TouchableHistory*)(preStepPoint->GetTouchable());
112  G4ThreeVector momWorld = aTrack.GetStep()->GetPreStepPoint()->GetMomentum();
113  G4ThreeVector mom = (*theTouchable->GetRotation())(momWorld);
114 
115  //----------------------------------------
116  // Get the momentum in the solid reference
117  // frame and rotate to the crystal reference frame
118  //----------------------------------------
119  aLCV->RotateToLattice(mom);
120 
121  //----------------------------------------
122  // Get the momentum in the crystal reference
123  // frame and rotate to the reference frame
124  // solidal to the bent planes
125  //----------------------------------------
126  if(isBent){
127  PosToLattice(preStepPoint,posPre);
128  G4ThreeVector axis010 = (*theTouchable->GetRotation())(k010);
129  mom.rotate(axis010,-posPre.z()/GetMatData(aTrack)->GetBR(posPre).x());
130  }
131 
132  //----------------------------------------
133  // Take the position stored in the track data.
134  // If the particle enters the crystal,
135  // the position in the channel is randomly
136  // generated using a uniform distribution
137  //----------------------------------------
138  G4ThreeVector pos;
139  if(GetTrackData(aTrack)->GetPosCh().x() == DBL_MAX){
140  G4double posX = G4UniformRand() * GetMatData(aTrack)->GetPot()->GetIntSp(0);
141  G4double posY = G4UniformRand() * GetMatData(aTrack)->GetPot()->GetIntSp(1);
142  pos = G4ThreeVector(posX,posY,0.);
143  }
144  else{
145  pos = GetTrackData(aTrack)->GetPosCh();
146  }
147 
148  G4double step=0., stepTot=0.;
149  G4double nud =0., eld =0.;
150  G4double efx =0., efy =0.;
151  G4double nud_temp =0., eld_temp =0.;
152 
153  G4double beta = aTrack.GetVelocity()/CLHEP::c_light;
154  G4double Z = GetParticleDefinition(aTrack)->GetPDGCharge();
155 
156  const G4double oneSixth = 1./6.;
157  G4ThreeVector posk1,posk2,posk3,posk4,posk5,posk6;
158  G4ThreeVector momk1,momk2,momk3,momk4,momk5,momk6;
159  G4ThreeVector pos_temp, efxy;
160 
161  do{
162  //----------------------------------------
163  // Limit the variable step length for the
164  // integration via the selected algorithm
165  // and update variables for the integration
166  //----------------------------------------
167 
168  UpdateIntegrationStep(aTrack,mom,step);
169  if(step + stepTot > integrationLimit){
170  step = integrationLimit - stepTot;
171  }
172 
173  //----------------------------------------
174  // Function integration algorithm
175  // 4th Order Runge-Kutta
176  //----------------------------------------
177 
178  GetEF(aTrack,pos,efxy);
179  posk1 = step / mom.z() * mom;
180  momk1 = step / beta * Z * efxy;
181  if(isBent) momk1.setX(momk1.x() - step * mom.z() * beta / (GetMatData(aTrack)->GetBR(pos)).x());
182 
183  GetEF(aTrack,pos_temp = pos + posk1 * 0.5,efxy);
184  posk2 = step / mom.z() * (mom + momk1 * 0.5);
185  momk2 = step / beta * Z * efxy;
186  if(isBent) momk2.setX(momk2.x() - step * mom.z() * beta / (GetMatData(aTrack)->GetBR(pos_temp)).x());
187 
188  GetEF(aTrack,pos_temp = pos + posk2 * 0.5,efxy);
189  posk3 = step / mom.z() * (mom + momk2 * 0.5);
190  momk3 = step / beta * Z * efxy;
191  if(isBent) momk3.setX(momk3.x() - step * mom.z() * beta / (GetMatData(aTrack)->GetBR(pos_temp)).x());
192 
193  GetEF(aTrack,pos_temp = pos + posk3,efxy);
194  posk4 = step / mom.z() * (mom + momk3);
195  momk4 = step / beta * Z * efxy;
196  if(isBent) momk4.setX(momk4.x() - step * mom.z() * beta / (GetMatData(aTrack)->GetBR(pos_temp)).x());
197 
198  pos = pos + oneSixth * (posk1 + 2.*posk2 + 2.*posk3 + posk4);
199  mom = mom + oneSixth * (momk1 + 2.*momk2 + 2.*momk3 + momk4);
200 
201  //----------------------------------------
202  // Function integration algorithm
203  // 2th Order Velocity-Verlet
204  //----------------------------------------
205 
206  /*
207  GetEF(aTrack,pos,efxy);
208  posk1 = pos + (step * 0.5 / mom.z()) * mom;
209  //momk1 = mom + step * 0.5 / betaZ * efxy;
210  momk1 = mom;
211  if(isBent) momk1.setX(momk1.x() - step * 0.5 * mom.z() * beta / (GetMatData(aTrack)->GetBR(pos)).x());
212 
213  GetEF(aTrack,posk1,efxy);
214  pos = pos + (step / momk1.z()) * momk1;
215  //mom = mom + step / betaZ * efxy;
216  mom = mom;
217  if(isBent) mom.setX(mom.x() - step * mom.z() * beta / (GetMatData(aTrack)->GetBR(posk1)).x());
218  */
219 
220  //----------------------------------------
221  // Update the total step and the electron
222  // and nuclei density experienced by
223  // the particle during its motion
224  //----------------------------------------
225 
226  stepTot += step;
227 
228  nud_temp = GetMatData(aTrack)->GetNuD()->GetEC(pos);
229  eld_temp = GetMatData(aTrack)->GetElD()->GetEC(pos);
230 
231  if(nud_temp < 0.) {nud_temp = 0.;}
232  if(eld_temp < 0.) {eld_temp = 0.;}
233 
234  nud += (step * nud_temp);
235  eld += (step * eld_temp);
236 
237  efx += (step * GetMatData(aTrack)->GetEFX()->GetEC(pos));
238  efy += (step * GetMatData(aTrack)->GetEFY()->GetEC(pos));
239 
240  } while(stepTot<integrationLimit);
241 
242  nud /= stepTot;
243  eld /= stepTot;
244 
245  if(nud < 1.E-10) {nud = 1.E-10;}
246  if(eld < 1.E-10) {eld = 1.E-10;}
247 
248  GetTrackData(aTrack)->SetNuD(nud);
249  GetTrackData(aTrack)->SetElD(eld);
250 
251  GetTrackData(aTrack)->SetEFX(efx);
252  GetTrackData(aTrack)->SetEFY(efy);
253 
254  GetTrackData(aTrack)->SetMomCh(mom);
255  GetTrackData(aTrack)->SetPosCh(pos);
256  return true;
257  }
258  else{
259  return false;
260  }
261 
262  return false;
263 }
264 
265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
266 
267 G4bool G4Channeling::
268 UpdateIntegrationStep(const G4Track& aTrack,
269  G4ThreeVector& mom,
270  G4double& step){
271 
272  if(mom.x() != 0.0 || mom.y() != 0.0){
273  double xy2 = mom.x() * mom.x() + mom.y()*mom.y();
274 
275  if(xy2!=0.){
276  step = std::fabs(fTransverseVariationMax * GetPre(aTrack)->GetKineticEnergy() / std::pow(xy2,0.5));
277  if(step < fTimeStepMin) step = fTimeStepMin;
278  else{
279  fTimeStepMax = std::sqrt( fTransverseVariationMax * GetPre(aTrack)->GetKineticEnergy()
280  / std::fabs(GetMatData(aTrack)->GetEFX()->GetMax()));
281 
282  if(step > fTimeStepMax) step = fTimeStepMax;
283  }
284  }
285  else{
286  step = fTimeStepMin;
287  }
288 
289  return true;
290  }
291  else{
292  step = fTimeStepMin;
293  }
294  return false;
295 }
296 
297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
298 
299 G4double G4Channeling::
300 GetMeanFreePath(const G4Track& aTrack,
301  G4double, // previousStepSize
302  G4ForceCondition* condition){
303 
304  //----------------------------------------
305  // the condition is forced to check if
306  // the volume has a lattice at each step.
307  // if it hasn't, return DBL_MAX
308  //----------------------------------------
309 
310  *condition = Forced;
311 
312  G4LogicalVolume* aLV = aTrack.GetVolume()->GetLogicalVolume();
313  G4LogicalVolume* aNLV = aTrack.GetNextVolume()->GetLogicalVolume();
314 
315  if(G4LogicalCrystalVolume::IsLattice(aLV) == true &&
316  G4LogicalCrystalVolume::IsLattice(aNLV) == true){
317  G4double osc_per = GetOscillationPeriod(aTrack);
318  fTimeStepMin = osc_per * 2.E-4;
319  return osc_per * 0.01;
320  }
321  else{
322  GetTrackData(aTrack)->Reset();
323  return DBL_MAX;
324  }
325 }
326 
327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
328 
329 G4VParticleChange* G4Channeling::
330 PostStepDoIt(const G4Track& aTrack,
331  const G4Step&){
332 
333  //----------------------------------------
334  // check if the volume has a lattice
335  // and if the particle is in channeling.
336  // If it is so, the particle is forced
337  // to follow the channeling plane
338  // direction. If the particle has
339  // dechanneled or exited the crystal,
340  // the outgoing angle is evaluated
341  //----------------------------------------
342 
343  aParticleChange.Initialize(aTrack);
344  G4LogicalVolume* aLV = aTrack.GetVolume()->GetLogicalVolume();
345  G4LogicalVolume* aNLV = aTrack.GetNextVolume()->GetLogicalVolume();
346 
347 
348  if(G4LogicalCrystalVolume::IsLattice(aLV) == true &&
349  G4LogicalCrystalVolume::IsLattice(aNLV) == true){
350 
351  G4bool bModifiedTraj = UpdateParameters(aTrack);
352 
353  if(bModifiedTraj==true){
354  //----------------------------------------
355  // Get the momentum in the reference frame
356  // solidal to the bent planes and rotate
357  // to the reference frame
358  //----------------------------------------
359  G4LogicalCrystalVolume* aLCV = (G4LogicalCrystalVolume*)(aTrack.GetVolume()->GetLogicalVolume());
360  G4ThreeVector momCh = GetTrackData(aTrack)->GetMomCh();
361 
362  G4StepPoint* postStepPoint = aTrack.GetStep()->GetPostStepPoint();
363  G4TouchableHistory* theTouchable = (G4TouchableHistory*)(postStepPoint->GetTouchable());
364 
365  if(GetMatData(aTrack)->IsBent()){
366  G4ThreeVector posPost = postStepPoint->GetPosition();
367  PosToLattice(postStepPoint,posPost);
368  G4ThreeVector axis010 = (*theTouchable->GetRotation())(k010);
369  momCh.rotate(axis010,posPost.z()/GetMatData(aTrack)->GetBR(posPost).x());
370  }
371 
372  //----------------------------------------
373  // Get the momentum in the crystal reference
374  // frame and rotate to the solid reference frame
375  //----------------------------------------
376 
377  aLCV->RotateToSolid(momCh);
378 
379  //----------------------------------------
380  // Get the momentum in the solid reference
381  // frame and rotate to the world reference frame
382  //----------------------------------------
383  G4ThreeVector mom = ((*theTouchable->GetRotation()).inverse())(momCh);
384 
385  aParticleChange.ProposeMomentumDirection(mom.unit());
386  aParticleChange.ProposePolarization(fSpin);
387  }
388  }
389  else{
390  // if the volume has no lattice it resets the density factors
391  GetTrackData(aTrack)->Reset();
392  }
393 
394  return &aParticleChange;
395 }
396 
397 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....