Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4VXTRenergyLoss.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 // $Id: G4VXTRenergyLoss.cc 108258 2018-01-26 13:49:05Z gcosmo $
28 //
29 // History:
30 // 2001-2002 R&D by V.Grichine
31 // 19.06.03 V. Grichine, modifications in BuildTable for the integration
32 // in respect of angle: range is increased, accuracy is
33 // improved
34 // 28.07.05, P.Gumplinger add G4ProcessType to constructor
35 // 28.09.07, V.Ivanchenko general cleanup without change of algorithms
36 //
37 
38 #include "G4VXTRenergyLoss.hh"
39 
40 #include "G4Timer.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 
44 #include "G4Poisson.hh"
45 #include "G4MaterialTable.hh"
46 #include "G4VDiscreteProcess.hh"
47 #include "G4VParticleChange.hh"
48 #include "G4VSolid.hh"
49 
50 #include "G4RotationMatrix.hh"
51 #include "G4ThreeVector.hh"
52 #include "G4AffineTransform.hh"
53 #include "G4SandiaTable.hh"
54 
55 #include "G4PhysicsVector.hh"
56 #include "G4PhysicsFreeVector.hh"
57 #include "G4PhysicsLinearVector.hh"
58 #include "G4EmProcessSubType.hh"
59 
61 //
62 // Constructor, destructor
63 
65  G4Material* foilMat,G4Material* gasMat,
66  G4double a, G4double b,
67  G4int n,const G4String& processName,
68  G4ProcessType type) :
69  G4VDiscreteProcess(processName, type),
70  fGammaCutInKineticEnergy(nullptr),
71  fGammaTkinCut(0.0),
72  fAngleDistrTable(nullptr),
73  fEnergyDistrTable(nullptr),
74  fPlatePhotoAbsCof(nullptr),
75  fGasPhotoAbsCof(nullptr),
76  fAngleForEnergyTable(nullptr)
77 {
78  verboseLevel = 1;
80 
81  fPtrGamma = nullptr;
84 
85  // Initialization of local constants
86  fTheMinEnergyTR = 1.0*keV;
87  fTheMaxEnergyTR = 100.0*keV;
88 
89  fTheMaxAngle = 1.0e-2; // 100 mrad
90 
91  fTheMinAngle = 2.5e-5;
92 
93  fBinTR = 200; // 100; // 50;
94 
95  fMinProtonTkin = 100.0*GeV;
96  fMaxProtonTkin = 100.0*TeV;
97  fTotBin = 50; // 100; //
98 
99  // Proton energy vector initialization
102  fTotBin );
103 
106  fBinTR );
107 
109 
111 
112  fEnvelope = anEnvelope;
113 
114  fPlateNumber = n;
115  if(verboseLevel > 0)
116  G4cout<<"### G4VXTRenergyLoss: the number of TR radiator plates = "
117  <<fPlateNumber<<G4endl;
118  if(fPlateNumber == 0)
119  {
120  G4Exception("G4VXTRenergyLoss::G4VXTRenergyLoss()","VXTRELoss01",
121  FatalException,"No plates in X-ray TR radiator");
122  }
123  // default is XTR dEdx, not flux after radiator
124 
125  fExitFlux = false;
126  fAngleRadDistr = true; // false;
127  fCompton = false;
128 
129  fLambda = DBL_MAX;
130  // Mean thicknesses of plates and gas gaps
131 
132  fPlateThick = a;
133  fGasThick = b;
135  if(verboseLevel > 0)
136  G4cout<<"total radiator thickness = "<<fTotalDist/cm<<" cm"<<G4endl;
137 
138  // index of plate material
139  fMatIndex1 = foilMat->GetIndex();
140  if(verboseLevel > 0)
141  G4cout<<"plate material = "<<foilMat->GetName()<<G4endl;
142 
143  // index of gas material
144  fMatIndex2 = gasMat->GetIndex();
145  if(verboseLevel > 0)
146  G4cout<<"gas material = "<<gasMat->GetName()<<G4endl;
147 
148  // plasma energy squared for plate material
149 
151  // fSigma1 = (20.9*eV)*(20.9*eV);
152  if(verboseLevel > 0)
153  G4cout<<"plate plasma energy = "<<std::sqrt(fSigma1)/eV<<" eV"<<G4endl;
154 
155  // plasma energy squared for gas material
156 
158  if(verboseLevel > 0)
159  G4cout<<"gas plasma energy = "<<std::sqrt(fSigma2)/eV<<" eV"<<G4endl;
160 
161  // Compute cofs for preparation of linear photo absorption
162 
165 
167 }
168 
170 
172 {
173  if(fEnvelope) delete fEnvelope;
174  delete fProtonEnergyVector;
175  delete fXTREnergyVector;
176  if(fEnergyDistrTable) {
178  delete fEnergyDistrTable;
179  }
180  if(fAngleRadDistr) {
182  delete fAngleDistrTable;
183  }
186  delete fAngleForEnergyTable;
187  }
188 }
189 
191 //
192 // Returns condition for application of the model depending on particle type
193 
194 
196 {
197  return ( particle.GetPDGCharge() != 0.0 );
198 }
199 
201 //
202 // Calculate step size for XTR process inside raaditor
203 
205  G4double, // previousStepSize,
207 {
208  G4int iTkin, iPlace;
209  G4double lambda, sigma, kinEnergy, mass, gamma;
210  G4double charge, chargeSq, massRatio, TkinScaled;
211  G4double E1,E2,W,W1,W2;
212 
213  *condition = NotForced;
214 
215  if( aTrack.GetVolume()->GetLogicalVolume() != fEnvelope ) lambda = DBL_MAX;
216  else
217  {
218  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
219  kinEnergy = aParticle->GetKineticEnergy();
220  mass = aParticle->GetDefinition()->GetPDGMass();
221  gamma = 1.0 + kinEnergy/mass;
222  if(verboseLevel > 1)
223  {
224  G4cout<<" gamma = "<<gamma<<"; fGamma = "<<fGamma<<G4endl;
225  }
226 
227  if ( std::fabs( gamma - fGamma ) < 0.05*gamma ) lambda = fLambda;
228  else
229  {
230  charge = aParticle->GetDefinition()->GetPDGCharge();
231  chargeSq = charge*charge;
232  massRatio = proton_mass_c2/mass;
233  TkinScaled = kinEnergy*massRatio;
234 
235  for(iTkin = 0; iTkin < fTotBin; iTkin++)
236  {
237  if( TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
238  }
239  iPlace = iTkin - 1;
240 
241  if(iTkin == 0) lambda = DBL_MAX; // Tkin is too small, neglect of TR photon generation
242  else // general case: Tkin between two vectors of the material
243  {
244  if(iTkin == fTotBin)
245  {
246  sigma = (*(*fEnergyDistrTable)(iPlace))(0)*chargeSq;
247  }
248  else
249  {
250  E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1);
252  W = 1.0/(E2 - E1);
253  W1 = (E2 - TkinScaled)*W;
254  W2 = (TkinScaled - E1)*W;
255  sigma = ( (*(*fEnergyDistrTable)(iPlace ))(0)*W1 +
256  (*(*fEnergyDistrTable)(iPlace+1))(0)*W2 )*chargeSq;
257 
258  }
259  if (sigma < DBL_MIN) lambda = DBL_MAX;
260  else lambda = 1./sigma;
261  fLambda = lambda;
262  fGamma = gamma;
263  if(verboseLevel > 1)
264  {
265  G4cout<<" lambda = "<<lambda/mm<<" mm"<<G4endl;
266  }
267  }
268  }
269  }
270  return lambda;
271 }
272 
274 //
275 // Interface for build table from physics list
276 
278 {
279  if(pd.GetPDGCharge() == 0.)
280  {
281  G4Exception("G4VXTRenergyLoss::BuildPhysicsTable", "Notification", JustWarning,
282  "XTR initialisation for neutral particle ?!" );
283  }
285 
286  if ( fAngleRadDistr )
287  {
288  if(verboseLevel > 0)
289  {
290  G4cout<<"Build angle for energy distribution according the current radiator"
291  <<G4endl;
292  }
294  }
295 }
296 
297 
299 //
300 // Build integral energy distribution of XTR photons
301 
303 {
304  G4int iTkin, iTR, iPlace;
305  G4double radiatorCof = 1.0; // for tuning of XTR yield
306  G4double energySum = 0.0;
307 
310 
311  fGammaTkinCut = 0.0;
312 
313  // setting of min/max TR energies
314 
317 
320 
322 
323  G4cout.precision(4);
324  G4Timer timer;
325  timer.Start();
326 
327  if(verboseLevel > 0)
328  {
329  G4cout<<G4endl;
330  G4cout<<"Lorentz Factor"<<"\t"<<"XTR photon number"<<G4endl;
331  G4cout<<G4endl;
332  }
333  for( iTkin = 0; iTkin < fTotBin; iTkin++ ) // Lorentz factor loop
334  {
336  fMaxEnergyTR,
337  fBinTR );
338 
339  fGamma = 1.0 + (fProtonEnergyVector->
340  GetLowEdgeEnergy(iTkin)/proton_mass_c2);
341 
342  fMaxThetaTR = 25.*2500.0/(fGamma*fGamma) ; // theta^2
343 
344  fTheMinAngle = 1.0e-3; // was 5.e-6, e-6 !!!, e-5, e-4
345 
348 
349  energySum = 0.0;
350 
351  energyVector->PutValue(fBinTR-1,energySum);
352 
353  for( iTR = fBinTR - 2; iTR >= 0; iTR-- )
354  {
355  // Legendre96 or Legendre10
356 
357  energySum += radiatorCof*fCofTR*integral.Legendre10(
359  energyVector->GetLowEdgeEnergy(iTR),
360  energyVector->GetLowEdgeEnergy(iTR+1) );
361 
362  energyVector->PutValue(iTR,energySum/fTotalDist);
363  }
364  iPlace = iTkin;
365  fEnergyDistrTable->insertAt(iPlace,energyVector);
366 
367  if(verboseLevel > 0)
368  {
369  G4cout
370  // <<iTkin<<"\t"
371  // <<"fGamma = "
372  <<fGamma<<"\t" // <<" fMaxThetaTR = "<<fMaxThetaTR
373  // <<"sumN = "
374  <<energySum // <<"; sumA = "<<angleSum
375  <<G4endl;
376  }
377  }
378  timer.Stop();
379  G4cout.precision(6);
380  if(verboseLevel > 0)
381  {
382  G4cout<<G4endl;
383  G4cout<<"total time for build X-ray TR energy loss tables = "
384  <<timer.GetUserElapsed()<<" s"<<G4endl;
385  }
386  fGamma = 0.;
387  return;
388 }
389 
391 //
392 // Bank of angle distributions for given energies (slow!)
393 
395 {
396  if( this->GetProcessName() == "TranspRegXTRadiator" ||
397  this->GetProcessName() == "TranspRegXTRmodel" ||
398  this->GetProcessName() == "RegularXTRadiator" ||
399  this->GetProcessName() == "RegularXTRmodel" )
400  {
401  BuildAngleTable();
402  return;
403  }
404  G4int i, iTkin, iTR;
405  G4double angleSum = 0.0;
406 
407 
408  fGammaTkinCut = 0.0;
409 
410  // setting of min/max TR energies
411 
414 
417 
419  fMaxEnergyTR,
420  fBinTR );
421 
423 
424  G4cout.precision(4);
425  G4Timer timer;
426  timer.Start();
427 
428  for( iTkin = 0; iTkin < fTotBin; iTkin++ ) // Lorentz factor loop
429  {
430 
431  fGamma = 1.0 + (fProtonEnergyVector->
432  GetLowEdgeEnergy(iTkin)/proton_mass_c2);
433 
434  fMaxThetaTR = 25*2500.0/(fGamma*fGamma) ; // theta^2
435 
436  fTheMinAngle = 1.0e-3; // was 5.e-6, e-6 !!!, e-5, e-4
437 
440 
442 
443  for( iTR = 0; iTR < fBinTR; iTR++ )
444  {
445  angleSum = 0.0;
446  fEnergy = energyVector->GetLowEdgeEnergy(iTR);
447  G4PhysicsLinearVector* angleVector = new G4PhysicsLinearVector(0.0,
448  fMaxThetaTR,
449  fBinTR );
450 
451  angleVector ->PutValue(fBinTR - 1, angleSum);
452 
453  for( i = fBinTR - 2; i >= 0; i-- )
454  {
455  // Legendre96 or Legendre10
456 
457  angleSum += integral.Legendre10(
459  angleVector->GetLowEdgeEnergy(i),
460  angleVector->GetLowEdgeEnergy(i+1) );
461 
462  angleVector ->PutValue(i, angleSum);
463  }
464  fAngleForEnergyTable->insertAt(iTR, angleVector);
465  }
466  fAngleBank.push_back(fAngleForEnergyTable);
467  }
468  timer.Stop();
469  G4cout.precision(6);
470  if(verboseLevel > 0)
471  {
472  G4cout<<G4endl;
473  G4cout<<"total time for build X-ray TR angle for energy loss tables = "
474  <<timer.GetUserElapsed()<<" s"<<G4endl;
475  }
476  fGamma = 0.;
477  delete energyVector;
478 }
479 
481 //
482 // Build XTR angular distribution at given energy based on the model
483 // of transparent regular radiator
484 
486 {
487  G4int iTkin, iTR;
489 
490  fGammaTkinCut = 0.0;
491 
492  // setting of min/max TR energies
493 
496 
499 
500  G4cout.precision(4);
501  G4Timer timer;
502  timer.Start();
503  if(verboseLevel > 0)
504  {
505  G4cout<<G4endl;
506  G4cout<<"Lorentz Factor"<<"\t"<<"XTR photon number"<<G4endl;
507  G4cout<<G4endl;
508  }
509  for( iTkin = 0; iTkin < fTotBin; iTkin++ ) // Lorentz factor loop
510  {
511 
512  fGamma = 1.0 + (fProtonEnergyVector->
513  GetLowEdgeEnergy(iTkin)/proton_mass_c2);
514 
515  fMaxThetaTR = 25*2500.0/(fGamma*fGamma); // theta^2
516 
517  fTheMinAngle = 1.0e-3; // was 5.e-6, e-6 !!!, e-5, e-4
518 
520  else
521  {
523  }
524 
526 
527  for( iTR = 0; iTR < fBinTR; iTR++ )
528  {
529  // energy = fMinEnergyTR*(iTR+1);
530 
531  energy = fXTREnergyVector->GetLowEdgeEnergy(iTR);
532 
533  G4PhysicsFreeVector* angleVector = GetAngleVector(energy,fBinTR);
534  // G4cout<<G4endl;
535 
536  fAngleForEnergyTable->insertAt(iTR,angleVector);
537  }
538 
539  fAngleBank.push_back(fAngleForEnergyTable);
540  }
541  timer.Stop();
542  G4cout.precision(6);
543  if(verboseLevel > 0)
544  {
545  G4cout<<G4endl;
546  G4cout<<"total time for build XTR angle for given energy tables = "
547  <<timer.GetUserElapsed()<<" s"<<G4endl;
548  }
549  fGamma = 0.;
550 
551  return;
552 }
553 
555 //
556 // Vector of angles and angle integral distributions
557 
559 {
560  G4double theta=0., result, tmp=0., cof1, cof2, cofMin, cofPHC, angleSum = 0.;
561  G4int iTheta, k, /*kMax,*/ kMin;
562 
563  G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(n);
564 
565  cofPHC = 4.*pi*hbarc;
566  tmp = (fSigma1 - fSigma2)/cofPHC/energy;
567  cof1 = fPlateThick*tmp;
568  cof2 = fGasThick*tmp;
569 
570  cofMin = energy*(fPlateThick + fGasThick)/fGamma/fGamma;
571  cofMin += (fPlateThick*fSigma1 + fGasThick*fSigma2)/energy;
572  cofMin /= cofPHC;
573 
574  kMin = G4int(cofMin);
575  if (cofMin > kMin) kMin++;
576 
577  //kMax = kMin + fBinTR -1;
578 
579  if(verboseLevel > 2)
580  {
581  G4cout<<"n-1 = "<<n-1<<"; theta = "
582  <<std::sqrt(fMaxThetaTR)*fGamma<<"; tmp = "
583  <<0.
584  <<"; angleSum = "<<angleSum<<G4endl;
585  }
586  // angleVector->PutValue(n-1,fMaxThetaTR, angleSum);
587 
588  for( iTheta = n - 1; iTheta >= 1; iTheta-- )
589  {
590 
591  k = iTheta - 1 + kMin;
592 
593  tmp = pi*fPlateThick*(k + cof2)/(fPlateThick + fGasThick);
594 
595  result = (k - cof1)*(k - cof1)*(k + cof2)*(k + cof2);
596 
597  tmp = std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result;
598 
599  if( k == kMin && kMin == G4int(cofMin) )
600  {
601  angleSum += 0.5*tmp; // 0.5*std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result;
602  }
603  else if(iTheta == n-1);
604  else
605  {
606  angleSum += tmp; // std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result;
607  }
608  theta = std::abs(k-cofMin)*cofPHC/energy/(fPlateThick + fGasThick);
609 
610  if(verboseLevel > 2)
611  {
612  G4cout<<"iTheta = "<<iTheta<<"; k = "<<k<<"; theta = "
613  <<std::sqrt(theta)*fGamma<<"; tmp = "
614  <<tmp // std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result
615  <<"; angleSum = "<<angleSum<<G4endl;
616  }
617  angleVector->PutValue( iTheta, theta, angleSum );
618  }
619  if (theta > 0.)
620  {
621  angleSum += 0.5*tmp;
622  theta = 0.;
623  }
624  if(verboseLevel > 2)
625  {
626  G4cout<<"iTheta = "<<iTheta<<"; theta = "
627  <<std::sqrt(theta)*fGamma<<"; tmp = "
628  <<tmp
629  <<"; angleSum = "<<angleSum<<G4endl;
630  }
631  angleVector->PutValue( iTheta, theta, angleSum );
632 
633  return angleVector;
634 }
635 
637 //
638 // Build XTR angular distribution based on the model of transparent regular radiator
639 
641 {
642  G4int iTkin, iTR, iPlace;
643  G4double radiatorCof = 1.0; // for tuning of XTR yield
644  G4double angleSum;
646 
647  fGammaTkinCut = 0.0;
648 
649  // setting of min/max TR energies
650 
653 
656 
657  G4cout.precision(4);
658  G4Timer timer;
659  timer.Start();
660  if(verboseLevel > 0) {
661  G4cout<<G4endl;
662  G4cout<<"Lorentz Factor"<<"\t"<<"XTR photon number"<<G4endl;
663  G4cout<<G4endl;
664  }
665  for( iTkin = 0; iTkin < fTotBin; iTkin++ ) // Lorentz factor loop
666  {
667 
668  fGamma = 1.0 + (fProtonEnergyVector->
669  GetLowEdgeEnergy(iTkin)/proton_mass_c2);
670 
671  fMaxThetaTR = 25.0/(fGamma*fGamma); // theta^2
672 
673  fTheMinAngle = 1.0e-3; // was 5.e-6, e-6 !!!, e-5, e-4
674 
676  else
677  {
679  }
680  G4PhysicsLinearVector* angleVector = new G4PhysicsLinearVector(0.0,
681  fMaxThetaTR,
682  fBinTR );
683 
684  angleSum = 0.0;
685 
687 
688 
689  angleVector->PutValue(fBinTR-1,angleSum);
690 
691  for( iTR = fBinTR - 2; iTR >= 0; iTR-- )
692  {
693 
694  angleSum += radiatorCof*fCofTR*integral.Legendre96(
696  angleVector->GetLowEdgeEnergy(iTR),
697  angleVector->GetLowEdgeEnergy(iTR+1) );
698 
699  angleVector ->PutValue(iTR,angleSum);
700  }
701  if(verboseLevel > 1) {
702  G4cout
703  // <<iTkin<<"\t"
704  // <<"fGamma = "
705  <<fGamma<<"\t" // <<" fMaxThetaTR = "<<fMaxThetaTR
706  // <<"sumN = "<<energySum // <<"; sumA = "
707  <<angleSum
708  <<G4endl;
709  }
710  iPlace = iTkin;
711  fAngleDistrTable->insertAt(iPlace,angleVector);
712  }
713  timer.Stop();
714  G4cout.precision(6);
715  if(verboseLevel > 0) {
716  G4cout<<G4endl;
717  G4cout<<"total time for build X-ray TR angle tables = "
718  <<timer.GetUserElapsed()<<" s"<<G4endl;
719  }
720  fGamma = 0.;
721 
722  return;
723 }
724 
725 
727 //
728 // The main function which is responsible for the treatment of a particle passage
729 // trough G4Envelope with discrete generation of G4Gamma
730 
732  const G4Step& aStep )
733 {
734  G4int iTkin /*, iPlace*/;
735  G4double energyTR, theta,theta2, phi, dirX, dirY, dirZ;
736 
737 
738  fParticleChange.Initialize(aTrack);
739 
740  if(verboseLevel > 1)
741  {
742  G4cout<<"Start of G4VXTRenergyLoss::PostStepDoIt "<<G4endl;
743  G4cout<<"name of current material = "
744  <<aTrack.GetVolume()->GetLogicalVolume()->GetMaterial()->GetName()<<G4endl;
745  }
746  if( aTrack.GetVolume()->GetLogicalVolume() != fEnvelope )
747  {
748  if(verboseLevel > 0)
749  {
750  G4cout<<"Go out from G4VXTRenergyLoss::PostStepDoIt: wrong volume "<<G4endl;
751  }
752  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
753  }
754  else
755  {
756  G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
757  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
758 
759  // Now we are ready to Generate one TR photon
760 
761  G4double kinEnergy = aParticle->GetKineticEnergy();
762  G4double mass = aParticle->GetDefinition()->GetPDGMass();
763  G4double gamma = 1.0 + kinEnergy/mass;
764 
765  if(verboseLevel > 1 )
766  {
767  G4cout<<"gamma = "<<gamma<<G4endl;
768  }
769  G4double massRatio = proton_mass_c2/mass;
770  G4double TkinScaled = kinEnergy*massRatio;
771  G4ThreeVector position = pPostStepPoint->GetPosition();
772  G4ParticleMomentum direction = aParticle->GetMomentumDirection();
773  G4double startTime = pPostStepPoint->GetGlobalTime();
774 
775  for( iTkin = 0; iTkin < fTotBin; iTkin++ )
776  {
777  if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
778  }
779  //iPlace = iTkin - 1;
780 
781  if(iTkin == 0) // Tkin is too small, neglect of TR photon generation
782  {
783  if( verboseLevel > 0)
784  {
785  G4cout<<"Go out from G4VXTRenergyLoss::PostStepDoIt:iTkin = "<<iTkin<<G4endl;
786  }
787  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
788  }
789  else // general case: Tkin between two vectors of the material
790  {
792 
793  energyTR = GetXTRrandomEnergy(TkinScaled,iTkin);
794 
795  if( verboseLevel > 1)
796  {
797  G4cout<<"energyTR = "<<energyTR/keV<<" keV"<<G4endl;
798  }
799  if ( fAngleRadDistr )
800  {
801  // theta = std::fabs(G4RandGauss::shoot(0.0,pi/gamma));
802 
803  theta2 = GetRandomAngle(energyTR,iTkin);
804 
805  if( theta2 > 0.) theta = std::sqrt(theta2);
806  else theta = 0.; // theta2;
807  }
808  else theta = std::fabs(G4RandGauss::shoot(0.0,pi/gamma));
809 
810  // theta = 0.; // check no spread
811 
812  if( theta >= 0.1 ) theta = 0.1;
813 
814  // G4cout<<" : theta = "<<theta<<endl;
815 
816  phi = twopi*G4UniformRand();
817 
818  dirX = std::sin(theta)*std::cos(phi);
819  dirY = std::sin(theta)*std::sin(phi);
820  dirZ = std::cos(theta);
821 
822  G4ThreeVector directionTR(dirX,dirY,dirZ);
823  directionTR.rotateUz(direction);
824  directionTR.unit();
825 
827  directionTR, energyTR);
828 
829  // A XTR photon is set on the particle track inside the radiator
830  // and is moved to the G4Envelope surface for standard X-ray TR models
831  // only. The case of fExitFlux=true
832 
833  if( fExitFlux )
834  {
835  const G4RotationMatrix* rotM = pPostStepPoint->GetTouchable()->GetRotation();
836  G4ThreeVector transl = pPostStepPoint->GetTouchable()->GetTranslation();
837  G4AffineTransform transform = G4AffineTransform(rotM,transl);
838  transform.Invert();
839  G4ThreeVector localP = transform.TransformPoint(position);
840  G4ThreeVector localV = transform.TransformAxis(directionTR);
841 
842  G4double distance = fEnvelope->GetSolid()->DistanceToOut(localP, localV);
843  if(verboseLevel > 1)
844  {
845  G4cout<<"distance to exit = "<<distance/mm<<" mm"<<G4endl;
846  }
847  position += distance*directionTR;
848  startTime += distance/c_light;
849  }
850  G4Track* aSecondaryTrack = new G4Track( aPhotonTR,
851  startTime, position );
852  aSecondaryTrack->SetTouchableHandle(
854  aSecondaryTrack->SetParentID( aTrack.GetTrackID() );
855 
856  fParticleChange.AddSecondary(aSecondaryTrack);
857  fParticleChange.ProposeEnergy(kinEnergy);
858  }
859  }
860  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
861 }
862 
864 //
865 // This function returns the spectral and angle density of TR quanta
866 // in X-ray energy region generated forward when a relativistic
867 // charged particle crosses interface between two materials.
868 // The high energy small theta approximation is applied.
869 // (matter1 -> matter2, or 2->1)
870 // varAngle =2* (1 - std::cos(theta)) or approximately = theta*theta
871 //
872 
874  G4double gamma,
875  G4double varAngle )
876 {
877  G4complex Z1 = GetPlateComplexFZ(energy,gamma,varAngle);
878  G4complex Z2 = GetGasComplexFZ(energy,gamma,varAngle);
879 
880  G4complex zOut = (Z1 - Z2)*(Z1 - Z2)
881  * (varAngle*energy/hbarc/hbarc);
882  return zOut;
883 
884 }
885 
886 
888 //
889 // For photon energy distribution tables. Integrate first over angle
890 //
891 
893 {
895  if(result < 0.0) result = 0.0;
896  return result;
897 }
898 
900 //
901 // For second integration over energy
902 
904 {
905  G4int i, iMax = 8;
906  G4double angleSum = 0.0;
907 
908  G4double lim[8] = { 0.0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0 };
909 
910  for( i = 0; i < iMax; i++ ) lim[i] *= fMaxThetaTR;
911 
913 
914  fEnergy = energy;
915  /*
916  if( fAngleRadDistr && ( fEnergy == fEnergyForAngle ) )
917  {
918  fAngleVector ->PutValue(fBinTR - 1, angleSum);
919 
920  for( i = fBinTR - 2; i >= 0; i-- )
921  {
922 
923  angleSum += integral.Legendre10(
924  this,&G4VXTRenergyLoss::SpectralAngleXTRdEdx,
925  fAngleVector->GetLowEdgeEnergy(i),
926  fAngleVector->GetLowEdgeEnergy(i+1) );
927 
928  fAngleVector ->PutValue(i, angleSum);
929  }
930  }
931  else
932  */
933  {
934  for( i = 0; i < iMax-1; i++ )
935  {
936  angleSum += integral.Legendre96(this,&G4VXTRenergyLoss::SpectralAngleXTRdEdx,
937  lim[i],lim[i+1]);
938  // result += integral.Legendre10(this,&G4VXTRenergyLoss::SpectralAngleXTRdEdx,
939  // lim[i],lim[i+1]);
940  }
941  }
942  return angleSum;
943 }
944 
946 //
947 // for photon angle distribution tables
948 //
949 
951 {
953  if(result < 0) result = 0.0;
954  return result;
955 }
956 
958 //
959 // The XTR angular distribution based on transparent regular radiator
960 
962 {
963  // G4cout<<"angle2 = "<<varAngle<<"; fGamma = "<<fGamma<<G4endl;
964 
966  G4double sum = 0., tmp1, tmp2, tmp=0., cof1, cof2, cofMin, cofPHC, energy1, energy2;
967  G4int k, kMax, kMin, i;
968 
969  cofPHC = twopi*hbarc;
970 
971  cof1 = (fPlateThick + fGasThick)*(1./fGamma/fGamma + varAngle);
973 
974  // G4cout<<"cof1 = "<<cof1<<"; cof2 = "<<cof2<<"; cofPHC = "<<cofPHC<<G4endl;
975 
976  cofMin = std::sqrt(cof1*cof2);
977  cofMin /= cofPHC;
978 
979  kMin = G4int(cofMin);
980  if (cofMin > kMin) kMin++;
981 
982  kMax = kMin + 9; // 9; // kMin + G4int(tmp);
983 
984  // G4cout<<"cofMin = "<<cofMin<<"; kMin = "<<kMin<<"; kMax = "<<kMax<<G4endl;
985 
986  for( k = kMin; k <= kMax; k++ )
987  {
988  tmp1 = cofPHC*k;
989  tmp2 = std::sqrt(tmp1*tmp1-cof1*cof2);
990  energy1 = (tmp1+tmp2)/cof1;
991  energy2 = (tmp1-tmp2)/cof1;
992 
993  for(i = 0; i < 2; i++)
994  {
995  if( i == 0 )
996  {
997  if (energy1 > fTheMaxEnergyTR || energy1 < fTheMinEnergyTR) continue;
998 
999  tmp1 = ( energy1*energy1*(1./fGamma/fGamma + varAngle) + fSigma1 )
1000  * fPlateThick/(4*hbarc*energy1);
1001  tmp2 = std::sin(tmp1);
1002  tmp = energy1*tmp2*tmp2;
1003  tmp2 = fPlateThick/(4.*tmp1);
1004  tmp1 = hbarc*energy1/( energy1*energy1*(1./fGamma/fGamma + varAngle) + fSigma2 );
1005  tmp *= (tmp1-tmp2)*(tmp1-tmp2);
1006  tmp1 = cof1/(4.*hbarc) - cof2/(4.*hbarc*energy1*energy1);
1007  tmp2 = std::abs(tmp1);
1008 
1009  if(tmp2 > 0.) tmp /= tmp2;
1010  else continue;
1011  }
1012  else
1013  {
1014  if (energy2 > fTheMaxEnergyTR || energy2 < fTheMinEnergyTR) continue;
1015 
1016  tmp1 = ( energy2*energy2*(1./fGamma/fGamma + varAngle) + fSigma1 )
1017  * fPlateThick/(4.*hbarc*energy2);
1018  tmp2 = std::sin(tmp1);
1019  tmp = energy2*tmp2*tmp2;
1020  tmp2 = fPlateThick/(4.*tmp1);
1021  tmp1 = hbarc*energy2/( energy2*energy2*(1./fGamma/fGamma + varAngle) + fSigma2 );
1022  tmp *= (tmp1-tmp2)*(tmp1-tmp2);
1023  tmp1 = cof1/(4.*hbarc) - cof2/(4.*hbarc*energy2*energy2);
1024  tmp2 = std::abs(tmp1);
1025 
1026  if(tmp2 > 0.) tmp /= tmp2;
1027  else continue;
1028  }
1029  sum += tmp;
1030  }
1031  // G4cout<<"k = "<<k<<"; energy1 = "<<energy1/keV<<" keV; energy2 = "<<energy2/keV
1032  // <<" keV; tmp = "<<tmp<<"; sum = "<<sum<<G4endl;
1033  }
1034  result = 4.*pi*fPlateNumber*sum*varAngle;
1035  result /= hbarc*hbarc;
1036 
1037  // old code based on general numeric integration
1038  // fVarAngle = varAngle;
1039  // G4Integrator<G4VXTRenergyLoss,G4double(G4VXTRenergyLoss::*)(G4double)> integral;
1040  // result = integral.Legendre10(this,&G4VXTRenergyLoss::AngleSpectralXTRdEdx,
1041  // fMinEnergyTR,fMaxEnergyTR);
1042  return result;
1043 }
1044 
1045 
1049 //
1050 // Calculates formation zone for plates. Omega is energy !!!
1051 
1053  G4double gamma ,
1054  G4double varAngle )
1055 {
1056  G4double cof, lambda;
1057  lambda = 1.0/gamma/gamma + varAngle + fSigma1/omega/omega;
1058  cof = 2.0*hbarc/omega/lambda;
1059  return cof;
1060 }
1061 
1063 //
1064 // Calculates complex formation zone for plates. Omega is energy !!!
1065 
1067  G4double gamma ,
1068  G4double varAngle )
1069 {
1070  G4double cof, length,delta, real_v, image_v;
1071 
1072  length = 0.5*GetPlateFormationZone(omega,gamma,varAngle);
1073  delta = length*GetPlateLinearPhotoAbs(omega);
1074  cof = 1.0/(1.0 + delta*delta);
1075 
1076  real_v = length*cof;
1077  image_v = real_v*delta;
1078 
1079  G4complex zone(real_v,image_v);
1080  return zone;
1081 }
1082 
1084 //
1085 // Computes matrix of Sandia photo absorption cross section coefficients for
1086 // plate material
1087 
1089 {
1090  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1091  const G4Material* mat = (*theMaterialTable)[fMatIndex1];
1093 
1094  return;
1095 }
1096 
1097 
1098 
1100 //
1101 // Returns the value of linear photo absorption coefficient (in reciprocal
1102 // length) for plate for given energy of X-ray photon omega
1103 
1105 {
1106  // G4int i;
1107  G4double omega2, omega3, omega4;
1108 
1109  omega2 = omega*omega;
1110  omega3 = omega2*omega;
1111  omega4 = omega2*omega2;
1112 
1113  const G4double* SandiaCof = fPlatePhotoAbsCof->GetSandiaCofForMaterial(omega);
1114  G4double cross = SandiaCof[0]/omega + SandiaCof[1]/omega2 +
1115  SandiaCof[2]/omega3 + SandiaCof[3]/omega4;
1116  return cross;
1117 }
1118 
1119 
1121 //
1122 // Calculates formation zone for gas. Omega is energy !!!
1123 
1125  G4double gamma ,
1126  G4double varAngle )
1127 {
1128  G4double cof, lambda;
1129  lambda = 1.0/gamma/gamma + varAngle + fSigma2/omega/omega;
1130  cof = 2.0*hbarc/omega/lambda;
1131  return cof;
1132 }
1133 
1134 
1136 //
1137 // Calculates complex formation zone for gas gaps. Omega is energy !!!
1138 
1140  G4double gamma ,
1141  G4double varAngle )
1142 {
1143  G4double cof, length,delta, real_v, image_v;
1144 
1145  length = 0.5*GetGasFormationZone(omega,gamma,varAngle);
1146  delta = length*GetGasLinearPhotoAbs(omega);
1147  cof = 1.0/(1.0 + delta*delta);
1148 
1149  real_v = length*cof;
1150  image_v = real_v*delta;
1151 
1152  G4complex zone(real_v,image_v);
1153  return zone;
1154 }
1155 
1156 
1157 
1159 //
1160 // Computes matrix of Sandia photo absorption cross section coefficients for
1161 // gas material
1162 
1164 {
1165  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1166  const G4Material* mat = (*theMaterialTable)[fMatIndex2];
1168  return;
1169 }
1170 
1172 //
1173 // Returns the value of linear photo absorption coefficient (in reciprocal
1174 // length) for gas
1175 
1177 {
1178  G4double omega2, omega3, omega4;
1179 
1180  omega2 = omega*omega;
1181  omega3 = omega2*omega;
1182  omega4 = omega2*omega2;
1183 
1184  const G4double* SandiaCof = fGasPhotoAbsCof->GetSandiaCofForMaterial(omega);
1185  G4double cross = SandiaCof[0]/omega + SandiaCof[1]/omega2 +
1186  SandiaCof[2]/omega3 + SandiaCof[3]/omega4;
1187  return cross;
1188 
1189 }
1190 
1192 //
1193 // Calculates the product of linear cof by formation zone for plate.
1194 // Omega is energy !!!
1195 
1197  G4double gamma ,
1198  G4double varAngle )
1199 {
1200  return GetPlateFormationZone(omega,gamma,varAngle)
1201  * GetPlateLinearPhotoAbs(omega);
1202 }
1204 //
1205 // Calculates the product of linear cof by formation zone for plate.
1206 // G4cout and output in file in some energy range.
1207 
1209 {
1210  std::ofstream outPlate("plateZmu.dat", std::ios::out );
1211  outPlate.setf( std::ios::scientific, std::ios::floatfield );
1212 
1213  G4int i;
1214  G4double omega, varAngle, gamma;
1215  gamma = 10000.;
1216  varAngle = 1/gamma/gamma;
1217  if(verboseLevel > 0)
1218  G4cout<<"energy, keV"<<"\t"<<"Zmu for plate"<<G4endl;
1219  for(i=0;i<100;i++)
1220  {
1221  omega = (1.0 + i)*keV;
1222  if(verboseLevel > 1)
1223  G4cout<<omega/keV<<"\t"<<GetPlateZmuProduct(omega,gamma,varAngle)<<"\t";
1224  if(verboseLevel > 0)
1225  outPlate<<omega/keV<<"\t\t"<<GetPlateZmuProduct(omega,gamma,varAngle)<<G4endl;
1226  }
1227  return;
1228 }
1229 
1231 //
1232 // Calculates the product of linear cof by formation zone for gas.
1233 // Omega is energy !!!
1234 
1236  G4double gamma ,
1237  G4double varAngle )
1238 {
1239  return GetGasFormationZone(omega,gamma,varAngle)*GetGasLinearPhotoAbs(omega);
1240 }
1242 //
1243 // Calculates the product of linear cof byformation zone for gas.
1244 // G4cout and output in file in some energy range.
1245 
1247 {
1248  std::ofstream outGas("gasZmu.dat", std::ios::out );
1249  outGas.setf( std::ios::scientific, std::ios::floatfield );
1250  G4int i;
1251  G4double omega, varAngle, gamma;
1252  gamma = 10000.;
1253  varAngle = 1/gamma/gamma;
1254  if(verboseLevel > 0)
1255  G4cout<<"energy, keV"<<"\t"<<"Zmu for gas"<<G4endl;
1256  for(i=0;i<100;i++)
1257  {
1258  omega = (1.0 + i)*keV;
1259  if(verboseLevel > 1)
1260  G4cout<<omega/keV<<"\t"<<GetGasZmuProduct(omega,gamma,varAngle)<<"\t";
1261  if(verboseLevel > 0)
1262  outGas<<omega/keV<<"\t\t"<<GetGasZmuProduct(omega,gamma,varAngle)<<G4endl;
1263  }
1264  return;
1265 }
1266 
1268 //
1269 // Computes Compton cross section for plate material in 1/mm
1270 
1272 {
1273  G4int i, numberOfElements;
1274  G4double xSection = 0., nowZ, sumZ = 0.;
1275 
1276  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1277  numberOfElements = (*theMaterialTable)[fMatIndex1]->GetNumberOfElements();
1278 
1279  for( i = 0; i < numberOfElements; i++ )
1280  {
1281  nowZ = (*theMaterialTable)[fMatIndex1]->GetElement(i)->GetZ();
1282  sumZ += nowZ;
1283  xSection += GetComptonPerAtom(omega,nowZ); // *nowZ;
1284  }
1285  xSection /= sumZ;
1286  xSection *= (*theMaterialTable)[fMatIndex1]->GetElectronDensity();
1287  return xSection;
1288 }
1289 
1290 
1292 //
1293 // Computes Compton cross section for gas material in 1/mm
1294 
1296 {
1297  G4int i, numberOfElements;
1298  G4double xSection = 0., nowZ, sumZ = 0.;
1299 
1300  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1301  numberOfElements = (*theMaterialTable)[fMatIndex2]->GetNumberOfElements();
1302 
1303  for( i = 0; i < numberOfElements; i++ )
1304  {
1305  nowZ = (*theMaterialTable)[fMatIndex2]->GetElement(i)->GetZ();
1306  sumZ += nowZ;
1307  xSection += GetComptonPerAtom(omega,nowZ); // *nowZ;
1308  }
1309  xSection /= sumZ;
1310  xSection *= (*theMaterialTable)[fMatIndex2]->GetElectronDensity();
1311  return xSection;
1312 }
1313 
1315 //
1316 // Computes Compton cross section per atom with Z electrons for gamma with
1317 // the energy GammaEnergy
1318 
1320 {
1321  G4double CrossSection = 0.0;
1322  if ( Z < 0.9999 ) return CrossSection;
1323  if ( GammaEnergy < 0.1*keV ) return CrossSection;
1324  if ( GammaEnergy > (100.*GeV/Z) ) return CrossSection;
1325 
1326  static const G4double a = 20.0 , b = 230.0 , c = 440.0;
1327 
1328  static const G4double
1329  d1= 2.7965e-1*barn, d2=-1.8300e-1*barn, d3= 6.7527 *barn, d4=-1.9798e+1*barn,
1330  e1= 1.9756e-5*barn, e2=-1.0205e-2*barn, e3=-7.3913e-2*barn, e4= 2.7079e-2*barn,
1331  f1=-3.9178e-7*barn, f2= 6.8241e-5*barn, f3= 6.0480e-5*barn, f4= 3.0274e-4*barn;
1332 
1333  G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z),
1334  p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z);
1335 
1336  G4double T0 = 15.0*keV;
1337  if (Z < 1.5) T0 = 40.0*keV;
1338 
1339  G4double X = std::max(GammaEnergy, T0) / electron_mass_c2;
1340  CrossSection = p1Z*std::log(1.+2.*X)/X
1341  + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
1342 
1343  // modification for low energy. (special case for Hydrogen)
1344 
1345  if (GammaEnergy < T0)
1346  {
1347  G4double dT0 = 1.*keV;
1348  X = (T0+dT0) / electron_mass_c2;
1349  G4double sigma = p1Z*std::log(1.+2.*X)/X
1350  + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
1351  G4double c1 = -T0*(sigma-CrossSection)/(CrossSection*dT0);
1352  G4double c2 = 0.150;
1353  if (Z > 1.5) c2 = 0.375-0.0556*std::log(Z);
1354  G4double y = std::log(GammaEnergy/T0);
1355  CrossSection *= std::exp(-y*(c1+c2*y));
1356  }
1357  // G4cout << "e= " << GammaEnergy << " Z= " << Z << " cross= " << CrossSection << G4endl;
1358  return CrossSection;
1359 }
1360 
1361 
1363 //
1364 // This function returns the spectral and angle density of TR quanta
1365 // in X-ray energy region generated forward when a relativistic
1366 // charged particle crosses interface between two materials.
1367 // The high energy small theta approximation is applied.
1368 // (matter1 -> matter2, or 2->1)
1369 // varAngle =2* (1 - std::cos(theta)) or approximately = theta*theta
1370 //
1371 
1372 G4double
1374  G4double varAngle ) const
1375 {
1376  G4double formationLength1, formationLength2;
1377  formationLength1 = 1.0/
1378  (1.0/(gamma*gamma)
1379  + fSigma1/(energy*energy)
1380  + varAngle);
1381  formationLength2 = 1.0/
1382  (1.0/(gamma*gamma)
1383  + fSigma2/(energy*energy)
1384  + varAngle);
1385  return (varAngle/energy)*(formationLength1 - formationLength2)
1386  *(formationLength1 - formationLength2);
1387 
1388 }
1389 
1391  G4double varAngle )
1392 {
1393  // return stack factor corresponding to one interface
1394 
1395  return std::real( OneInterfaceXTRdEdx(energy,gamma,varAngle) );
1396 }
1397 
1399 //
1400 // For photon energy distribution tables. Integrate first over angle
1401 //
1402 
1404 {
1405  return OneBoundaryXTRNdensity(fEnergy,fGamma,varAngle)*
1406  GetStackFactor(fEnergy,fGamma,varAngle);
1407 }
1408 
1410 //
1411 // For second integration over energy
1412 
1414 {
1415  fEnergy = energy;
1418  0.0,0.2*fMaxThetaTR) +
1420  0.2*fMaxThetaTR,fMaxThetaTR);
1421 }
1422 
1424 //
1425 // for photon angle distribution tables
1426 //
1427 
1429 {
1430  return OneBoundaryXTRNdensity(energy,fGamma,fVarAngle)*
1432 }
1433 
1435 //
1436 //
1437 
1439 {
1440  fVarAngle = varAngle;
1444 }
1445 
1447 //
1448 // Check number of photons for a range of Lorentz factors from both energy
1449 // and angular tables
1450 
1452 {
1453  G4int iTkin;
1454  G4double gamma, numberE;
1455 
1456  std::ofstream outEn("numberE.dat", std::ios::out );
1457  outEn.setf( std::ios::scientific, std::ios::floatfield );
1458 
1459  std::ofstream outAng("numberAng.dat", std::ios::out );
1460  outAng.setf( std::ios::scientific, std::ios::floatfield );
1461 
1462  for(iTkin=0;iTkin<fTotBin;iTkin++) // Lorentz factor loop
1463  {
1464  gamma = 1.0 + (fProtonEnergyVector->
1465  GetLowEdgeEnergy(iTkin)/proton_mass_c2);
1466  numberE = (*(*fEnergyDistrTable)(iTkin))(0);
1467  // numberA = (*(*fAngleDistrTable)(iTkin))(0);
1468  if(verboseLevel > 1)
1469  G4cout<<gamma<<"\t\t"<<numberE<<"\t" // <<numberA
1470  <<G4endl;
1471  if(verboseLevel > 0)
1472  outEn<<gamma<<"\t\t"<<numberE<<G4endl;
1473  }
1474  return;
1475 }
1476 
1478 //
1479 // Returns randon energy of a X-ray TR photon for given scaled kinetic energy
1480 // of a charged particle
1481 
1483 {
1484  G4int iTransfer, iPlace;
1485  G4double transfer = 0.0, position, E1, E2, W1, W2, W;
1486 
1487  iPlace = iTkin - 1;
1488 
1489  // G4cout<<"iPlace = "<<iPlace<<endl;
1490 
1491  if(iTkin == fTotBin) // relativistic plato, try from left
1492  {
1493  position = (*(*fEnergyDistrTable)(iPlace))(0)*G4UniformRand();
1494 
1495  for(iTransfer=0;;iTransfer++)
1496  {
1497  if(position >= (*(*fEnergyDistrTable)(iPlace))(iTransfer)) break;
1498  }
1499  transfer = GetXTRenergy(iPlace,position,iTransfer);
1500  }
1501  else
1502  {
1503  E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1);
1505  W = 1.0/(E2 - E1);
1506  W1 = (E2 - scaledTkin)*W;
1507  W2 = (scaledTkin - E1)*W;
1508 
1509  position =( (*(*fEnergyDistrTable)(iPlace))(0)*W1 +
1510  (*(*fEnergyDistrTable)(iPlace+1))(0)*W2 )*G4UniformRand();
1511 
1512  // G4cout<<position<<"\t";
1513 
1514  for(iTransfer=0;;iTransfer++)
1515  {
1516  if( position >=
1517  ( (*(*fEnergyDistrTable)(iPlace))(iTransfer)*W1 +
1518  (*(*fEnergyDistrTable)(iPlace+1))(iTransfer)*W2) ) break;
1519  }
1520  transfer = GetXTRenergy(iPlace,position,iTransfer);
1521 
1522  }
1523  // G4cout<<"XTR transfer = "<<transfer/keV<<" keV"<<endl;
1524  if(transfer < 0.0 ) transfer = 0.0;
1525  return transfer;
1526 }
1527 
1529 //
1530 // Returns approximate position of X-ray photon energy during random sampling
1531 // over integral energy distribution
1532 
1534  G4double /* position */,
1535  G4int iTransfer )
1536 {
1537  G4double x1, x2, y1, y2, result;
1538 
1539  if(iTransfer == 0)
1540  {
1541  result = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1542  }
1543  else
1544  {
1545  y1 = (*(*fEnergyDistrTable)(iPlace))(iTransfer-1);
1546  y2 = (*(*fEnergyDistrTable)(iPlace))(iTransfer);
1547 
1548  x1 = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
1549  x2 = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1550 
1551  if ( x1 == x2 ) result = x2;
1552  else
1553  {
1554  if ( y1 == y2 ) result = x1 + (x2 - x1)*G4UniformRand();
1555  else
1556  {
1557  // result = x1 + (position - y1)*(x2 - x1)/(y2 - y1);
1558  result = x1 + (x2 - x1)*G4UniformRand();
1559  }
1560  }
1561  }
1562  return result;
1563 }
1564 
1566 //
1567 // Get XTR photon angle at given energy and Tkin
1568 
1570 {
1571  G4int iTR, iAngle;
1573 
1574  if (iTkin == fTotBin) iTkin--;
1575 
1577 
1578  for( iTR = 0; iTR < fBinTR; iTR++ )
1579  {
1580  if( energyXTR < fXTREnergyVector->GetLowEdgeEnergy(iTR) ) break;
1581  }
1582  if (iTR == fBinTR) iTR--;
1583 
1584  position = (*(*fAngleForEnergyTable)(iTR))(0)*G4UniformRand();
1585 
1586  for( iAngle = 0;; iAngle++)
1587  {
1588  if( position >= (*(*fAngleForEnergyTable)(iTR))(iAngle) ) break;
1589  }
1590  angle = GetAngleXTR(iTR,position,iAngle);
1591  return angle;
1592 }
1593 
1595 //
1596 // Returns approximate position of X-ray photon angle at given energy during random sampling
1597 // over integral energy distribution
1598 
1600  G4double position,
1601  G4int iTransfer )
1602 {
1603  G4double x1, x2, y1, y2, result;
1604 
1605  if( iTransfer == 0 )
1606  {
1607  result = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1608  }
1609  else
1610  {
1611  y1 = (*(*fAngleForEnergyTable)(iPlace))(iTransfer-1);
1612  y2 = (*(*fAngleForEnergyTable)(iPlace))(iTransfer);
1613 
1614  x1 = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
1615  x2 = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1616 
1617  if ( x1 == x2 ) result = x2;
1618  else
1619  {
1620  if ( y1 == y2 ) result = x1 + (x2 - x1)*G4UniformRand();
1621  else
1622  {
1623  result = x1 + (position - y1)*(x2 - x1)/(y2 - y1);
1624  }
1625  }
1626  }
1627  return result;
1628 }
1629 
1630 
1631 //
1632 //
1634 
G4ParticleChange fParticleChange
G4double GetGasCompton(G4double)
Float_t f4
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double OneBoundaryXTRNdensity(G4double energy, G4double gamma, G4double varAngle) const
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const
const G4VTouchable * GetTouchable() const
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
virtual G4double SpectralXTRdEdx(G4double energy)
G4LogicalVolume * GetLogicalVolume() const
G4double GetLowEdgeEnergy(size_t binNumber) const
Float_t y1[n_points_granero]
Definition: compare.C:5
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:593
G4complex GetGasComplexFZ(G4double, G4double, G4double)
void Start()
static constexpr double keV
Definition: G4SIunits.hh:216
G4ParticleDefinition * fPtrGamma
static constexpr double mm
Definition: G4SIunits.hh:115
Float_t x1[n_points_granero]
Definition: compare.C:5
size_t GetIndex() const
Definition: G4Material.hh:262
#define G4endl
Definition: G4ios.hh:61
Float_t y
Definition: compare.C:6
const G4ThreeVector & GetMomentumDirection() const
static constexpr double hbarc
G4double GetXTRenergy(G4int iPlace, G4double position, G4int iTransfer)
G4double Legendre96(T &typeT, F f, G4double a, G4double b)
G4int GetTrackID() const
G4double GetPDGCharge() const
virtual const G4ThreeVector & GetTranslation(G4int depth=0) const =0
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
virtual const G4RotationMatrix * GetRotation(G4int depth=0) const =0
Float_t tmp
Float_t f2
void insertAt(size_t, G4PhysicsVector *)
G4double condition(const G4ErrorSymMatrix &m)
G4double GetPDGMass() const
const G4TouchableHandle & GetTouchableHandle() const
void AddSecondary(G4Track *aSecondary)
G4Material * GetMaterial() const
G4PhysicsLogVector * fXTREnergyVector
Float_t y2[n_points_geant4]
Definition: compare.C:26
G4double AngleXTRdEdx(G4double varAngle)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
static constexpr double proton_mass_c2
const G4String & GetName() const
Definition: G4Material.hh:179
Float_t Z
static constexpr double TeV
Definition: G4SIunits.hh:218
G4double SpectralAngleXTRdEdx(G4double varAngle)
Float_t f3
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4PhysicsTable * fEnergyDistrTable
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4ParticleDefinition * GetDefinition() const
#define position
Definition: xmlparse.cc:622
static const G4double d2
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:283
TCanvas * c2
Definition: plot_hist.C:75
virtual void Initialize(const G4Track &)
G4double GetGasFormationZone(G4double, G4double, G4double)
static const G4double T0[78]
G4ProcessType
static constexpr double electron_mass_c2
const G4ThreeVector & GetPosition() const
static const G4double d1
virtual G4bool IsApplicable(const G4ParticleDefinition &) override
#define G4UniformRand()
Definition: Randomize.hh:53
static constexpr double twopi
Definition: G4SIunits.hh:76
double energy
Definition: plottest35.C:25
G4SandiaTable * fGasPhotoAbsCof
G4AffineTransform & Invert()
void SetTouchableHandle(const G4TouchableHandle &apValue)
void Stop()
G4double GetAngleXTR(G4int iTR, G4double position, G4int iAngle)
const G4String & GetProcessName() const
Definition: G4VProcess.hh:411
G4double GetUserElapsed() const
Definition: G4Timer.cc:145
G4VXTRenergyLoss(G4LogicalVolume *anEnvelope, G4Material *, G4Material *, G4double, G4double, G4int, const G4String &processName="XTRenergyLoss", G4ProcessType type=fElectromagnetic)
static constexpr double eV
Definition: G4SIunits.hh:215
G4double GetGlobalTime() const
Float_t f1
virtual G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep) override
Definition: G4Step.hh:76
G4StepPoint * GetPostStepPoint() const
void ProposeEnergy(G4double finalEnergy)
Float_t mat
G4double GetPlateLinearPhotoAbs(G4double)
G4double GetPlateFormationZone(G4double, G4double, G4double)
G4double GetSandiaCofForMaterial(G4int, G4int) const
G4PhysicsLogVector * fProtonEnergyVector
ThreeVector shoot(const G4int Ap, const G4int Af)
std::vector< G4PhysicsTable * > fAngleBank
G4double GetXTRrandomEnergy(G4double scaledTkin, G4int iTkin)
Hep3Vector unit() const
G4double G4ParticleHPJENDLHEData::G4double result
G4complex GetPlateComplexFZ(G4double, G4double, G4double)
std::vector< G4Material * > G4MaterialTable
Double_t Z2
Double_t Z1
G4double XTRNAngleSpectralDensity(G4double energy)
G4double GetComptonPerAtom(G4double, G4double)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
G4double GetRandomAngle(G4double energyXTR, G4int iTkin)
int G4int
Definition: G4Types.hh:78
G4PhysicsTable * fAngleDistrTable
static constexpr double c_light
G4int verboseLevel
Definition: G4VProcess.hh:371
void SetParentID(const G4int aValue)
static constexpr double barn
Definition: G4SIunits.hh:105
virtual G4double GetStackFactor(G4double energy, G4double gamma, G4double varAngle)
G4PhysicsFreeVector * GetAngleVector(G4double energy, G4int n)
G4double GetPlateCompton(G4double)
#define DBL_MIN
Definition: templates.hh:75
virtual G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition) override
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:435
G4ForceCondition
G4double GetKineticEnergy() const
static constexpr double cm
Definition: G4SIunits.hh:119
G4GLOB_DLL std::ostream G4cout
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
G4VSolid * GetSolid() const
void clearAndDestroy()
G4SandiaTable * GetSandiaTable() const
Definition: G4Material.hh:230
Char_t n[5]
virtual void BuildPhysicsTable(const G4ParticleDefinition &) override
void PutValue(size_t index, G4double energy, G4double dataValue)
static constexpr double pi
Definition: G4SIunits.hh:75
G4double XTRNSpectralDensity(G4double energy)
static constexpr double fine_structure_const
Float_t x2[n_points_geant4]
Definition: compare.C:26
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4double GetGasLinearPhotoAbs(G4double)
G4PhysicsTable * fAngleForEnergyTable
Float_t X
#define DBL_MAX
Definition: templates.hh:83
G4double AngleSpectralXTRdEdx(G4double energy)
G4double XTRNSpectralAngleDensity(G4double varAngle)
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
static constexpr double GeV
Definition: G4SIunits.hh:217
const G4DynamicParticle * GetDynamicParticle() const
G4double GetElectronDensity() const
Definition: G4Material.hh:218
G4LogicalVolume * fEnvelope
void PutValue(size_t index, G4double theValue)
virtual ~G4VXTRenergyLoss()
G4VPhysicalVolume * GetVolume() const
G4complex OneInterfaceXTRdEdx(G4double energy, G4double gamma, G4double varAngle)
G4double XTRNAngleDensity(G4double varAngle)
void SetNumberOfSecondaries(G4int totSecondaries)
G4SandiaTable * fPlatePhotoAbsCof