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