Changeset 961 for trunk/source/processes/electromagnetic/standard/src/G4PAIModel.cc

Timestamp:

Apr 6, 2009, 12:21:12 PM (15 years ago)

Author:

garnier

Message:

update processes

File:

• trunk/source/processes/electromagnetic/standard/src/G4PAIModel.cc (modified) (9 diffs)

trunk/source/processes/electromagnetic/standard/src/G4PAIModel.cc

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4PAIModel.cc,v 1.46 2009/02/19 19:17:50 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class
 // File name:     G4PAIModel.cc
 //
@@ -177 +183 @@
       fCutCouple = theCoupleTable->GetMaterialCutsCouple( fMaterial, 
                                           curReg->GetProductionCuts() );
+      //G4cout << "Reg <" <<curReg->GetName() << ">  mat <" 
+      //             << fMaterial->GetName() << ">  fCouple= " 
+      //             << fCutCouple<<G4endl;
       if( fCutCouple ) {
         fMaterialCutsCoupleVector.push_back(fCutCouple);
@@ -197 +206 @@
   }
 }
+
+//////////////////////////////////////////////////////////////////
+
+void G4PAIModel::InitialiseMe(const G4ParticleDefinition*) 
+{}
 
 //////////////////////////////////////////////////////////////////
@@ -393 +407 @@
   {
     //  if ( x1 == x2  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
-    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    //    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*0.5 ;
     else             dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;      
   }
@@ -435 +450 @@
   {
     //  if ( x1 == x2  ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
-    if ( std::abs(x1-x2) <= eV  ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    //    if ( std::abs(x1-x2) <= eV  ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    if ( std::abs(x1-x2) <= eV  ) dEdxCut = y1 + (y2 - y1)*0.5 ;
     else             dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;      
   }
@@ -444 +460 @@
 //////////////////////////////////////////////////////////////////////////////
 
-G4double G4PAIModel::ComputeDEDX(const G4MaterialCutsCouple* matCC,
-                                 const G4ParticleDefinition* p,
-                                       G4double kineticEnergy,
-                                       G4double cutEnergy)
+G4double G4PAIModel::ComputeDEDXPerVolume(const G4Material*,
+                                          const G4ParticleDefinition* p,
+                                          G4double kineticEnergy,
+                                          G4double cutEnergy)
 {
   G4int iTkin,iPlace;
   size_t jMat;
+  
+  //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
+  G4double cut = cutEnergy;
+
   G4double massRatio  = fMass/p->GetPDGMass();
   G4double scaledTkin = kineticEnergy*massRatio;
   G4double charge     = p->GetPDGCharge();
-  G4double charge2    = charge*charge, dEdx;
+  G4double charge2    = charge*charge;
+  const G4MaterialCutsCouple* matCC = CurrentCouple();
 
   for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
@@ -470 +491 @@
   iPlace = iTkin - 1;
   if(iPlace < 0) iPlace = 0;
-  dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cutEnergy) ) ;  
-
+  G4double dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) );
   if( dEdx < 0.) dEdx = 0.;
   return dEdx;
@@ -478 +498 @@
 /////////////////////////////////////////////////////////////////////////
 
-G4double G4PAIModel::CrossSection( const G4MaterialCutsCouple* matCC,
-                                   const G4ParticleDefinition* p,
-                                         G4double kineticEnergy,
-                                         G4double cutEnergy,
-                                         G4double maxEnergy  ) 
+G4double G4PAIModel::CrossSectionPerVolume( const G4Material*,
+                                            const G4ParticleDefinition* p,
+                                            G4double kineticEnergy,
+                                            G4double cutEnergy,
+                                            G4double maxEnergy  ) 
 {
   G4int iTkin,iPlace;
   size_t jMat;
-  G4double tmax = min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
+  G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
+  if(tmax <= cutEnergy) return 0.0;
   G4double massRatio  = fMass/p->GetPDGMass();
   G4double scaledTkin = kineticEnergy*massRatio;
   G4double charge     = p->GetPDGCharge();
   G4double charge2    = charge*charge, cross, cross1, cross2;
+  const G4MaterialCutsCouple* matCC = CurrentCouple();
 
   for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
@@ -935 +957 @@
 }
 
+/////////////////////////////////////////////////////////////////////
+
+G4double G4PAIModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,
+                                         G4double kinEnergy) 
+{
+  G4double tmax = kinEnergy;
+  if(p == fElectron) tmax *= 0.5;
+  else if(p != fPositron) { 
+    G4double mass = p->GetPDGMass();
+    G4double ratio= electron_mass_c2/mass;
+    G4double gamma= kinEnergy/mass + 1.0;
+    tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
+                  (1. + 2.0*gamma*ratio + ratio*ratio);
+  }
+  return tmax;
+}
+
+///////////////////////////////////////////////////////////////
+
+void G4PAIModel::DefineForRegion(const G4Region* r) 
+{
+  fPAIRegionVector.push_back(r);
+}
 
 //