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