Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ExExChProcessChanneling.cc
Go to the documentation of this file.
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 
31 
32 #include "G4Step.hh"
33 #include "G4StepPoint.hh"
34 #include "G4VParticleChange.hh"
35 #include "Randomize.hh"
36 #include "G4RandomDirection.hh"
37 #include "G4RandomTools.hh"
38 
39 
40 #include "G4Box.hh"
41 #include "G4Tubs.hh"
42 
43 #include "G4Material.hh"
44 #include "G4NistManager.hh"
46 #include "G4Navigator.hh"
47 #include "G4GeometryTolerance.hh"
48 
49 #include "G4SystemOfUnits.hh"
50 #include "G4PhysicalConstants.hh"
51 
52 #include "XLatticeManager3.hh"
53 
54 #include "XLogicalAtomicLattice.hh"
56 #include "XLogicalBase.hh"
57 #include "XUnitCell.hh"
58 
60 
63 G4VDiscreteProcess(aName){
64  fLatticeManager = XLatticeManager3::GetXLatticeManager();
65 
68  G4cout<<"\n ExExChProcessChanneling::Constructor:";
69  G4cout<<"Geometry surface tolerance is: ";
70  G4cout<< kCarTolerance / mm << " mm"<<G4endl;
71  if(verboseLevel>1) {
72  G4cout << GetProcessName() << " is created "<< G4endl;
73  }
74 
75  fFileCharacteristicsName = "";
76  fTimeStepMax = 0.;
77  fTimeStepMin = 2.E2 * CLHEP::angstrom;
78  fTransverseVariationMax = 2.E-2 * CLHEP::angstrom;
79 
80  fPointYPre = -1.;
81  fPointYPost = -1.;
82 }
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90 
93 G4VDiscreteProcess(right){
94  ;
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 
100  return fPotentialEnergy;
101 }
102 
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104 
107  fPotentialEnergy = vPotential;
108 }
109 
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 
113  return fElectricField;
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117 
120  fElectricField = vElectricField;
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124 
127  return fIntegratedDensity;
128 }
129 
130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
131 
134  fIntegratedDensity = vIntegratedDensity;
135 }
136 
137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138 
140  return fNucleiDensity;
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
144 
147  fNucleiDensity = vNucleiDensity;
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151 
153  return fElectronDensity;
154 }
155 
156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157 
160  fElectronDensity = vElectronDensity;
161 }
162 
163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
164 
165 G4bool ExExChProcessChanneling::UpdateInitialParameters(const G4Track& aTrack){
166 
167  if(GetInfo(aTrack)->GetPositionChanneled().x() == DBL_MAX){
168  // when the particle enters the crystal the position in the channel
169  //is randomly generated using a uniform distribution
170  G4double vXposition = G4UniformRand() *
171  GetXPL(aTrack)->ComputeInterplanarPeriod();
172 
173  //vXposition = 1.0 * CLHEP::angstrom;
174 
175  //initial position in the channel is stored
176  GetInfo(aTrack)->SetPositionChanneled(G4ThreeVector(vXposition,
177  0.,
178  0.));
179  GetInfo(aTrack)->SetPositionChanneledInitial(G4ThreeVector(vXposition,
180  0.,
181  0.));
182  }
183 
184  if(GetInfo(aTrack)->GetMomentumChanneledInitial().x() == DBL_MAX){
185  // the first time it enter the crystal we take the momentum
186  // for the post step which is the only one in the crystal
187  G4ThreeVector vMomentum =
188  ComputeMomentum(aTrack,aTrack.GetStep()->GetPostStepPoint());
189 
190  GetInfo(aTrack)->SetMomentumChanneled(vMomentum);
191  GetInfo(aTrack)->SetMomentumChanneledInitial(vMomentum);
192  return true;
193  }
194 
195  return false;
196 }
197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
198 
199 void ExExChProcessChanneling::UpdateParameters(const G4Track& aTrack){
200 
201  if(fIntegratedDensity->HasBeenInitialized(GetXPL(aTrack))
202  == false){
203  ComputeCrystalCharacteristic(aTrack);
204  G4cout << "ChannelingProcess::UpdatePositionMomentumDensity::";
205  G4cout<<"fIntegratedDensity->Initialized" << G4endl;
206  }
207 
208  if(UpdateInitialParameters(aTrack) == false){
209  G4ThreeVector vMomentumNew =
210  ComputeMomentum(aTrack,
211  aTrack.GetStep()->GetPreStepPoint());
212  GetInfo(aTrack)->SetMomentumChanneled(vMomentumNew);
213  }
214  G4ThreeVector vPositionPost =
215  ComputePositionInTheCrystal(aTrack.GetStep()->GetPostStepPoint(),aTrack);
216  G4ThreeVector vPositionPre =
217  ComputePositionInTheCrystal(aTrack.GetStep()->GetPreStepPoint(),aTrack);
218 
219  fHasToComputeTrajectory = true;
220 
221  if(vPositionPost.y() == fPointYPost &&
222  vPositionPre.y() == fPointYPre){
223  fHasToComputeTrajectory = false;
224  }
225  else{
226  if(GetXPL(aTrack)->IsBent()){
227  fPointYPre = vPositionPre.y();
228  fPointYPost = vPositionPost.y();
229  }
230  else{
231  fPointYPre = vPositionPre.z();
232  fPointYPost = vPositionPost.z();
233  }
234  }
235 
236  G4ThreeVector fMom = GetInfo(aTrack)->GetMomentumChanneled();
237  G4ThreeVector fPos = GetInfo(aTrack)->GetPositionChanneled();
238  G4ThreeVector fMomHalf = GetInfo(aTrack)->GetMomentumChanneled();
239  G4ThreeVector fPosHalf = GetInfo(aTrack)->GetPositionChanneled();
240 
241  if(GetXPL(aTrack)->IsBent()){
242  fIntegrationPeriod = (vPositionPost.phi() - vPositionPre.phi())*
243  GetXPL(aTrack)->GetCurvatureRadius().x();
244  fIntegrationPeriod = vPositionPost.y() - vPositionPre.y();
245  }
246  else{
247  fIntegrationPeriod = vPositionPost.z() - vPositionPre.z();
248  }
249 
250  fTimeStepTotal = 0.;
251 
252  G4double vNucleiDensity=0.;
253  G4double vElectronDensity=0.;
254 
255  if(fIntegrationPeriod>0. && fHasToComputeTrajectory==true){
256  G4double kBeta = 0.;
257  G4double kPos = 0.;
258  G4double kMom = 0.;
259  G4double kBR = 0.;
260  G4double Z = 0.;
261  do{
262  UpdateIntegrationStep(aTrack,fMom);
263 
264  fPosHalf = fPos;
265  fMomHalf = fMom;
266 
267  kBeta = aTrack.GetVelocity()/c_light;
268  if(fMom.z()!=0.){
269  kPos = fTimeStep / fMom.z();
270  }
271  else{
272  kPos = fTimeStep / 1.E-20;
273  }
274  kMom = fTimeStep / kBeta;
275  kBR = fTimeStep * (fMom.z() * kBeta);;
276  Z = GetParticleDefinition(aTrack)->GetPDGCharge();
277 
278  fPosHalf += (fMom * kPos * 0.5);
279  fMomHalf +=
280  (GetElectricField()->GetEC(fPos,GetXPL(aTrack))
281  * Z * kMom * 0.5);
282 
283  if(GetXPL(aTrack)->IsBent()){
284  G4double temp =
285  fMomHalf.x() + kBR * 0.5 /
286  (GetXPL(aTrack)->GetCurvatureRadius()).x();
287  fMomHalf.setX(temp);
288  }
289 
290  fPos += (fMomHalf * kPos);
291  fMom +=
292  (GetElectricField()->GetEC(fPosHalf,GetXPL(aTrack))
293  * Z * kMom );
294 
295  if(GetXPL(aTrack)->IsBent()){
296  G4double temp =
297  fMom.x() + kBR /
298  (GetXPL(aTrack)->GetCurvatureRadius()).x();
299  fMom.setX(temp);
300  }
301 
302  fTimeStepTotal += fTimeStep;
303 
304  vNucleiDensity +=
305  (fTimeStep *
306  (GetNucleiDensity()->GetEC(fPos,GetXPL(aTrack)).x()
307  +GetNucleiDensity()->GetEC(fPos,GetXPL(aTrack)).x()
308  ) * 0.5);
309  vElectronDensity +=
310  (fTimeStep * (
311  GetElectronDensity()->GetEC(fPos,GetXPL(aTrack)).x() +
312  GetElectronDensity()->GetEC(fPos,GetXPL(aTrack)).x()) * 0.5);
313 
314 
315  } while(fTimeStepTotal<fIntegrationPeriod);
316 
317  vNucleiDensity /= fIntegrationPeriod;
318  vElectronDensity /= fIntegrationPeriod;
319  GetInfo(aTrack)->SetNucleiDensity(vNucleiDensity);
320  GetInfo(aTrack)->SetElectronDensity(vElectronDensity);
321  }
322  else{
323  ResetDensity(aTrack);
324  }
325 
326  GetInfo(aTrack)->SetMomentumChanneled(fMom);
327  GetInfo(aTrack)->SetPositionChanneled(fPos);
328 
329 }
330 
331 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
332 
333 void ExExChProcessChanneling::ResetDensity(const G4Track& aTrack){
334  GetInfo(aTrack)->SetNucleiDensity(1.);
335  GetInfo(aTrack)->SetElectronDensity(1.);
336 }
337 
338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
339 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
340 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
341 
342 
343 G4ThreeVector ExExChProcessChanneling::
344 ComputePositionInTheCrystal(G4StepPoint* vStep,const G4Track& aTrack){
345 
346  G4StepPoint* vStepVol = CheckStepPointLatticeForVolume(vStep,aTrack);
347  G4StepPoint* vStepPos = CheckStepPointLatticeForPosition(vStep,aTrack);
348 
349  G4TouchableHistory* theTouchable =
350  (G4TouchableHistory*)(vStepVol->GetTouchable());
351  G4ThreeVector vWorldPos = vStepPos->GetPosition();
352  G4ThreeVector vLocalPos =
353  theTouchable->GetHistory()->GetTopTransform().TransformPoint(vWorldPos);
354 
355 
356  if(GetXPL(aTrack)->IsBent() == false){
357  G4Box* vXtalSolid =
358  (G4Box*) vStepVol->GetPhysicalVolume()
359  ->GetLogicalVolume()->GetSolid();
360  vLocalPos += G4ThreeVector(vXtalSolid->GetXHalfLength(),
361  vXtalSolid->GetYHalfLength(),
362  vXtalSolid->GetZHalfLength());
363  }
364  else{
365  G4Tubs* vXtalSolid =
366  (G4Tubs*) vStepVol->GetPhysicalVolume()
367  ->GetLogicalVolume()->GetSolid();
368  vLocalPos.rotateZ(-vXtalSolid->GetStartPhiAngle());
369  }
370  return vLocalPos;
371 }
372 
373 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
374 
375 G4StepPoint* ExExChProcessChanneling::
376 CheckStepPointLatticeForVolume(G4StepPoint* vStep, const G4Track& aTrack){
377 
378  G4StepPoint* vStepPre = aTrack.GetStep()->GetPreStepPoint();
379  G4StepPoint* vStepPost = aTrack.GetStep()->GetPostStepPoint();
380 
381  if( fLatticeManager->HasLattice(vStep->GetPhysicalVolume()) ) {
382  return vStep;
383  }
384  else if(fLatticeManager->HasLattice(vStepPost->GetPhysicalVolume())
385  == false
386  && fLatticeManager->HasLattice(vStepPre->GetPhysicalVolume())
387  == true
388  && vStep == vStepPost &&
389  vStepPost->GetStepStatus() == fGeomBoundary) {
390  return vStepPre;
391  }
392  else if(fLatticeManager->HasLattice(vStepPre->GetPhysicalVolume())
393  == false &&
394  fLatticeManager->HasLattice(vStepPost->GetPhysicalVolume())
395  == true &&
396  vStep == vStepPre &&
397  vStepPre->GetStepStatus() == fGeomBoundary) {
398  return vStepPost;
399  }
400  else{
401  return vStep;
402  }
403 }
404 
405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
406 
407 G4StepPoint* ExExChProcessChanneling::
408 CheckStepPointLatticeForPosition(G4StepPoint* vStep, const G4Track& aTrack){
409 
410  G4StepPoint* vStepPre = aTrack.GetStep()->GetPreStepPoint();
411  G4StepPoint* vStepPost = aTrack.GetStep()->GetPostStepPoint();
412 
413  if( fLatticeManager->HasLattice(vStep->GetPhysicalVolume()) ) {
414  return vStep;
415  }
416  else if(fLatticeManager->HasLattice(vStepPost->GetPhysicalVolume())
417  == false &&
418  fLatticeManager->HasLattice(vStepPre->GetPhysicalVolume())
419  == true &&
420  vStep == vStepPost &&
421  vStepPost->GetStepStatus() == fGeomBoundary) {
422  return vStepPost;
423  }
424  else if(fLatticeManager->HasLattice(vStepPre->GetPhysicalVolume())
425  == false &&
426  fLatticeManager->HasLattice(vStepPost->GetPhysicalVolume())
427  == true &&
428  vStep == vStepPre &&
429  vStepPre->GetStepStatus() == fGeomBoundary) {
430  return vStepPre;
431  }
432  else{
433  return vStep;
434  }
435 }
436 
437 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
438 
439 G4bool ExExChProcessChanneling::
440 UpdateIntegrationStep(const G4Track& aTrack,G4ThreeVector& vMom){
441 
442  if(vMom.x() != 0.0 || vMom.y() != 0.0){
443  double xy2 = vMom.x() * vMom.x() + vMom.y()*vMom.y();
444 
445  if(xy2!=0.){
446  fTimeStep =
447  std::fabs(fTransverseVariationMax *
448  aTrack.GetStep()->GetPreStepPoint()->GetTotalEnergy() /
449  std::pow(xy2,0.5));
450  if(fTimeStep < fTimeStepMin) fTimeStep = fTimeStepMin;
451  else{
452  fTimeStepMax = std::sqrt( fTransverseVariationMax *
453  aTrack.GetStep()->GetPreStepPoint()->GetTotalEnergy() /
454  std::fabs(fElectricField->GetMaximum(GetXPL(aTrack))));
455 
456  if(fTimeStep > fTimeStepMax) fTimeStep = fTimeStepMax;
457  }
458  }
459  else{
460  fTimeStep = fTimeStepMin;
461  }
462 
463  if(fTimeStep + fTimeStepTotal > fIntegrationPeriod){
464  fTimeStep = fIntegrationPeriod - fTimeStepTotal;
465  }
466 
467  return true;
468  }
469  else{
470  fTimeStep = fTimeStepMin;
471  }
472  return false;
473 }
474 
475 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
476 
477 G4double ExExChProcessChanneling::
478 GetChannelingMeanFreePath(const G4Track& aTrack){
479  //----------------------------------------
480  // return the channeling MFP
481  //----------------------------------------
482 
483 
484  G4double vMFP = 0.1 * ComputeOscillationPeriod(aTrack);
485 
486  if(HasLatticeOnBoundaryPre(aTrack) == true){
487  vMFP = 0.001 * ComputeOscillationPeriod(aTrack);
488  }
489 
490  return vMFP;
491 }
492 
493 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
494 
496 GetMeanFreePath(const G4Track& aTrack,
497  G4double, // previousStepSize
499 
500  //----------------------------------------
501  // the condition is forced to check if
502  // the volume has a lattice at each step.
503  // if it hasn't, return DBL_MAX
504  //----------------------------------------
505 
506  *condition = Forced;
507 
508  if(HasLattice(aTrack)){
509  GetInfo(aTrack)->SetInTheCrystal(true);
510  return GetChannelingMeanFreePath(aTrack);
511  }
512  else{
513  GetInfo(aTrack)->SetInTheCrystal(false);
514  return DBL_MAX;
515  }
516 }
517 
518 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
519 
521 PostStepDoIt(const G4Track& aTrack,
522  const G4Step&){
523 
524  //----------------------------------------
525  // check if the volume has a lattice
526  // and if the particle is in channeling.
527  // If it is so, the particle is forced
528  // to follow the channeling plane
529  // direction. If the particle has
530  // dechanneled or exited the crystal,
531  // the outgoing angle is evaluated
532  //----------------------------------------
533 
534  aParticleChange.Initialize(aTrack);
535 
536  GetInfo(aTrack)->StoreDensityPreviousStep();
537 
538  if((HasLattice(aTrack) == true) &&
539  (HasLatticeOnBoundaryPost(aTrack) == false)){
540  UpdateParameters(aTrack);
541 
542  G4ThreeVector vMomentum =
543  GetInfo(aTrack)->GetMomentumChanneled().unit();
544 
545  G4ThreeVector vPosition;
546  vPosition =
547  ComputePositionInTheCrystal(aTrack.GetStep()->GetPostStepPoint(),
548  aTrack);
549 
550  GetXPL(aTrack)->ProjectMomentumVectorFromLatticeToWorld(vMomentum,
551  vPosition);
552 
554  }
555  else{
556  // if the volume has no lattice it resets the density factors
557  ResetDensity(aTrack);
558  GetInfo(aTrack)->SetMomentumChanneled(G4ThreeVector(DBL_MAX,
559  DBL_MAX,
560  DBL_MAX));
561 
562  GetInfo(aTrack)->SetPositionChanneled(G4ThreeVector(DBL_MAX,
563  DBL_MAX,
564  DBL_MAX));
565  }
566 
567  return &aParticleChange;
568 }
569 
570 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
571 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
572 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
573 
574 G4ThreeVector ExExChProcessChanneling::
575 ComputeTransverseEnergy(const G4Track& aTrack){
576  //----------------------------------------
577  // compute the particle transverse energy
578  // in the crystal reference system
579  //----------------------------------------
580 
581  G4ThreeVector vTransverseEnergy = ComputePotentialEnergy(aTrack)
582  + ComputeKineticEnergy(aTrack);
583  return vTransverseEnergy;
584 }
585 
586 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
587 G4ThreeVector ExExChProcessChanneling::
588 ComputeKineticEnergy(const G4Track& aTrack){
589  //----------------------------------------
590  // compute the particle kinetic energy
591  // in the crystal reference system
592  //----------------------------------------
593 
594  G4double vTotalEnergy =
595  aTrack.GetStep()->GetPostStepPoint()->GetTotalEnergy();
596 
597  G4ThreeVector vMom = GetInfo(aTrack)->GetMomentumChanneled();
598 
599  G4ThreeVector vKineticEnergy = 0.5 *
600  G4ThreeVector((vMom.x() * vMom.x()) / vTotalEnergy,
601  0.,
602  (vMom.z() * vMom.z()) / vTotalEnergy);
603 
604  return vKineticEnergy;
605 }
606 
607 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
608 
609 G4ThreeVector ExExChProcessChanneling::
610 ComputePotentialEnergy(const G4Track& aTrack){
611  //----------------------------------------
612  // compute the particle transverse energy
613  // in the crystal reference system
614  //----------------------------------------
615 
616 
617  G4ThreeVector vPotentialEnergy =
618  fPotentialEnergy->GetEC(GetInfo(aTrack)->GetPositionChanneled(),
619  GetXPL(aTrack));
620 
621  vPotentialEnergy *= GetParticleDefinition(aTrack)->GetPDGCharge();
622 
623  return vPotentialEnergy;
624 }
625 
626 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
627 
628 G4ThreeVector ExExChProcessChanneling::
629 ComputeCentrifugalEnergy(const G4Track& aTrack,G4ThreeVector vPosition){
630  //----------------------------------------
631  // compute the transverse energy variation
632  // in the crystal reference system
633  //----------------------------------------
634 
635  G4double vTotalEnergy =
636  aTrack.GetStep()->GetPreStepPoint()->GetTotalEnergy();
637 
638  G4double vPositionX = vPosition.x();
639 
640  if(ParticleIsNegative(aTrack) && false){
641  vPositionX -=
642  GetXPL(aTrack)->ComputeInterplanarPeriod() * 0.5;
643  }
644  G4ThreeVector vEnergyVariation = G4ThreeVector();;
645  if(GetXPL(aTrack)->IsBent()){
646  vEnergyVariation =
647  G4ThreeVector(vTotalEnergy * vPositionX /
648  GetXPL(aTrack)->GetCurvatureRadius().x(),
649  0.,
650  0.);
651  }
652 
653  return vEnergyVariation;
654 }
655 
656 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
657 
658 G4ThreeVector ExExChProcessChanneling::
659 ComputeMomentum(const G4Track& aTrack,G4StepPoint* vStep){
660  //----------------------------------------
661  // compute the particle momentum
662  // in the crystal reference system
663  //----------------------------------------
664 
665  G4ThreeVector vPosition = ComputePositionInTheCrystal(vStep,aTrack);
666 
667  G4ThreeVector vMomentum = aTrack.GetMomentum();
668 
669  GetXPL(aTrack)->
670  ProjectMomentumVectorFromWorldToLattice(vMomentum,vPosition);
671 
672  return vMomentum;
673 }
674 
675 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
676 
677 G4ThreeVector ExExChProcessChanneling::
678 ComputeCentrifugalEnergyMaximumVariation(const G4Track& aTrack){
679  //----------------------------------------
680  // compute the transverse energy variation
681  // in the crystal reference system
682  //----------------------------------------
683 
684  G4double vTotalEnergy =
685  aTrack.GetStep()->GetPreStepPoint()->GetTotalEnergy();
686 
687  G4ThreeVector vEnergyVariation = G4ThreeVector();
688  if(GetXPL(aTrack)->IsBent()){
689  vEnergyVariation = G4ThreeVector(vTotalEnergy *
690  GetXPL(aTrack)->ComputeInterplanarPeriod() /
691  GetXPL(aTrack)->GetCurvatureRadius().x(),
692  0.,
693  0.);
694  }
695 
696  return vEnergyVariation;
697 }
698 
699 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
700 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
701 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
702 
703 
704 G4double ExExChProcessChanneling::
705 ComputeCriticalEnergyMaximum(const G4Track& aTrack){
706  //----------------------------------------
707  // compute the critical energy
708  // for channeling
709  //----------------------------------------
710 
711  G4double vCriticalEnergy = 0.;
712 
713  if(ParticleIsNegative(aTrack)){
714  vCriticalEnergy =
715  - fPotentialEnergy->GetMinimum(GetXPL(aTrack));
716  }
717  else{
718  vCriticalEnergy =
719  + fPotentialEnergy->GetMaximum(GetXPL(aTrack));
720  }
721 
722  vCriticalEnergy *= std::fabs(GetParticleDefinition(aTrack)->GetPDGCharge());
723 
724  return vCriticalEnergy;
725 }
726 
727 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
728 
729 G4double ExExChProcessChanneling::
730 ComputeCriticalEnergyMinimum(const G4Track& aTrack){
731  //----------------------------------------
732  // compute the critical energy minimum
733  // for channeling
734  //----------------------------------------
735 
736  G4double vCriticalEnergy = 0.;
737 
738  if(ParticleIsNegative(aTrack)){
739  vCriticalEnergy =
740  - fPotentialEnergy->GetMaximum(GetXPL(aTrack));
741  }
742  else{
743  vCriticalEnergy =
744  + fPotentialEnergy->GetMinimum(GetXPL(aTrack));
745  }
746 
747  vCriticalEnergy *= std::fabs(GetParticleDefinition(aTrack)->GetPDGCharge());
748 
749  return vCriticalEnergy;
750 }
751 
752 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
753 
754 G4double ExExChProcessChanneling::
755 ComputeCriticalAngle(const G4Track& aTrack){
756  //----------------------------------------
757  // compute the critical angle
758  // for chenneling
759  //----------------------------------------
760 
761  G4double vTotalEnergy =
762  aTrack.GetStep()->GetPostStepPoint()->GetTotalEnergy();
763  G4double vCriticalAngle = std::pow( 2.0 *
764  std::fabs( ( ComputeCriticalEnergyMaximum(aTrack)
765  - ComputeCriticalEnergyMinimum(aTrack) )
766  / vTotalEnergy ) , 0.5);
767  return vCriticalAngle;
768 }
769 
770 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
771 
772 G4double ExExChProcessChanneling::
773 ComputeOscillationPeriod(const G4Track& aTrack){
774  //----------------------------------------
775  // compute the particle oscillation
776  // period in the crystal channel
777  //----------------------------------------
778 
779  G4double vInterplanarPeriod =
780  GetXPL(aTrack)->ComputeInterplanarPeriod();
781  G4double vOscillationPeriod =
782  CLHEP::pi * vInterplanarPeriod / ComputeCriticalAngle(aTrack);
783  return vOscillationPeriod;
784 }
785 
786 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
787 
788 G4double ExExChProcessChanneling::
789 ComputeCriticalRadius(const G4Track& aTrack){
790  //----------------------------------------
791  // compute the critical radius
792  // for channeling
793  //----------------------------------------
794 
795  G4double vTotalEnergy =
796  aTrack.GetStep()->GetPreStepPoint()->GetTotalEnergy();
797 
798  G4double vCriticalRadius = 1.E-20;
799  if(fElectricField->GetMaximum(GetXPL(aTrack)) != 0.){
800  vCriticalRadius = vTotalEnergy
801  / fElectricField->GetMaximum(GetXPL(aTrack));
802  }
803  return vCriticalRadius;
804 }
805 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
806 
807 G4ThreeVector ExExChProcessChanneling::
808 ComputePotentialWellCentre(const G4Track& aTrack){
809  //----------------------------------------
810  // compute the central point of
811  // the potential well for channeling
812  //----------------------------------------
813 
814  G4double vInterplanarPeriodHalf =
815  0.5 * GetXPL(aTrack)->ComputeInterplanarPeriod();
816 
817  G4double vCentreX = vInterplanarPeriodHalf;
818 
819  if(GetXPL(aTrack)->IsBent()){
820  G4double vTotalEnergy =
821  aTrack.GetStep()->GetPreStepPoint()->GetTotalEnergy();
822 
823  G4double vPotentialWellDepth =
824  ComputeCriticalEnergyMaximum(aTrack)
825  - (ComputeCriticalEnergyMinimum(aTrack));
826 
827  vCentreX *= (1. - 0.5 * vTotalEnergy /
828  vPotentialWellDepth /
829  GetXPL(aTrack)->GetCurvatureRadius().x() *
830  vInterplanarPeriodHalf );
831  }
832 
833  G4ThreeVector vCentre = G4ThreeVector(vCentreX,0.,0.);
834 
835 
836  return vCentre;
837 }
838 
839 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
840 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
841 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
842 
843 
846  return(aPD.GetPDGCharge() != 0.);
847 }
848 
849 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
850 
853 
854 }
855 
856 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
857 
858 XPhysicalLattice* ExExChProcessChanneling::
859 GetXPL(const G4Track& aTrack){
860  if(fLatticeManager->HasLattice(
861  aTrack.GetStep()->GetPostStepPoint()->GetPhysicalVolume())){
862  return fLatticeManager->GetXPhysicalLattice(aTrack.GetStep()
864  }
865  else if(fLatticeManager->HasLattice(aTrack.GetStep()->GetPreStepPoint()
866  ->GetPhysicalVolume()) &&
867  aTrack.GetStep()->GetPostStepPoint()->GetStepStatus()
868  == fGeomBoundary) {
869  return fLatticeManager->GetXPhysicalLattice(aTrack.GetStep()
871  }
872  else{
873  G4cout << "LATTICE NOT FOUND: ERROR" << G4endl;
874  return NULL;
875  }
876 }
877 
878 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
879 
880 G4bool ExExChProcessChanneling::HasLattice(const G4Track& aTrack){
881  if(fLatticeManager->HasLattice(aTrack.GetStep()
883  return true;
884  }
885  else{
886  return false;
887  }
888 }
889 
890 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
891 
892 G4bool ExExChProcessChanneling::
893 HasLatticeOnBoundary(const G4Track& aTrack){
894  if(HasLatticeOnBoundaryPost(aTrack) || HasLatticeOnBoundaryPre(aTrack)) {
895  return true;
896  }
897  else{
898  return false;
899  }
900 }
901 
902 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
903 
904 G4bool ExExChProcessChanneling::
905 HasLatticeOnBoundaryPre(const G4Track& aTrack){
906  if(fLatticeManager->HasLattice(aTrack.GetStep()->GetPreStepPoint()->
907  GetPhysicalVolume()) &&
909  return true;
910  }
911  else{
912  return false;
913  }
914 }
915 
916 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
917 
918 G4bool ExExChProcessChanneling::
919 HasLatticeOnBoundaryPost(const G4Track& aTrack){
920  if(fLatticeManager->HasLattice(aTrack.GetStep()->
921  GetPostStepPoint()->GetPhysicalVolume()) &&
923  return true;
924  }
925  else{
926  return false;
927  }
928 }
929 
930 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
931 
932 G4bool ExExChProcessChanneling::ParticleIsNegative(const G4Track& aTrack){
933  if(GetParticleDefinition(aTrack)->GetPDGCharge() < 0.) {
934  return true;
935  }
936  else{
937  return false;
938  }
939 }
940 
941 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
942 
943 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
944 
945 G4bool ExExChProcessChanneling::
946 ParticleIsNotOnBoundaryPre(const G4Track& aTrack){
947  if(aTrack.GetStep()->GetPreStepPoint()->GetStepStatus() != fGeomBoundary){
948  return true;
949  }
950  else{
951  return false;
952  }
953 }
954 
955 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
956 
957 G4bool ExExChProcessChanneling::
958 ParticleIsNotOnBoundaryPost(const G4Track& aTrack){
959  if(aTrack.GetStep()->GetPostStepPoint()->GetStepStatus() != fGeomBoundary){
960  return true;
961  }
962  else{
963  return false;
964  }
965 }
966 
967 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
968 
969 G4bool ExExChProcessChanneling::
970 ParticleIsNotOnBoundary(const G4Track& aTrack){
971  if(ParticleIsNotOnBoundaryPost(aTrack) ||
972  ParticleIsNotOnBoundaryPre(aTrack)){
973  return true;
974  }
975  else{
976  return false;
977  }
978 }
979 
980 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
981 
982 ExExChParticleUserInfo* ExExChProcessChanneling::
983 GetInfo(const G4Track& aTrack){
984  return (ExExChParticleUserInfo*) aTrack.GetUserInformation();
985 }
986 
987 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
988 
989 G4ParticleDefinition* ExExChProcessChanneling::
990 GetParticleDefinition(const G4Track& aTrack){
991  return const_cast<G4ParticleDefinition*>(aTrack.GetParticleDefinition());
992 }
993 
994 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
995 
996 void ExExChProcessChanneling::
997 ComputeCrystalCharacteristic(const G4Track& aTrack){
998  fIntegratedDensity->SetXPhysicalLattice(GetXPL(aTrack));
999  fIntegratedDensity->InitializeTables();
1000 
1001  if(fFileCharacteristicsName != ""){
1002  G4String filename;
1003  fElectronDensity->InitializePhysicalLattice(GetXPL(aTrack));
1004  fNucleiDensity->InitializePhysicalLattice(GetXPL(aTrack));
1005  fPotentialEnergy->InitializePhysicalLattice(GetXPL(aTrack));
1006  fElectricField->InitializePhysicalLattice(GetXPL(aTrack));
1007 
1008  filename = fFileCharacteristicsName + "_pot.txt";
1009  G4cout << filename << G4endl;
1010  fPotentialEnergy->ReadFromECHARM(filename,CLHEP::eV);
1011  fPotentialEnergy->PrintOnFile("temppot.dat",GetXPL(aTrack));
1012 
1013  filename = fFileCharacteristicsName + "_efx.txt";
1014  G4cout << filename << G4endl;
1015  fElectricField->ReadFromECHARM(filename,CLHEP::eV/CLHEP::m);
1016  fElectricField->PrintOnFile("tempefx.dat",GetXPL(aTrack));
1017 
1018  filename = fFileCharacteristicsName + "_atd.txt";
1019  G4cout << filename << G4endl;
1020  fNucleiDensity->ReadFromECHARM(filename);
1021 
1022  filename = fFileCharacteristicsName + "_eld.txt";
1023  G4cout << filename << G4endl;
1024  fElectronDensity->ReadFromECHARM(filename);
1025  fIntegratedDensity->ReadFromFiles(fFileCharacteristicsName);
1026  }
1027  else{
1028  fPotentialEnergy->InitializePhysicalLattice(
1029  GetXPL(aTrack));
1030  fElectricField->InitializePhysicalLattice(GetXPL(aTrack));
1031  }
1032 }
1033 
1034 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
G4double GetTotalEnergy() const
G4double condition(const G4ErrorSymMatrix &m)
static constexpr double m
G4double GetXHalfLength() const
static XLatticeManager3 * GetXLatticeManager()
static constexpr double mm
Definition: G4SIunits.hh:115
G4int verboseLevel
Definition: G4VProcess.hh:368
void SetMomentumChanneledInitial(G4ThreeVector)
CLHEP::Hep3Vector G4ThreeVector
G4double GetVelocity() const
XVCrystalCharacteristic * GetNucleiDensity()
double x() const
Definition of the XLogicalBase class.
Definition: G4Box.hh:64
Definition of the XLatticeManager3 class.
G4StepStatus GetStepStatus() const
Definition: G4Tubs.hh:85
G4double GetSurfaceTolerance() const
G4VSolid * GetSolid() const
ExExChProcessChanneling(const G4String &processName="ExExChProcessChanneling")
G4ThreeVector GetCurvatureRadius()
tuple x
Definition: test.py:50
Definition of the XLogicalAtomicLatticeDiamond class.
const G4VTouchable * GetTouchable() const
const G4Step * GetStep() const
G4double GetZHalfLength() const
virtual G4double GetMinimum(XPhysicalLattice *)
double z() const
void setX(double)
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
G4VUserTrackInformation * GetUserInformation() const
G4bool HasBeenInitialized(XPhysicalLattice *)
G4StepPoint * GetPreStepPoint() const
Definition of the ExExChParticleUserInfo class.
void SetXPhysicalLattice(XPhysicalLattice *)
void SetMomentumChanneled(G4ThreeVector)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
XVCrystalCharacteristic * GetElectronDensity()
G4VPhysicalVolume * GetPhysicalVolume() const
const G4ThreeVector & GetPosition() const
Definition of the ExExChProcessChanneling class.
void SetPositionChanneled(G4ThreeVector)
bool G4bool
Definition: G4Types.hh:79
XVCrystalCharacteristic * GetElectricField()
void SetPotential(XVCrystalCharacteristic *)
Hep3Vector & rotateZ(double)
Definition: ThreeVector.cc:110
const G4ParticleDefinition * GetParticleDefinition() const
G4double GetYHalfLength() const
virtual void PrintOnFile(const G4String &, XPhysicalLattice *, G4double=1)=0
const G4double kCarTolerance
Definition: G4Step.hh:76
G4double GetStartPhiAngle() const
double phi() const
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
static constexpr double eV
G4ThreeVector GetEC(G4ThreeVector, XPhysicalLattice *)
void SetIntegratedDensity(XCrystalIntegratedDensityHub *)
G4ThreeVector ProjectMomentumVectorFromLatticeToWorld(G4ThreeVector &, G4ThreeVector &)
G4ThreeVector GetMomentum() const
virtual void Initialize(const G4Track &)
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
void SetElectricField(XVCrystalCharacteristic *)
void InitializePhysicalLattice(XPhysicalLattice *)
G4LogicalVolume * GetLogicalVolume() const
XPhysicalLattice * GetXPhysicalLattice(G4VPhysicalVolume *)
virtual G4double GetMaximum(XPhysicalLattice *)
XVCrystalCharacteristic * GetPotential()
Hep3Vector unit() const
virtual G4double GetMeanFreePath(const G4Track &, G4double, G4ForceCondition *)
G4StepPoint * GetPostStepPoint() const
XCrystalIntegratedDensityHub * GetIntegratedDensity()
void SetNucleiDensity(XVCrystalCharacteristic *)
double y() const
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
Definition of the XLogicalAtomicLattice class.
void SetElectronDensity(XVCrystalCharacteristic *)
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
#define G4endl
Definition: G4ios.hh:61
const G4AffineTransform & GetTopTransform() const
const G4NavigationHistory * GetHistory() const
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
G4ForceCondition
G4double GetPDGCharge() const
#define DBL_MAX
Definition: templates.hh:83
Definition of the XUnitCell class.
G4double ComputeInterplanarPeriod()
void SetPositionChanneledInitial(G4ThreeVector)
virtual G4bool IsApplicable(const G4ParticleDefinition &)
float c_light
Definition: hepunit.py:257
static constexpr double pi
Definition: SystemOfUnits.h:54
static G4GeometryTolerance * GetInstance()
static constexpr double angstrom
Definition: SystemOfUnits.h:82
virtual void ReadFromECHARM(const G4String &, G4double=1)=0
bool HasLattice(G4VPhysicalVolume *)