Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4Channeling.cc
이 파일의 문서화 페이지로 가기
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 
27 #include "G4Channeling.hh"
28 
29 #include "Randomize.hh"
30 
31 #include "G4ChannelingTrackData.hh"
32 #include "G4TouchableHistory.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4LambdacPlus.hh"
35 
37 G4VDiscreteProcess("channeling"),
38 fChannelingID(-1),
39 fTimeStepMin(0.),
40 fTimeStepMax(0.),
41 fTransverseVariationMax(2.E-2 * CLHEP::angstrom),
42 k010(G4ThreeVector(0.,1.,0.)){
44  if(fChannelingID == -1){
46  }
47  fSpin = G4ThreeVector(0.,0.,0.);
48 }
49 
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 
53 
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55 
57  G4ChannelingTrackData* trackdata =
59  if(trackdata == nullptr){
60  trackdata = new G4ChannelingTrackData();
62  }
63  return trackdata;
64 }
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67 
68 void G4Channeling::GetEF(const G4Track& aTrack,
70  G4ThreeVector& out){
71  out = G4ThreeVector((GetMatData(aTrack)->GetEFX()->GetEC(pos)),
72  (GetMatData(aTrack)->GetEFY()->GetEC(pos)),
73  0.);
74 }
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 
79  G4TouchableHistory* theTouchable = (G4TouchableHistory*)(step->GetTouchable());
80 
81  pos -= theTouchable->GetTranslation();
82  pos = ((*theTouchable->GetRotation()).inverse())(pos);
83 }
84 
85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86 
88 
90 
91  G4StepPoint* postStepPoint = aTrack.GetStep()->GetPostStepPoint();
92  G4StepPoint* preStepPoint = aTrack.GetStep()->GetPreStepPoint();
93 
94  G4ThreeVector posPost = postStepPoint->GetPosition();
95  aLCV->RotateToLattice(posPost);
96  G4ThreeVector posPre = preStepPoint->GetPosition();
97  aLCV->RotateToLattice(posPre);
98 
99  G4double integrationLimit = std::fabs(posPost.z() - posPre.z());
100 
101  if(integrationLimit > 0.){
102  //----------------------------------------
103  // Check if the crystal is bent
104  //----------------------------------------
105  G4bool isBent = GetMatData(aTrack)->IsBent();
106 
107  //----------------------------------------
108  // Get the momentum in the world reference
109  // frame and rotate to the solid reference frame
110  //----------------------------------------
111  G4TouchableHistory* theTouchable = (G4TouchableHistory*)(preStepPoint->GetTouchable());
112  G4ThreeVector momWorld = aTrack.GetStep()->GetPreStepPoint()->GetMomentum();
113  G4ThreeVector mom = (*theTouchable->GetRotation())(momWorld);
114 
115  //----------------------------------------
116  // Get the momentum in the solid reference
117  // frame and rotate to the crystal reference frame
118  //----------------------------------------
119  aLCV->RotateToLattice(mom);
120 
121  //----------------------------------------
122  // Get the momentum in the crystal reference
123  // frame and rotate to the reference frame
124  // solidal to the bent planes
125  //----------------------------------------
126  if(isBent){
127  PosToLattice(preStepPoint,posPre);
128  G4ThreeVector axis010 = (*theTouchable->GetRotation())(k010);
129  mom.rotate(axis010,-posPre.z()/GetMatData(aTrack)->GetBR(posPre).x());
130  }
131 
132  //----------------------------------------
133  // Take the position stored in the track data.
134  // If the particle enters the crystal,
135  // the position in the channel is randomly
136  // generated using a uniform distribution
137  //----------------------------------------
139  if(GetTrackData(aTrack)->GetPosCh().x() == DBL_MAX){
140  G4double posX = G4UniformRand() * GetMatData(aTrack)->GetPot()->GetIntSp(0);
141  G4double posY = G4UniformRand() * GetMatData(aTrack)->GetPot()->GetIntSp(1);
142  pos = G4ThreeVector(posX,posY,0.);
143  }
144  else{
145  pos = GetTrackData(aTrack)->GetPosCh();
146  }
147 
148  G4double step=0., stepTot=0.;
149  G4double nud =0., eld =0.;
150  G4double efx =0., efy =0.;
151  G4double nud_temp =0., eld_temp =0.;
152 
155 
156  const G4double oneSixth = 1./6.;
157  G4ThreeVector posk1,posk2,posk3,posk4,posk5,posk6;
158  G4ThreeVector momk1,momk2,momk3,momk4,momk5,momk6;
159  G4ThreeVector pos_temp, efxy;
160 
161  do{
162  //----------------------------------------
163  // Limit the variable step length for the
164  // integration via the selected algorithm
165  // and update variables for the integration
166  //----------------------------------------
167 
168  UpdateIntegrationStep(aTrack,mom,step);
169  if(step + stepTot > integrationLimit){
170  step = integrationLimit - stepTot;
171  }
172 
173  //----------------------------------------
174  // Function integration algorithm
175  // 4th Order Runge-Kutta
176  //----------------------------------------
177 
178  GetEF(aTrack,pos,efxy);
179  posk1 = step / mom.z() * mom;
180  momk1 = step / beta * Z * efxy;
181  if(isBent) momk1.setX(momk1.x() - step * mom.z() * beta / (GetMatData(aTrack)->GetBR(pos)).x());
182 
183  GetEF(aTrack,pos_temp = pos + posk1 * 0.5,efxy);
184  posk2 = step / mom.z() * (mom + momk1 * 0.5);
185  momk2 = step / beta * Z * efxy;
186  if(isBent) momk2.setX(momk2.x() - step * mom.z() * beta / (GetMatData(aTrack)->GetBR(pos_temp)).x());
187 
188  GetEF(aTrack,pos_temp = pos + posk2 * 0.5,efxy);
189  posk3 = step / mom.z() * (mom + momk2 * 0.5);
190  momk3 = step / beta * Z * efxy;
191  if(isBent) momk3.setX(momk3.x() - step * mom.z() * beta / (GetMatData(aTrack)->GetBR(pos_temp)).x());
192 
193  GetEF(aTrack,pos_temp = pos + posk3,efxy);
194  posk4 = step / mom.z() * (mom + momk3);
195  momk4 = step / beta * Z * efxy;
196  if(isBent) momk4.setX(momk4.x() - step * mom.z() * beta / (GetMatData(aTrack)->GetBR(pos_temp)).x());
197 
198  pos = pos + oneSixth * (posk1 + 2.*posk2 + 2.*posk3 + posk4);
199  mom = mom + oneSixth * (momk1 + 2.*momk2 + 2.*momk3 + momk4);
200 
201  //----------------------------------------
202  // Function integration algorithm
203  // 2th Order Velocity-Verlet
204  //----------------------------------------
205 
206  /*
207  GetEF(aTrack,pos,efxy);
208  posk1 = pos + (step * 0.5 / mom.z()) * mom;
209  //momk1 = mom + step * 0.5 / betaZ * efxy;
210  momk1 = mom;
211  if(isBent) momk1.setX(momk1.x() - step * 0.5 * mom.z() * beta / (GetMatData(aTrack)->GetBR(pos)).x());
212 
213  GetEF(aTrack,posk1,efxy);
214  pos = pos + (step / momk1.z()) * momk1;
215  //mom = mom + step / betaZ * efxy;
216  mom = mom;
217  if(isBent) mom.setX(mom.x() - step * mom.z() * beta / (GetMatData(aTrack)->GetBR(posk1)).x());
218  */
219 
220  //----------------------------------------
221  // Update the total step and the electron
222  // and nuclei density experienced by
223  // the particle during its motion
224  //----------------------------------------
225 
226  stepTot += step;
227 
228  nud_temp = GetMatData(aTrack)->GetNuD()->GetEC(pos);
229  eld_temp = GetMatData(aTrack)->GetElD()->GetEC(pos);
230 
231  if(nud_temp < 0.) {nud_temp = 0.;}
232  if(eld_temp < 0.) {eld_temp = 0.;}
233 
234  nud += (step * nud_temp);
235  eld += (step * eld_temp);
236 
237  efx += (step * GetMatData(aTrack)->GetEFX()->GetEC(pos));
238  efy += (step * GetMatData(aTrack)->GetEFY()->GetEC(pos));
239 
240  } while(stepTot<integrationLimit);
241 
242  nud /= stepTot;
243  eld /= stepTot;
244 
245  if(nud < 1.E-10) {nud = 1.E-10;}
246  if(eld < 1.E-10) {eld = 1.E-10;}
247 
248  GetTrackData(aTrack)->SetNuD(nud);
249  GetTrackData(aTrack)->SetElD(eld);
250 
251  GetTrackData(aTrack)->SetEFX(efx);
252  GetTrackData(aTrack)->SetEFY(efy);
253 
254  GetTrackData(aTrack)->SetMomCh(mom);
255  GetTrackData(aTrack)->SetPosCh(pos);
256  return true;
257  }
258  else{
259  return false;
260  }
261 
262  return false;
263 }
264 
265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
266 
269  G4ThreeVector& mom,
270  G4double& step){
271 
272  if(mom.x() != 0.0 || mom.y() != 0.0){
273  double xy2 = mom.x() * mom.x() + mom.y()*mom.y();
274 
275  if(xy2!=0.){
276  step = std::fabs(fTransverseVariationMax * GetPre(aTrack)->GetKineticEnergy() / std::pow(xy2,0.5));
277  if(step < fTimeStepMin) step = fTimeStepMin;
278  else{
279  fTimeStepMax = std::sqrt( fTransverseVariationMax * GetPre(aTrack)->GetKineticEnergy()
280  / std::fabs(GetMatData(aTrack)->GetEFX()->GetMax()));
281 
282  if(step > fTimeStepMax) step = fTimeStepMax;
283  }
284  }
285  else{
286  step = fTimeStepMin;
287  }
288 
289  return true;
290  }
291  else{
292  step = fTimeStepMin;
293  }
294  return false;
295 }
296 
297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
298 
300 GetMeanFreePath(const G4Track& aTrack,
301  G4double, // previousStepSize
303 
304  //----------------------------------------
305  // the condition is forced to check if
306  // the volume has a lattice at each step.
307  // if it hasn't, return DBL_MAX
308  //----------------------------------------
309 
310  *condition = Forced;
311 
312  G4LogicalVolume* aLV = aTrack.GetVolume()->GetLogicalVolume();
313  G4LogicalVolume* aNLV = aTrack.GetNextVolume()->GetLogicalVolume();
314 
315  if(G4LogicalCrystalVolume::IsLattice(aLV) == true &&
316  G4LogicalCrystalVolume::IsLattice(aNLV) == true){
317  G4double osc_per = GetOscillationPeriod(aTrack);
318  fTimeStepMin = osc_per * 2.E-4;
319  return osc_per * 0.01;
320  }
321  else{
322  GetTrackData(aTrack)->Reset();
323  return DBL_MAX;
324  }
325 }
326 
327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
328 
330 PostStepDoIt(const G4Track& aTrack,
331  const G4Step&){
332 
333  //----------------------------------------
334  // check if the volume has a lattice
335  // and if the particle is in channeling.
336  // If it is so, the particle is forced
337  // to follow the channeling plane
338  // direction. If the particle has
339  // dechanneled or exited the crystal,
340  // the outgoing angle is evaluated
341  //----------------------------------------
342 
343  aParticleChange.Initialize(aTrack);
344  G4LogicalVolume* aLV = aTrack.GetVolume()->GetLogicalVolume();
345  G4LogicalVolume* aNLV = aTrack.GetNextVolume()->GetLogicalVolume();
346 
347 
348  if(G4LogicalCrystalVolume::IsLattice(aLV) == true &&
349  G4LogicalCrystalVolume::IsLattice(aNLV) == true){
350 
351  G4bool bModifiedTraj = UpdateParameters(aTrack);
352 
353  if(bModifiedTraj==true){
354  //----------------------------------------
355  // Get the momentum in the reference frame
356  // solidal to the bent planes and rotate
357  // to the reference frame
358  //----------------------------------------
360  G4ThreeVector momCh = GetTrackData(aTrack)->GetMomCh();
361 
362  G4StepPoint* postStepPoint = aTrack.GetStep()->GetPostStepPoint();
363  G4TouchableHistory* theTouchable = (G4TouchableHistory*)(postStepPoint->GetTouchable());
364 
365  if(GetMatData(aTrack)->IsBent()){
366  G4ThreeVector posPost = postStepPoint->GetPosition();
367  PosToLattice(postStepPoint,posPost);
368  G4ThreeVector axis010 = (*theTouchable->GetRotation())(k010);
369  momCh.rotate(axis010,posPost.z()/GetMatData(aTrack)->GetBR(posPost).x());
370  }
371 
372  //----------------------------------------
373  // Get the momentum in the crystal reference
374  // frame and rotate to the solid reference frame
375  //----------------------------------------
376 
377  aLCV->RotateToSolid(momCh);
378 
379  //----------------------------------------
380  // Get the momentum in the solid reference
381  // frame and rotate to the world reference frame
382  //----------------------------------------
383  G4ThreeVector mom = ((*theTouchable->GetRotation()).inverse())(momCh);
384 
387  }
388  }
389  else{
390  // if the volume has no lattice it resets the density factors
391  GetTrackData(aTrack)->Reset();
392  }
393 
394  return &aParticleChange;
395 }
396 
397 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Float_t x
Definition: compare.C:6
static G4bool IsLattice(G4LogicalVolume *aLV)
const G4ThreeVector & GetTranslation(G4int depth=0) const
G4double GetIntSp(G4int index)
G4double fTimeStepMin
const G4VTouchable * GetTouchable() const
CLHEP::Hep3Vector G4ThreeVector
void ProposePolarization(G4double Px, G4double Py, G4double Pz)
G4double GetMax(const G4ToolsBaseHisto &baseHisto, G4int dimension)
static const G4double pos
G4LogicalVolume * GetLogicalVolume() const
const G4ThreeVector k010
const G4ThreeVector & RotateToLattice(G4ThreeVector &dir) const
G4StepPoint * GetPreStepPoint() const
void PosToLattice(G4StepPoint *step, G4ThreeVector &)
Definition: G4Channeling.cc:78
G4ParticleDefinition * GetParticleDefinition(const G4Track &aTrack)
Definition: G4Channeling.hh:59
Float_t posY
Definition: plotDend.C:21
void GetEF(const G4Track &, G4ThreeVector &, G4ThreeVector &)
Definition: G4Channeling.cc:68
G4double GetPDGCharge() const
double z() const
G4ThreeVector fSpin
void SetMomCh(G4ThreeVector a3vec)
Double_t beta
G4double condition(const G4ErrorSymMatrix &m)
G4VPhysicalVolume * GetNextVolume() const
void setX(double)
G4double GetEC(G4ThreeVector &)
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
void SetEFX(G4double aDouble)
static G4int GetIndex(const G4String &)
void SetPosCh(G4ThreeVector a3vec)
Float_t Z
G4double GetOscillationPeriod(const G4Track &aTrack)
Definition: G4Channeling.hh:87
G4double fTimeStepMax
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
virtual G4double GetMeanFreePath(const G4Track &, G4double, G4ForceCondition *)
G4ChannelingTrackData * GetTrackData(const G4Track &)
Definition: G4Channeling.cc:56
const G4RotationMatrix * GetRotation(G4int depth=0) const
virtual G4ThreeVector GetBR(G4ThreeVector &v3)
virtual void Initialize(const G4Track &)
const G4ThreeVector & GetPosition() const
#define G4UniformRand()
Definition: Randomize.hh:53
virtual ~G4Channeling()
Definition: G4Channeling.cc:52
Definition: G4Step.hh:76
G4StepPoint * GetPostStepPoint() const
const G4Step * GetStep() const
G4double GetVelocity() const
Hep3Vector & rotate(double, const Hep3Vector &)
Definition: ThreeVectorR.cc:28
Hep3Vector unit() const
G4StepPoint * GetPre(const G4Track &aTrack)
Definition: G4Channeling.hh:63
void SetNuD(G4double aDouble)
G4VAuxiliaryTrackInformation * GetAuxiliaryTrackInformation(G4int idx) const
Definition: G4Track.cc:348
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
void SetEFY(G4double aDouble)
static constexpr double c_light
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void SetElD(G4double aDouble)
G4ForceCondition
double x() const
G4ChannelingMaterialData * GetMatData(const G4Track &aTrack)
Definition: G4Channeling.hh:68
void SetAuxiliaryTrackInformation(G4int idx, G4VAuxiliaryTrackInformation *info) const
Definition: G4Track.cc:332
static constexpr double angstrom
Definition: G4SIunits.hh:102
G4int fChannelingID
Definition: G4Channeling.hh:95
G4ThreeVector GetMomentum() const
double y() const
Float_t posX
Definition: plotDend.C:21
#define DBL_MAX
Definition: templates.hh:83
G4double fTransverseVariationMax
G4bool UpdateParameters(const G4Track &)
Definition: G4Channeling.cc:87
G4bool UpdateIntegrationStep(const G4Track &, G4ThreeVector &, G4double &)
const G4ThreeVector & RotateToSolid(G4ThreeVector &dir) const
G4VPhysicalVolume * GetVolume() const
static G4int Register(const G4String &)