Back to home page

EIC code displayed by LXR

 
 

    


Warning, file /geant4/examples/advanced/eRosita/physics/src/G4RDAtomicDeexcitation.cc was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 //
0028 // Authors: Elena Guardincerri (Elena.Guardincerri@ge.infn.it)
0029 //          Alfonso Mantero (Alfonso.Mantero@ge.infn.it)
0030 //
0031 // History:
0032 // -----------
0033 //  
0034 //  16 Sept 2001  First committed to cvs
0035 //  12 Sep  2003  Bug in auger production fixed
0036 //
0037 // -------------------------------------------------------------------
0038 
0039 #include "G4RDAtomicDeexcitation.hh"
0040 #include "Randomize.hh"
0041 #include "G4PhysicalConstants.hh"
0042 #include "G4SystemOfUnits.hh"
0043 #include "G4Gamma.hh"
0044 #include "G4Electron.hh"
0045 #include "G4RDAtomicTransitionManager.hh"
0046 #include "G4RDFluoTransition.hh"
0047 
0048 G4RDAtomicDeexcitation::G4RDAtomicDeexcitation():
0049   minGammaEnergy(100.*eV),
0050   minElectronEnergy(100.*eV),
0051   fAuger(false)
0052 {}
0053 
0054 G4RDAtomicDeexcitation::~G4RDAtomicDeexcitation()
0055 {}
0056 
0057 std::vector<G4DynamicParticle*>* G4RDAtomicDeexcitation::GenerateParticles(G4int Z,G4int givenShellId)
0058 { 
0059 
0060   std::vector<G4DynamicParticle*>* vectorOfParticles;
0061   
0062   vectorOfParticles = new std::vector<G4DynamicParticle*>;
0063   G4DynamicParticle* aParticle;
0064   G4int provShellId = 0;
0065   G4int counter = 0;
0066   
0067   // The aim of this loop is to generate more than one fluorecence photon 
0068   // from the same ionizing event 
0069   do
0070     {
0071       if (counter == 0) 
0072     // First call to GenerateParticles(...):
0073     // givenShellId is given by the process
0074     {
0075       provShellId = SelectTypeOfTransition(Z, givenShellId);
0076       //std::cout << "AtomicDeexcitation::Generate counter 0 - provShellId = "
0077       //<< provShellId << std::endl;
0078 
0079       if  ( provShellId >0) 
0080         {
0081           aParticle = GenerateFluorescence(Z,givenShellId,provShellId);  
0082           //std::cout << "AtomicDeexcitation::Generate Fluo counter 0 " << std::endl;
0083         }
0084       else if ( provShellId == -1)
0085         {
0086           aParticle = GenerateAuger(Z, givenShellId);
0087           //std::cout << "AtomicDeexcitation::Generate Auger counter 0 " << std::endl;
0088         }
0089       else
0090         {
0091           G4Exception("G4RDAtomicDeexcitation::GenerateParticles()",
0092                           "InvalidSetup", FatalException,
0093                           "Starting shell uncorrect: check it!");
0094         }
0095     }
0096       else 
0097     // Following calls to GenerateParticles(...):
0098     // newShellId is given by GenerateFluorescence(...)
0099     {
0100       provShellId = SelectTypeOfTransition(Z,newShellId);
0101       //std::cout << "AtomicDeexcitation::Generate counter 0 - provShellId = "
0102       //<< provShellId << ", new ShellId = "<< newShellId
0103       //<< std::endl;
0104 
0105 
0106       if  (provShellId >0)
0107         {
0108           aParticle = GenerateFluorescence(Z,newShellId,provShellId);
0109           //std::cout << "AtomicDeexcitation::Generate Fluo " << std::endl;
0110         }
0111       else if ( provShellId == -1)
0112         {
0113           aParticle = GenerateAuger(Z, newShellId);
0114           //std::cout << "AtomicDeexcitation::Generate Auger " << std::endl;
0115         }
0116       else
0117         {
0118           G4Exception("G4RDAtomicDeexcitation::GenerateParticles()",
0119                           "InvalidSetup", FatalException,
0120                           "Starting shell uncorrect: check it!");
0121         }
0122     }
0123       counter++;
0124       if (aParticle != 0) {vectorOfParticles->push_back(aParticle);}
0125       else {provShellId = -2;}
0126     }
0127   
0128   // Look this in a particular way: only one auger emitted! //
0129   while (provShellId > -2); 
0130   
0131   return vectorOfParticles;
0132 }
0133 
0134 G4int G4RDAtomicDeexcitation::SelectTypeOfTransition(G4int Z, G4int shellId)
0135 {
0136   if (shellId <=0 ) 
0137     {
0138       G4Exception("G4RDAtomicDeexcitation::SelectTypeOfTransition()",
0139                   "InvalidCondition", FatalException,
0140                   "Zero or negative shellId!");
0141     }
0142 
0143   const G4RDAtomicTransitionManager*  transitionManager = 
0144         G4RDAtomicTransitionManager::Instance();
0145   G4int provShellId = -1;
0146   G4int shellNum = 0;
0147   G4int maxNumOfShells = transitionManager->NumberOfReachableShells(Z);  
0148   
0149   //std::cout << "AtomicDeexcitation::SelectType -  NumberOfReachableShells = "
0150   //<< maxNumOfShells<< std::endl;
0151 
0152   const G4RDFluoTransition* refShell = transitionManager->ReachableShell(Z,maxNumOfShells-1);
0153 
0154   // This loop gives shellNum the value of the index of shellId
0155   // in the vector storing the list of the shells reachable through
0156   // a radiative transition
0157   if ( shellId <= refShell->FinalShellId())
0158     {
0159       while (shellId != transitionManager->ReachableShell(Z,shellNum)->FinalShellId())
0160     {
0161       if(shellNum ==maxNumOfShells-1)
0162         {
0163           break;
0164         }
0165       shellNum++;
0166     }
0167       G4int transProb = 0; //AM change 29/6/07 was 1
0168    
0169       G4double partialProb = G4UniformRand();      
0170       G4double partSum = 0;
0171       const G4RDFluoTransition* aShell = transitionManager->ReachableShell(Z,shellNum);      
0172       G4int trSize =  (aShell->TransitionProbabilities()).size();
0173     
0174       // Loop over the shells wich can provide an electron for a 
0175       // radiative transition towards shellId:
0176       // in every loop the partial sum of the first transProb shells
0177       // is calculated and compared with a random number [0,1].
0178       // If the partial sum is greater, the shell whose index is transProb
0179       // is chosen as the starting shell for a radiative transition
0180       // and its identity is returned
0181       // Else, terminateded the loop, -1 is returned
0182       while(transProb < trSize){
0183     
0184      partSum += aShell->TransitionProbability(transProb);
0185 
0186      if(partialProb <= partSum)
0187        {
0188          provShellId = aShell->OriginatingShellId(transProb);
0189          break;
0190        }
0191      transProb++;
0192       }
0193 
0194       // here provShellId is the right one or is -1.
0195       // if -1, the control is passed to the Auger generation part of the package 
0196     }
0197 
0198 
0199 
0200   else 
0201     {
0202  
0203      provShellId = -1;
0204 
0205     }
0206   return provShellId;
0207 }
0208 
0209 G4DynamicParticle* G4RDAtomicDeexcitation::GenerateFluorescence(G4int Z, 
0210                                   G4int shellId,
0211                                   G4int provShellId )
0212 { 
0213 
0214 
0215   const G4RDAtomicTransitionManager*  transitionManager = G4RDAtomicTransitionManager::Instance();
0216   //  G4int provenienceShell = provShellId;
0217 
0218   //isotropic angular distribution for the outcoming photon
0219   G4double newcosTh = 1.-2.*G4UniformRand();
0220   G4double  newsinTh = std::sqrt(1.-newcosTh*newcosTh);
0221   G4double newPhi = twopi*G4UniformRand();
0222   
0223   G4double xDir =  newsinTh*std::sin(newPhi);
0224   G4double yDir = newsinTh*std::cos(newPhi);
0225   G4double zDir = newcosTh;
0226   
0227   G4ThreeVector newGammaDirection(xDir,yDir,zDir);
0228   
0229   G4int shellNum = 0;
0230   G4int maxNumOfShells = transitionManager->NumberOfReachableShells(Z);
0231   
0232   // find the index of the shell named shellId
0233   while (shellId != transitionManager->
0234      ReachableShell(Z,shellNum)->FinalShellId())
0235     {
0236       if(shellNum == maxNumOfShells-1)
0237     {
0238       break;
0239     }
0240       shellNum++;
0241     }
0242   // number of shell from wich an electron can reach shellId
0243   size_t transitionSize = transitionManager->
0244     ReachableShell(Z,shellNum)->OriginatingShellIds().size();
0245   
0246   size_t index = 0;
0247   
0248   // find the index of the shell named provShellId in the vector
0249   // storing the shells from which shellId can be reached 
0250   while (provShellId != transitionManager->
0251      ReachableShell(Z,shellNum)->OriginatingShellId(index))
0252     {
0253       if(index ==  transitionSize-1)
0254     {
0255       break;
0256     }
0257       index++;
0258     }
0259   // energy of the gamma leaving provShellId for shellId
0260   G4double transitionEnergy = transitionManager->
0261     ReachableShell(Z,shellNum)->TransitionEnergy(index);
0262   
0263   // This is the shell where the new vacancy is: it is the same
0264   // shell where the electron came from
0265   newShellId = transitionManager->
0266     ReachableShell(Z,shellNum)->OriginatingShellId(index);
0267   
0268   
0269   G4DynamicParticle* newPart = new G4DynamicParticle(G4Gamma::Gamma(), 
0270                              newGammaDirection,
0271                              transitionEnergy);
0272   return newPart;
0273 }
0274 
0275 G4DynamicParticle* G4RDAtomicDeexcitation::GenerateAuger(G4int Z, G4int shellId)
0276 {
0277   if(!fAuger) return 0;
0278   
0279 
0280   const G4RDAtomicTransitionManager*  transitionManager = 
0281         G4RDAtomicTransitionManager::Instance();
0282 
0283 
0284 
0285   if (shellId <=0 ) 
0286     {
0287       G4Exception("G4RDAtomicDeexcitation::GenerateAuger()",
0288                   "InvalidCondition", FatalException,
0289                   "Zero or negative shellId!");
0290     }
0291   
0292   // G4int provShellId = -1;
0293   G4int maxNumOfShells = transitionManager->NumberOfReachableAugerShells(Z);  
0294   
0295   const G4RDAugerTransition* refAugerTransition = 
0296         transitionManager->ReachableAugerShell(Z,maxNumOfShells-1);
0297 
0298 
0299   // This loop gives to shellNum the value of the index of shellId
0300   // in the vector storing the list of the vacancies in the variuos shells 
0301   // that can originate a NON-radiative transition
0302   
0303   // ---- MGP ---- Next line commented out to remove compilation warning
0304   // G4int p = refAugerTransition->FinalShellId();
0305 
0306   G4int shellNum = 0;
0307 
0308 
0309   if ( shellId <= refAugerTransition->FinalShellId() ) 
0310     //"FinalShellId" is final from the point of view of the elctron who makes the transition, 
0311     // being the Id of the shell in which there is a vacancy
0312     {
0313       G4int pippo = transitionManager->ReachableAugerShell(Z,shellNum)->FinalShellId();
0314       if (shellId  != pippo ) {
0315     do { 
0316       shellNum++;
0317       if(shellNum == maxNumOfShells)
0318         {
0319 //            G4cout << "G4RDAtomicDeexcitation warning: No Auger transition found" <<  G4endl;
0320 //        G4cout << "Absorbed enrgy deposited locally" << G4endl;
0321           return 0;
0322 //        //  G4Exception("G4RDAtomicDeexcitation: No Auger transition found");
0323         }
0324     }
0325     while (shellId != (transitionManager->ReachableAugerShell(Z,shellNum)->FinalShellId()) ) ;
0326       }
0327       /*    {
0328 
0329       if(shellNum == maxNumOfShells-1)
0330         {
0331           G4Exception("G4RDAtomicDeexcitation: No Auger tramsition found");
0332         }
0333       shellNum++;
0334       }*/
0335     
0336 
0337 
0338 
0339       // Now we have that shellnum is the shellIndex of the shell named ShellId
0340 
0341       //      G4cout << " the index of the shell is: "<<shellNum<<G4endl;
0342 
0343       // But we have now to select two shells: one for the transition, 
0344       // and another for the auger emission.
0345 
0346       G4int transitionLoopShellIndex = 0;      
0347       G4double partSum = 0;
0348       const G4RDAugerTransition* anAugerTransition = 
0349             transitionManager->ReachableAugerShell(Z,shellNum);
0350 
0351       //      G4cout << " corresponding to the ID: "<< anAugerTransition->FinalShellId() << G4endl;
0352 
0353 
0354       G4int transitionSize = 
0355             (anAugerTransition->TransitionOriginatingShellIds())->size();
0356       while (transitionLoopShellIndex < transitionSize) {
0357 
0358         std::vector<G4int>::const_iterator pos = 
0359                anAugerTransition->TransitionOriginatingShellIds()->begin();
0360 
0361         G4int transitionLoopShellId = *(pos+transitionLoopShellIndex);
0362         G4int numberOfPossibleAuger = 
0363               (anAugerTransition->AugerTransitionProbabilities(transitionLoopShellId))->size();
0364         G4int augerIndex = 0;
0365         //      G4int partSum2 = 0;
0366 
0367 
0368     if (augerIndex < numberOfPossibleAuger) {
0369       
0370       do 
0371         {
0372           G4double thisProb = anAugerTransition->AugerTransitionProbability(augerIndex, 
0373                                         transitionLoopShellId);
0374           partSum += thisProb;
0375           augerIndex++;
0376           
0377         } while (augerIndex < numberOfPossibleAuger);
0378         }
0379         transitionLoopShellIndex++;
0380       }
0381       
0382 
0383 
0384       // Now we have the entire probability of an auger transition for the vacancy 
0385       // located in shellNum (index of shellId) 
0386 
0387       // AM *********************** F I X E D **************************** AM
0388       // Here we duplicate the previous loop, this time looking to the sum of the probabilities 
0389       // to be under the random number shoot by G4 UniformRdandom. This could have been done in the 
0390       // previuos loop, while integrating the probabilities. There is a bug that will be fixed 
0391       // 5 minutes from now: a line:
0392       // G4int numberOfPossibleAuger = (anAugerTransition->
0393       // AugerTransitionProbabilities(transitionLoopShellId))->size();
0394       // to be inserted.
0395       // AM *********************** F I X E D **************************** AM
0396 
0397       // Remains to get the same result with a single loop.
0398 
0399       // AM *********************** F I X E D **************************** AM
0400       // Another Bug: in EADL Auger Transition are normalized to all the transitions deriving from 
0401       // a vacancy in one shell, but not all of these are present in data tables. So if a transition 
0402       // doesn't occur in the main one a local energy deposition must occur, instead of (like now) 
0403       // generating the last transition present in EADL data.
0404       // AM *********************** F I X E D **************************** AM
0405 
0406 
0407       G4double totalVacancyAugerProbability = partSum;
0408 
0409 
0410       //And now we start to select the right auger transition and emission
0411       G4int transitionRandomShellIndex = 0;
0412       G4int transitionRandomShellId = 1;
0413       G4int augerIndex = 0;
0414       partSum = 0; 
0415       G4double partialProb = G4UniformRand();
0416       // G4int augerOriginatingShellId = 0;
0417       
0418       G4int numberOfPossibleAuger = 
0419       (anAugerTransition->AugerTransitionProbabilities(transitionRandomShellId))->size();
0420       G4bool foundFlag = false;
0421 
0422       while (transitionRandomShellIndex < transitionSize) {
0423 
0424         std::vector<G4int>::const_iterator pos = 
0425                anAugerTransition->TransitionOriginatingShellIds()->begin();
0426 
0427         transitionRandomShellId = *(pos+transitionRandomShellIndex);
0428         
0429     augerIndex = 0;
0430     numberOfPossibleAuger = (anAugerTransition-> 
0431                  AugerTransitionProbabilities(transitionRandomShellId))->size();
0432 
0433         while (augerIndex < numberOfPossibleAuger) {
0434       G4double thisProb =anAugerTransition->AugerTransitionProbability(augerIndex, 
0435                                        transitionRandomShellId);
0436 
0437           partSum += thisProb;
0438           
0439           if (partSum >= (partialProb*totalVacancyAugerProbability) ) { // was /
0440         foundFlag = true;
0441         break;
0442       }
0443           augerIndex++;
0444         }
0445         if (partSum >= (partialProb*totalVacancyAugerProbability) ) {break;} // was /
0446         transitionRandomShellIndex++;
0447       }
0448 
0449       // Now we have the index of the shell from wich comes the auger electron (augerIndex), 
0450       // and the id of the shell, from which the transition e- come (transitionRandomShellid)
0451       // If no Transition has been found, 0 is returned.  
0452 
0453       if (!foundFlag) {return 0;}      
0454       
0455       // Isotropic angular distribution for the outcoming e-
0456       G4double newcosTh = 1.-2.*G4UniformRand();
0457       G4double  newsinTh = std::sqrt(1.-newcosTh*newcosTh);
0458       G4double newPhi = twopi*G4UniformRand();
0459       
0460       G4double xDir =  newsinTh*std::sin(newPhi);
0461       G4double yDir = newsinTh*std::cos(newPhi);
0462       G4double zDir = newcosTh;
0463       
0464       G4ThreeVector newElectronDirection(xDir,yDir,zDir);
0465       
0466       // energy of the auger electron emitted
0467       
0468       
0469       G4double transitionEnergy = anAugerTransition->AugerTransitionEnergy(augerIndex, transitionRandomShellId);
0470       /*
0471     G4cout << "AUger TransitionId " << anAugerTransition->FinalShellId() << G4endl;
0472     G4cout << "augerIndex: " << augerIndex << G4endl;
0473     G4cout << "transitionShellId: " << transitionRandomShellId << G4endl;
0474       */
0475       
0476       // This is the shell where the new vacancy is: it is the same
0477       // shell where the electron came from
0478       newShellId = transitionRandomShellId;
0479       
0480       
0481       G4DynamicParticle* newPart = new G4DynamicParticle(G4Electron::Electron(), 
0482                              newElectronDirection,
0483                              transitionEnergy);
0484       return newPart;
0485 
0486     }
0487   else 
0488     {
0489       //G4Exception("G4RDAtomicDeexcitation: no auger transition found");
0490       return 0;
0491     }
0492   
0493 }
0494 
0495 void G4RDAtomicDeexcitation::SetCutForSecondaryPhotons(G4double cut)
0496 {
0497   minGammaEnergy = cut;
0498 }
0499 
0500 void G4RDAtomicDeexcitation::SetCutForAugerElectrons(G4double cut)
0501 {
0502   minElectronEnergy = cut;
0503 }
0504 
0505 void G4RDAtomicDeexcitation::ActivateAugerElectronProduction(G4bool val)
0506 {
0507   fAuger = val;
0508 }
0509 
0510 
0511 
0512 
0513 
0514 
0515