Geant4  10.01.p03
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 }
74 
75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 
78 }
79 
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 
84 G4VDiscreteProcess(right){
85  ;
86 }
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 
91  return fPotentialEnergy;
92 }
93 
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
95 
98  fPotentialEnergy = vPotential;
99 }
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102 
104  return fElectricField;
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
108 
111  fElectricField = vElectricField;
112 }
113 
114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
115 
118  return fIntegratedDensity;
119 }
120 
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
122 
125  fIntegratedDensity = vIntegratedDensity;
126 }
127 
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
129 
131  return fNucleiDensity;
132 }
133 
134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135 
138  fNucleiDensity = vNucleiDensity;
139 }
140 
141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142 
144  return fElectronDensity;
145 }
146 
147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
148 
151  fElectronDensity = vElectronDensity;
152 }
153 
154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155 
157 
159  == false){
161  G4cout << "ChannelingProcess::UpdatePositionMomentumDensity::";
162  G4cout<<"fIntegratedDensity->Initialized" << G4endl;
163  }
164 
165  UpdatePosition(aTrack);
166  UpdateMomentum(aTrack);
167  UpdateDensity(aTrack);
168 }
169 
170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
171 
173 
174  if(GetInfo(aTrack)->GetPositionChanneledInitial().x() == DBL_MAX ||
175  HasLatticeOnBoundaryPost(aTrack)){
176  // when the particle enters the crystal the position in the channel
177  //is randomly generated using a uniform distribution
178  G4double vXposition = G4UniformRand() *
180 
181  //initial position in the channel is stored
183  0.,
184  0.));
185 
186  //initial position for the measurement of apparent centrifugal force
187  //is zero at crystal entrance
188  GetInfo(aTrack)->SetPositionChanneled(G4ThreeVector(0.,0.,0.));
189  }
190  else{
191  G4double vPositionX = GetInfo(aTrack)->GetPositionChanneled().x();
192 
193  //if the particle has been under channeling the position
194  //for the measurement of the apparent centrifugal force is reset
195  if(GetInfo(aTrack)->HasBeenUnderCoherentEffect() == 1){
196  GetInfo(aTrack)->SetPositionChanneled(G4ThreeVector(0.,0.,0.));
197  }
198  else{
199  //if the particle has not been under channeling the distance
200  //between the new and the old point is computed for the evaluation
201  //of the centrifugal potential acting on the particle
202  vPositionX += (ComputePositionInTheCrystal(
203  aTrack.GetStep()->GetPostStepPoint(),
204  aTrack).x() - ComputePositionInTheCrystal(
205  aTrack.GetStep()->GetPreStepPoint(),
206  aTrack).x());
207  GetInfo(aTrack)->SetPositionChanneled(G4ThreeVector(vPositionX,
208  0.,
209  0.));
210  }
211  }
212 
213 }
214 
215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
216 
218  if(GetInfo(aTrack)->GetMomentumChanneledInitial().x() == DBL_MAX){
219  // the first time it enter the crystal we take the momentum
220  // for the post step which is the only one in the crystal
221  G4ThreeVector vMomentum =
222  ComputeMomentum(aTrack,aTrack.GetStep()->GetPostStepPoint());
223 
224  GetInfo(aTrack)->SetMomentumChanneled(vMomentum);
225 
227  GetInfo(aTrack)->GetMomentumChanneled());
228  }
229  else{
230  // we take the PREVIOUS step point to compare,
231  // otherwise the momentum is not computed correctly
232  G4ThreeVector vMomentum =
233  G4ThreeVector(GetInfo(aTrack)->GetMomentumChanneled().x(),
234  GetInfo(aTrack)->GetMomentumChanneled().y(),0.);
235 
236  vMomentum+=ComputeMomentum(aTrack,aTrack.GetStep()->GetPreStepPoint());
237 
238  GetInfo(aTrack)->SetMomentumChanneled(vMomentum);
239  }
240 }
241 
242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
243 
245 
246 
247  if(GetInfo(aTrack)->HasBeenUnderCoherentEffect() == 1){
248 
249  G4double vTransverseEnergy = ComputeTransverseEnergy(aTrack).x();
250 
251  if(GetXPhysicalLattice(aTrack)->IsBent()){
252  if(GetInfo(aTrack)->HasBeenUnderCoherentEffect() == 1){
253  if(ParticleIsNegative(aTrack)){
254  vTransverseEnergy -=
256  }
257  else{
258  vTransverseEnergy +=
260  }
261  }
262  }
263 
264  G4double vCharge = GetParticleDefinition(aTrack)->GetPDGCharge();
265  G4double vNucleiDensity =
267  GetXPhysicalLattice(aTrack),
268  G4int(vCharge));
269  G4double vElectronDensity =
271  GetXPhysicalLattice(aTrack),
272  G4int(vCharge));
273 
274  G4double vLowerBoundNegative = 1.;
275  G4double vLowerBoundPositive = 0.01;
276 
277  if(ParticleIsNegative(aTrack)){
278  if(vNucleiDensity < vLowerBoundNegative)
279  {vNucleiDensity = vLowerBoundNegative;}
280  if(vElectronDensity < vLowerBoundNegative)
281  {vElectronDensity = vLowerBoundNegative;}
282  }
283  else{
284  if(vNucleiDensity < vLowerBoundPositive)
285  {vNucleiDensity = vLowerBoundPositive;}
286  if(vElectronDensity < vLowerBoundPositive)
287  {vElectronDensity = vLowerBoundPositive;}
288  }
289 
290  GetInfo(aTrack)->SetNucleiDensity(vNucleiDensity);
291  GetInfo(aTrack)->SetElectronDensity(vElectronDensity);
292  }
293  else{
294  ResetDensity(aTrack);
295  }
296 }
297 
298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
299 
301  GetInfo(aTrack)->SetNucleiDensity(1.);
302  GetInfo(aTrack)->SetElectronDensity(1.);
303 }
304 
305 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
306 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
308 
311 
312  G4StepPoint* vStepPre = aTrack.GetStep()->GetPreStepPoint();
313  G4StepPoint* vStepPost = aTrack.GetStep()->GetPostStepPoint();
314 
315  G4double vTotalEnergy = vStepPre->GetTotalEnergy();
316 
317  G4double vTransverseEnergyX = std::fabs(ComputeTransverseEnergy(aTrack).x());
318  G4double vTransverseEnergyY = std::fabs(ComputeTransverseEnergy(aTrack).y());
319  double vPotentialEnergyX = 0.;
320  double vPotentialEnergyY = 0.;
321 
322  bool bExit = false;
323  do{
324  G4double vXposition = G4UniformRand() *
326 
328  0.,
329  0.));
330 
331  vPotentialEnergyX = ComputePotentialEnergy(aTrack).x();
332  vPotentialEnergyY = ComputePotentialEnergy(aTrack).y();
333  if(vPotentialEnergyX<=vTransverseEnergyX &&
334  vPotentialEnergyY<=vTransverseEnergyY){
335  bExit = true;
336  }
337  } while(bExit == false);
338 
339  vTransverseEnergyX-=vPotentialEnergyX;
340  vTransverseEnergyY-=vPotentialEnergyY;
341 
342  G4double vChAngleX = std::pow(+ 2. * std::fabs(vTransverseEnergyX)
343  / vTotalEnergy , 0.5);
344  G4double vChAngleY = std::pow(+ 2. * std::fabs(vTransverseEnergyY)
345  / vTotalEnergy , 0.5);
346 
347  G4double vPhi = 2. * ( G4UniformRand() - 0.5) * vChAngleX;
348  G4double vTheta = 2. * ( G4UniformRand() - 0.5) * vChAngleY;
349 
350  G4ThreeVector vNewMomentum =
351  G4ThreeVector(0.,0.,1.).rotate(G4ThreeVector(0,1,0),- vPhi)
352  .rotate(G4ThreeVector(1,0,0),- vTheta);
353  G4ThreeVector vPosition = ComputePositionInTheCrystal(vStepPost,aTrack);
354 
355  return GetXPhysicalLattice(aTrack)->
356  ProjectMomentumVectorFromLatticeToWorld(vNewMomentum,vPosition);
357 }
358 
359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
360 
363 
364  G4StepPoint* vStep = aTrack.GetStep()->GetPostStepPoint();
365 
366  G4double vVrAngle = 0.;
367 
368  if(GetXPhysicalLattice(aTrack)->IsBent()) {
369  G4double vRadiusX = GetXPhysicalLattice(aTrack)->
370  GetCurvatureRadius().x();
371 
372  G4double vTotalEnergy = vStep->GetTotalEnergy();
373 
374  G4double vEnergyMax =
375  std::fabs(ComputeCriticalEnergyMaximum(aTrack)
376  - ComputeCriticalEnergyMinimum(aTrack));
377 
378  G4double vEnergyRMS =
379  std::fabs(ComputeCentrifugalEnergyMaximumVariation(aTrack).x());
380 
381  G4double vTransverseEnergy =
382  vEnergyMax + (G4UniformRand() * std::fabs(vEnergyRMS) );
383 
384  vVrAngle = - std::fabs(vRadiusX)/vRadiusX *
385  std::pow(+ 2. * std::fabs(vTransverseEnergy) / vTotalEnergy , 0.5);
386 
387  if(ParticleIsNegative(aTrack)){
388  vVrAngle *= 0.8; // = see PLB 681 (2009) 233
389  }
390  else{
391  vVrAngle *= 1.4;
392  }
393 
394  G4ThreeVector vMomentumChanneled =
395  GetInfo(aTrack)->GetMomentumChanneled();
396  G4double vAngleRatio =
397  (vMomentumChanneled.x()/vTotalEnergy)/ComputeCriticalAngle(aTrack);
398 
399  if(std::fabs(vAngleRatio)<1.5){
400  vVrAngle *= (-(std::fabs(vAngleRatio) - 1.5)/3.);
401  }
402  }
403 
404  G4double vOmega = GetXPhysicalLattice(aTrack)->GetLatticeAngles().y();
405  G4double vPhi = vVrAngle * std::cos(vOmega);
406  G4double vTheta = vVrAngle * std::sin(vOmega);
407 
408  G4ThreeVector vNewMomentum =
409  aTrack.GetMomentum().unit()
410  .rotate(G4ThreeVector(0.,1.,0.), - vPhi)
411  .rotate(G4ThreeVector(1.,0.,0.), -vTheta);
412 
413  return vNewMomentum.unit();
414 }
415 
416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
418 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
419 
420 
423 
424  G4StepPoint* vStepVol = CheckStepPointLatticeForVolume(vStep,aTrack);
425  G4StepPoint* vStepPos = CheckStepPointLatticeForPosition(vStep,aTrack);
426 
427  G4TouchableHistory* theTouchable =
428  (G4TouchableHistory*)(vStepVol->GetTouchable());
429  G4ThreeVector vWorldPos = vStepPos->GetPosition();
430  G4ThreeVector vLocalPos =
431  theTouchable->GetHistory()->GetTopTransform().TransformPoint(vWorldPos);
432 
433 
434  if(GetXPhysicalLattice(aTrack)->IsBent() == false){
435  G4Box* vXtalSolid =
436  (G4Box*) vStepVol->GetPhysicalVolume()
437  ->GetLogicalVolume()->GetSolid();
438  vLocalPos += G4ThreeVector(vXtalSolid->GetXHalfLength(),
439  vXtalSolid->GetYHalfLength(),
440  vXtalSolid->GetZHalfLength());
441  }
442  else{
443  G4Tubs* vXtalSolid =
444  (G4Tubs*) vStepVol->GetPhysicalVolume()
445  ->GetLogicalVolume()->GetSolid();
446  vLocalPos.rotateZ(-vXtalSolid->GetStartPhiAngle());
447  }
448  return vLocalPos;
449 }
450 
451 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
452 
455  G4ThreeVector vPositionPre =
457  G4ThreeVector vMomentumPre =
458  GetXPhysicalLattice(aTrack)->
459  ProjectMomentumVectorFromWorldToLattice(
460  aTrack.GetStep()->GetPreStepPoint()->GetMomentum(),
461  vPositionPre);
462 
463  G4double vDeltaProportion = 1.;
464 
465  if((vMomentumPre.x())!=0.){
466  vDeltaProportion = std::fabs(vMomentumPre.unit().x());
467  }
468  G4double vDeltaPosition = vDeltaProportion*
470 
471  return std::abs(vDeltaPosition);
472 }
473 
474 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
475 
478 
479  G4StepPoint* vStepPre = aTrack.GetStep()->GetPreStepPoint();
480  G4StepPoint* vStepPost = aTrack.GetStep()->GetPostStepPoint();
481 
482  if( fLatticeManager->HasLattice(vStep->GetPhysicalVolume()) ) {
483  return vStep;
484  }
485  else if(fLatticeManager->HasLattice(vStepPost->GetPhysicalVolume())
486  == false
488  == true
489  && vStep == vStepPost &&
490  vStepPost->GetStepStatus() == fGeomBoundary) {
491  return vStepPre;
492  }
493  else if(fLatticeManager->HasLattice(vStepPre->GetPhysicalVolume())
494  == false &&
496  == true &&
497  vStep == vStepPre &&
498  vStepPre->GetStepStatus() == fGeomBoundary) {
499  return vStepPost;
500  }
501  else{
502  return vStep;
503  }
504 }
505 
506 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
507 
510 
511  G4StepPoint* vStepPre = aTrack.GetStep()->GetPreStepPoint();
512  G4StepPoint* vStepPost = aTrack.GetStep()->GetPostStepPoint();
513 
514  if( fLatticeManager->HasLattice(vStep->GetPhysicalVolume()) ) {
515  return vStep;
516  }
517  else if(fLatticeManager->HasLattice(vStepPost->GetPhysicalVolume())
518  == false &&
520  == true &&
521  vStep == vStepPost &&
522  vStepPost->GetStepStatus() == fGeomBoundary) {
523  return vStepPost;
524  }
525  else if(fLatticeManager->HasLattice(vStepPre->GetPhysicalVolume())
526  == false &&
528  == true &&
529  vStep == vStepPre &&
530  vStepPre->GetStepStatus() == fGeomBoundary) {
531  return vStepPre;
532  }
533  else{
534  return vStep;
535  }
536 }
537 
538 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
539 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
540 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
541 
544  //----------------------------------------
545  // check if the particle momentum
546  // transverse to the (h,k,l) plane
547  // is small enough to permit channeling
548  //----------------------------------------
549 
550  UpdateParameters(aTrack);
551 
552  G4double vEnergyMax = ComputeCriticalEnergyMaximum(aTrack);
553 
554  G4double vTransverseEnergy = ComputeTransverseEnergy(aTrack).x();
555 
556  if(GetXPhysicalLattice(aTrack)->IsBent() == false){
557  if(vTransverseEnergy <= vEnergyMax){
558  GetInfo(aTrack)->SetCoherentEffect(1);
559  // the particle is in channeling
560  return true;
561  }
562  }
563  else{
564  G4ThreeVector vPositionInTheCrystal =
565  GetInfo(aTrack)->GetPositionChanneled()
567  vTransverseEnergy += std::fabs(ComputeCentrifugalEnergy(aTrack,
568  vPositionInTheCrystal).x());
569  if(vTransverseEnergy <= vEnergyMax){
570  // the particle is in channeling
571  GetInfo(aTrack)->SetCoherentEffect(1);
572  return true;
573  }
574  else{
575  G4bool bNotBoundary = ParticleIsNotOnBoundary(aTrack);
576  G4bool bTangentToPlane = ParticleIsTangentToBentPlane(aTrack);
577 
578  if(bTangentToPlane == true &&
579  bNotBoundary == true &&
580  GetInfo(aTrack)->HasBeenUnderCoherentEffect() != 2){
581  // the particle is in volume reflection
582  GetInfo(aTrack)->SetCoherentEffect(2);
583  return true;
584  }
585  }
586  }
587 
588  // the particle is not under coherent effect
589  GetInfo(aTrack)->SetCoherentEffect(0);
590  return false;
591 }
592 
593 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
596  //----------------------------------------
597  // return the channeling MFP
598  //----------------------------------------
599 
600  G4double vMFPosc = ComputeOscillationPeriod(aTrack);
601 
602  if(GetInfo(aTrack)->GetNucleiDensity() < 1.){
603  vMFPosc /= GetInfo(aTrack)->GetNucleiDensity();
604  }
605 
606  G4double vMFP = vMFPosc * 2.;
607 
608  if(GetXPhysicalLattice(aTrack)->IsBent()){
609  G4double vMFPVR =
611 
612  if((std::fabs(vMFPVR) < vMFP) && (std::fabs(vMFPVR) > (0.5 * vMFPosc))){
613  vMFP = vMFPVR;
614  }
615  }
616 
617 
618  return vMFP;
619 }
620 
621 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
622 
624 GetMeanFreePath(const G4Track& aTrack,
625  G4double, // previousStepSize
627 
628  //----------------------------------------
629  // the condition is forced to check if
630  // the volume has a lattice at each step.
631  // if it hasn't, return DBL_MAX
632  //----------------------------------------
633 
634  *condition = Forced;
635 
636  if(HasLattice(aTrack)){
637  return GetChannelingMeanFreePath(aTrack);
638  }
639  else{
640  return DBL_MAX;
641  }
642 }
643 
644 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
645 
647 PostStepDoIt(const G4Track& aTrack,
648  const G4Step&){
649 
650  //----------------------------------------
651  // check if the volume has a lattice
652  // and if the particle is in channeling.
653  // If it is so, the particle is forced
654  // to follow the channeling plane
655  // direction. If the particle has
656  // dechanneled or exited the crystal,
657  // the outgoing angle is evaluated
658  //----------------------------------------
659 
660  aParticleChange.Initialize(aTrack);
661 
663 
664  G4bool bIsUnderCoherentEffect = false;
665 
666  if((HasLattice(aTrack) == true) &&
667  (ParticleIsNotOnBoundaryPost(aTrack) == true)){
668  bIsUnderCoherentEffect = IsUnderCoherentEffect(aTrack);
669  if(bIsUnderCoherentEffect){
670  if(GetInfo(aTrack)->HasBeenUnderCoherentEffect() == 1){
671  // if the particle is in channeling it gives the direction
672  // of the lattice to the particle momentum
673  G4ThreeVector vPosition =
675  aTrack.GetStep()->GetPostStepPoint(),aTrack);
677  GetXPhysicalLattice(aTrack)->
678  GetLatticeDirection(vPosition).unit());
679  }
680  else if(GetInfo(aTrack)->HasBeenUnderCoherentEffect() == 2){
681  // if the particle is in VR it gives a kick
682  // to the opposite site of the bending to the particle
685  GetInfo(aTrack)->SetCoherentEffect(0);
686  ResetDensity(aTrack);
687  }
688  }
689  }
690  else{
691  // if the volume has no lattice it resets the density factors
692  ResetDensity(aTrack);
693  }
694 
695  if( (bIsUnderCoherentEffect == false && (HasLattice(aTrack) == true) )
696  || (HasLatticeOnBoundaryPre(aTrack) == true) ) {
697  // if has been under coherent effect but now it is not,
698  // the outgoing momentum is evaluated starting from the current position
699  if(GetInfo(aTrack)->HasBeenUnderCoherentEffect() == 1){
702  }
703 
704  // If is not under coherent effect sets coherent effect to zero
705  // and resets the density factors after the outgoing angle
706  // has been evaluated
707  GetInfo(aTrack)->SetCoherentEffect(0);
708  ResetDensity(aTrack);
709  }
710 
711  return &aParticleChange;
712 }
713 
714 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
715 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
716 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
717 
720  //----------------------------------------
721  // compute the particle transverse energy
722  // in the crystal reference system
723  //----------------------------------------
724 
725  G4ThreeVector vTransverseEnergy = ComputePotentialEnergy(aTrack)
726  + ComputeKineticEnergy(aTrack);
727  //while(!getchar());
728  return vTransverseEnergy;
729 }
730 
733  //----------------------------------------
734  // compute the particle kinetic energy
735  // in the crystal reference system
736  //----------------------------------------
737 
738  G4double vTotalEnergy =
739  aTrack.GetStep()->GetPostStepPoint()->GetTotalEnergy();
740 
741  G4ThreeVector vMom = GetInfo(aTrack)->GetMomentumChanneled();
742 
743  G4ThreeVector vKineticEnergy = 0.5 *
744  G4ThreeVector((vMom.x() * vMom.x()) / vTotalEnergy,
745  0.,
746  (vMom.z() * vMom.z()) / vTotalEnergy);
747 
748  return vKineticEnergy;
749 }
750 
751 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
752 
755  //----------------------------------------
756  // compute the particle transverse energy
757  // in the crystal reference system
758  //----------------------------------------
759 
760 
761  G4ThreeVector vPotentialEnergy =
762  fPotentialEnergy->GetEC(GetInfo(aTrack)->GetPositionChanneledInitial(),
763  GetXPhysicalLattice(aTrack));
764 
765  vPotentialEnergy *= GetParticleDefinition(aTrack)->GetPDGCharge();
766 
767  return vPotentialEnergy;
768 }
769 
770 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
771 
774  //----------------------------------------
775  // compute the transverse energy variation
776  // in the crystal reference system
777  //----------------------------------------
778 
779  G4double vTotalEnergy =
780  aTrack.GetStep()->GetPreStepPoint()->GetTotalEnergy();
781 
782  G4double vPositionX = vPosition.x();
783 
784  G4ThreeVector vEnergyVariation =
785  G4ThreeVector(vTotalEnergy * vPositionX /
786  GetXPhysicalLattice(aTrack)->GetCurvatureRadius().x(),
787  0.,
788  0.);
789 
790  return vEnergyVariation;
791 }
792 
793 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
794 
796 ComputeMomentum(const G4Track& aTrack,G4StepPoint* vStep){
797  //----------------------------------------
798  // compute the particle momentum
799  // in the crystal reference system
800  //----------------------------------------
801 
802  G4ThreeVector vPosition = ComputePositionInTheCrystal(vStep,aTrack);
803 
804  G4ThreeVector vMomentum =
805  GetXPhysicalLattice(aTrack)->
806  ProjectMomentumVectorFromWorldToLattice(aTrack.GetMomentum(),
807  vPosition);
808 
809  return vMomentum;
810 }
811 
812 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
813 
816  //----------------------------------------
817  // compute the transverse energy variation
818  // in the crystal reference system
819  //----------------------------------------
820 
821  G4double vTotalEnergy =
822  aTrack.GetStep()->GetPreStepPoint()->GetTotalEnergy();
823 
824  G4ThreeVector vEnergyVariation = G4ThreeVector(vTotalEnergy *
825  GetXPhysicalLattice(aTrack)->ComputeInterplanarPeriod() /
826  GetXPhysicalLattice(aTrack)->GetCurvatureRadius().x(),
827  0.,
828  0.);
829 
830  return vEnergyVariation;
831 }
832 
833 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
834 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
835 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
836 
837 
840  //----------------------------------------
841  // compute the critical energy
842  // for channeling
843  //----------------------------------------
844 
845  G4double vCriticalEnergy = 0.;
846 
847  if(ParticleIsNegative(aTrack)){
848  vCriticalEnergy =
850  }
851  else{
852  vCriticalEnergy =
854  }
855 
856  vCriticalEnergy *= std::fabs(GetParticleDefinition(aTrack)->GetPDGCharge());
857 
858  return vCriticalEnergy;
859 }
860 
861 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
862 
865  //----------------------------------------
866  // compute the critical energy minimum
867  // for channeling
868  //----------------------------------------
869 
870  G4double vCriticalEnergy = 0.;
871 
872  if(ParticleIsNegative(aTrack)){
873  vCriticalEnergy =
875  }
876  else{
877  vCriticalEnergy =
879  }
880 
881  vCriticalEnergy *= std::fabs(GetParticleDefinition(aTrack)->GetPDGCharge());
882 
883  return vCriticalEnergy;
884 }
885 
886 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
887 
890  //----------------------------------------
891  // compute the critical angle
892  // for chenneling
893  //----------------------------------------
894 
895  G4double vTotalEnergy =
896  aTrack.GetStep()->GetPostStepPoint()->GetTotalEnergy();
897  G4double vCriticalAngle = std::pow( 2.0 *
898  std::fabs( ( ComputeCriticalEnergyMaximum(aTrack)
899  - ComputeCriticalEnergyMinimum(aTrack) )
900  / vTotalEnergy ) , 0.5);
901  return vCriticalAngle;
902 }
903 
904 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
905 
908  //----------------------------------------
909  // compute the particle oscillation
910  // period in the crystal channel
911  //----------------------------------------
912 
913  G4double vInterplanarPeriod =
915  G4double vOscillationPeriod =
916  CLHEP::pi * vInterplanarPeriod / ComputeCriticalAngle(aTrack);
917  return vOscillationPeriod;
918 }
919 
920 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
921 
924  //----------------------------------------
925  // compute the critical radius
926  // for channeling
927  //----------------------------------------
928 
929  G4double vTotalEnergy =
930  aTrack.GetStep()->GetPreStepPoint()->GetTotalEnergy();
931  G4double vCriticalRadius =
932  vTotalEnergy / fElectricField->GetMaximum(GetXPhysicalLattice(aTrack));
933  return vCriticalRadius;
934 }
935 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
936 
939  //----------------------------------------
940  // compute the central point of
941  // the potential well for channeling
942  //----------------------------------------
943 
944  G4double vInterplanarPeriodHalf =
946 
947  G4double vCentreX = vInterplanarPeriodHalf;
948 
949  if(GetXPhysicalLattice(aTrack)->IsBent()){
950  G4double vTotalEnergy =
951  aTrack.GetStep()->GetPreStepPoint()->GetTotalEnergy();
952 
953  G4double vPotentialWellDepth =
955  - (ComputeCriticalEnergyMinimum(aTrack));
956 
957  vCentreX *= (1. - 0.5 * vTotalEnergy /
958  vPotentialWellDepth /
959  GetXPhysicalLattice(aTrack)->GetCurvatureRadius().x() *
960  vInterplanarPeriodHalf );
961  }
962 
963  G4ThreeVector vCentre = G4ThreeVector(vCentreX,0.,0.);
964 
965 
966  return vCentre;
967 }
968 
969 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
970 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
971 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
972 
973 
976  return(aPD.GetPDGCharge() != 0.);
977 }
978 
979 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
980 
983 
984 }
985 
986 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
987 
991  aTrack.GetStep()->GetPostStepPoint()->GetPhysicalVolume())){
994  }
996  ->GetPhysicalVolume()) &&
997  aTrack.GetStep()->GetPostStepPoint()->GetStepStatus()
998  == fGeomBoundary) {
1001  }
1002  else{
1003  G4cout << "LATTICE NOT FOUND: ERROR" << G4endl;
1004  return NULL;
1005  }
1006 }
1007 
1008 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1009 
1011  if(fLatticeManager->HasLattice(aTrack.GetStep()
1013  return true;
1014  }
1015  else{
1016  return false;
1017  }
1018 }
1019 
1020 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1021 
1024  if(HasLatticeOnBoundaryPost(aTrack) || HasLatticeOnBoundaryPre(aTrack)) {
1025  return true;
1026  }
1027  else{
1028  return false;
1029  }
1030 }
1031 
1032 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1033 
1037  GetPhysicalVolume()) &&
1039  return true;
1040  }
1041  else{
1042  return false;
1043  }
1044 }
1045 
1046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1047 
1050  if(fLatticeManager->HasLattice(aTrack.GetStep()->
1051  GetPostStepPoint()->GetPhysicalVolume()) &&
1052  aTrack.GetStep()->GetPreStepPoint()->GetStepStatus() == fGeomBoundary) {
1053  return true;
1054  }
1055  else{
1056  return false;
1057  }
1058 }
1059 
1060 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1061 
1063  if(GetParticleDefinition(aTrack)->GetPDGCharge() < 0.) {
1064  return true;
1065  }
1066  else{
1067  return false;
1068  }
1069 }
1070 
1071 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1072 
1075  G4ThreeVector vPositionPre =
1077  G4ThreeVector vMomentumPre =
1078  GetXPhysicalLattice(aTrack)->
1079  ProjectMomentumVectorFromWorldToLattice(aTrack.GetStep()->
1080  GetPreStepPoint()->GetMomentum(),vPositionPre);
1081 
1082  G4ThreeVector vPositionPost =
1084  GetPostStepPoint(),aTrack);
1085  G4ThreeVector vMomentumPost = GetXPhysicalLattice(aTrack)->
1086  ProjectMomentumVectorFromWorldToLattice(aTrack.GetStep()->
1087  GetPostStepPoint()->GetMomentum(),vPositionPost);
1088 
1089  if(vMomentumPost.x()<0. &&
1090  vMomentumPre.x()>0. &&
1091  GetXPhysicalLattice(aTrack)->GetCurvatureRadius().x() < 0.){
1092  return true;
1093  }
1094  if(vMomentumPost.x()>0. &&
1095  vMomentumPre.x()<0. &&
1096  GetXPhysicalLattice(aTrack)->GetCurvatureRadius().x() > 0.){
1097  return true;
1098  }
1099 
1100  return false;
1101 }
1102 
1103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1104 
1107  if(aTrack.GetStep()->GetPreStepPoint()->GetStepStatus() != fGeomBoundary){
1108  return true;
1109  }
1110  else{
1111  return false;
1112  }
1113 }
1114 
1115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1116 
1119  if(aTrack.GetStep()->GetPostStepPoint()->GetStepStatus() != fGeomBoundary){
1120  return true;
1121  }
1122  else{
1123  return false;
1124  }
1125 }
1126 
1127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1128 
1131  if(ParticleIsNotOnBoundaryPost(aTrack) ||
1132  ParticleIsNotOnBoundaryPre(aTrack)){
1133  return true;
1134  }
1135  else{
1136  return false;
1137  }
1138 }
1139 
1140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1141 
1143 GetInfo(const G4Track& aTrack){
1144  return (ExExChParticleUserInfo*) aTrack.GetUserInformation();
1145 }
1146 
1147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1148 
1151  return const_cast<G4ParticleDefinition*>(aTrack.GetParticleDefinition());
1152 }
1153 
1154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1155 
1160 
1161  if(fFileCharacteristicsName != ""){
1162  G4String filename;
1163 
1164  fPotentialEnergy->ReadFromECHARM(filename =
1165  fFileCharacteristicsName + "_pot.txt");
1166  //fElectricField->ReadFromECHARM("efx.txt");
1169  }
1170  else{
1172  GetXPhysicalLattice(aTrack));
1174  }
1175 }
1176 
1177 //....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
G4double ComputeDistanceWhereParticleTangentToBentPlane(const G4Track &)
void SetMomentumChanneledInitial(G4ThreeVector)
G4ThreeVector GetLatticeAngles()
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector ComputeChannelingOutgoingMomentum(const G4Track &)
XVCrystalCharacteristic * GetNucleiDensity()
Definition: G4Box.hh:64
G4bool ParticleIsNegative(const G4Track &)
G4StepStatus GetStepStatus() const
const G4double pi
Definition: G4Tubs.hh:85
G4double GetSurfaceTolerance() const
G4ThreeVector ComputeMomentum(const G4Track &, G4StepPoint *)
G4bool IsUnderCoherentEffect(const G4Track &)
ExExChProcessChanneling(const G4String &processName="ExExChProcessChanneling")
G4ThreeVector GetCurvatureRadius()
G4ThreeVector ComputeKineticEnergy(const G4Track &)
G4ThreeVector GetMomentum() const
void UpdateDensity(const G4Track &)
const G4VTouchable * GetTouchable() const
G4ThreeVector GetPositionChanneledInitial()
const G4Step * GetStep() const
int G4int
Definition: G4Types.hh:78
G4double GetZHalfLength() const
XVCrystalIntegratedDensity * GetIntegratedDensityNuclei(G4int)
XVCrystalIntegratedDensity * GetIntegratedDensityElectron(G4int)
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 *)
G4ThreeVector ComputeVolumeReflectionOutgoingMomentum(const G4Track &)
G4double ComputeCriticalEnergyMaximum(const G4Track &)
void SetMomentumChanneled(G4ThreeVector)
G4bool ParticleIsNotOnBoundary(const G4Track &)
#define G4UniformRand()
Definition: Randomize.hh:93
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
XPhysicalLattice * GetXPhysicalLattice(const G4Track &)
void SetPotential(XVCrystalCharacteristic *)
const G4ParticleDefinition * GetParticleDefinition() const
G4double GetYHalfLength() const
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 *)
G4ThreeVector ComputeTransverseEnergy(const G4Track &)
G4ParticleDefinition * GetParticleDefinition(const G4Track &aTrack)
G4ThreeVector GetMomentum() const
G4bool ParticleIsNotOnBoundaryPost(const G4Track &)
virtual void Initialize(const G4Track &)
G4ThreeVector ComputePotentialWellCentre(const G4Track &)
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
void SetElectricField(XVCrystalCharacteristic *)
void UpdateParameters(const G4Track &)
void InitializePhysicalLattice(XPhysicalLattice *)
G4LogicalVolume * GetLogicalVolume() const
XPhysicalLattice * GetXPhysicalLattice(G4VPhysicalVolume *)
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 *)
G4bool HasLatticeOnBoundaryPost(const G4Track &)
void SetCoherentEffect(G4int flag)
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
ExExChParticleUserInfo * GetInfo(const G4Track &)
void SetElectronDensity(XVCrystalCharacteristic *)
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
const G4NavigationHistory * GetHistory() const
void ComputeCrystalCharacteristic(const G4Track &)
void UpdateMomentum(const G4Track &)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4bool ParticleIsTangentToBentPlane(const G4Track &)
double G4double
Definition: G4Types.hh:76
G4ForceCondition
G4double GetPDGCharge() const
#define DBL_MAX
Definition: templates.hh:83
static const double mm
Definition: G4SIunits.hh:102
void ResetDensity(const G4Track &)
G4double ComputeInterplanarPeriod()
void SetPositionChanneledInitial(G4ThreeVector)
G4bool ParticleIsNotOnBoundaryPre(const G4Track &)
virtual G4bool IsApplicable(const G4ParticleDefinition &)
void UpdatePosition(const G4Track &)
XVCrystalCharacteristic * fPotentialEnergy
static G4GeometryTolerance * GetInstance()
G4VSolid * GetSolid() const
virtual void ReadFromECHARM(const G4String &, G4double=1)=0
bool HasLattice(G4VPhysicalVolume *)