G4PiMinusAbsorptionAtRest.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //      File name:     G4PiMinusAbsorptionAtRest.hh 
00027 //
00028 //      Author:        Maria Grazia Pia (pia@genova.infn.it)
00029 // 
00030 //      Creation date: 8 May 1998
00031 //
00032 //      Modifications: 
00033 //      MGP            4 Jul 1998 Changed excitation energy calculation
00034 //      MGP           14 Sep 1998 Fixed excitation energy calculation
00035 //
00036 // -------------------------------------------------------------------
00037 
00038 #include "G4ios.hh"
00039 
00040 #include "G4PiMinusAbsorptionAtRest.hh"
00041 
00042 #include "G4SystemOfUnits.hh"
00043 #include "G4PiMinusStopLi.hh"
00044 #include "G4PiMinusStopC.hh"
00045 #include "G4PiMinusStopN.hh"
00046 #include "G4PiMinusStopO.hh"
00047 #include "G4PiMinusStopAl.hh"
00048 #include "G4PiMinusStopCu.hh"
00049 #include "G4PiMinusStopCo.hh"
00050 #include "G4PiMinusStopTa.hh"
00051 #include "G4PiMinusStopPb.hh"
00052 #include "G4StopTheoDeexcitation.hh"
00053 #include "G4StopDummyDeexcitation.hh"
00054 #include "G4DynamicParticle.hh"
00055 #include "G4DynamicParticleVector.hh"
00056 #include "Randomize.hh"
00057 #include "G4ThreeVector.hh"
00058 #include "G4LorentzVector.hh"
00059 #include "G4HadronicProcessStore.hh"
00060 #include "G4HadronicDeprecate.hh"
00061 
00062 
00063 // Constructor
00064 
00065 G4PiMinusAbsorptionAtRest::G4PiMinusAbsorptionAtRest(const G4String& processName,
00066                                                      G4ProcessType aType) :
00067   G4VRestProcess (processName, aType)
00068 {
00069   G4HadronicDeprecate("G4PiMinusAbsorptionAtRest");
00070 
00071   SetProcessSubType(fHadronAtRest);
00072 
00073   _indexDeexcitation = 0;
00074 
00075   if (verboseLevel>0) 
00076     { G4cout << GetProcessName() << " is created "<< G4endl; }
00077   G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
00078 }
00079 
00080 
00081 // Destructor
00082 
00083 G4PiMinusAbsorptionAtRest::~G4PiMinusAbsorptionAtRest()
00084 {
00085   G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
00086 }
00087 
00088 void G4PiMinusAbsorptionAtRest::PreparePhysicsTable(const G4ParticleDefinition& p) 
00089 {
00090   G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p);
00091 }
00092 
00093 void G4PiMinusAbsorptionAtRest::BuildPhysicsTable(const G4ParticleDefinition& p) 
00094 {
00095   G4HadronicProcessStore::Instance()->PrintInfo(&p);
00096 }
00097 
00098 G4VParticleChange* G4PiMinusAbsorptionAtRest::AtRestDoIt(const G4Track& track, const G4Step& )
00099 {
00100   const G4DynamicParticle* stoppedHadron = track.GetDynamicParticle();
00101       
00102   // Check applicability
00103   if (! IsApplicable(*(stoppedHadron->GetDefinition())))
00104     {
00105       G4cerr  
00106       << "G4PiMinusAbsorptionAtRest: ERROR, particle must be a pion minus!" 
00107       << G4endl;
00108       return NULL;
00109     }
00110 
00111   // Get the current material
00112   const G4Material* material = track.GetMaterial();
00113 
00114   G4double A=-1;
00115   G4double Z=-1;
00116   G4double random = G4UniformRand();
00117   const G4ElementVector* theElementVector = material->GetElementVector();
00118   unsigned int i;
00119   G4double sum = 0;
00120   G4double totalsum=0;
00121   for(i=0; i<material->GetNumberOfElements(); ++i)
00122     {
00123       if((*theElementVector)[i]->GetZ()!=1) totalsum+=material->GetFractionVector()[i];
00124     }
00125   for (i = 0; i<material->GetNumberOfElements(); ++i)
00126     {
00127       if((*theElementVector)[i]->GetZ()!=1) sum += material->GetFractionVector()[i];
00128       if ( sum/totalsum > random )  
00129         { 
00130           A = (*theElementVector)[i]->GetA()*mole/g;
00131           Z = (*theElementVector)[i]->GetZ();
00132           break;
00133         }
00134     }
00135 
00136   // Do the interaction with the nucleon cluster
00137   
00138   G4PiMinusStopMaterial* algorithm = LoadAlgorithm(static_cast<G4int>(Z));
00139   G4PiMinusStopAbsorption stopAbsorption(algorithm,Z,A);
00140   stopAbsorption.SetVerboseLevel(verboseLevel);
00141 
00142   G4DynamicParticleVector* absorptionProducts = stopAbsorption.DoAbsorption();
00143 
00144   // Deal with the leftover nucleus
00145 
00146   G4double pionEnergy = stoppedHadron->GetTotalEnergy();
00147   G4double excitation = pionEnergy - stopAbsorption.Energy();
00148   if (excitation < 0.) 
00149   {
00150     G4Exception("G4PiMinusAbsorptionAtRest::AtRestDoIt()", "HAD_STOP_0000",
00151                 FatalException, "Excitation energy < 0");
00152   }
00153   if (verboseLevel>0) { G4cout << " excitation " << excitation << G4endl; }
00154 
00155   G4StopDeexcitationAlgorithm* nucleusAlgorithm = LoadNucleusAlgorithm();
00156   G4StopDeexcitation stopDeexcitation(nucleusAlgorithm);
00157 
00158   G4double newZ = Z - stopAbsorption.NProtons();
00159   G4double newN = A - Z - stopAbsorption.NNeutrons();
00160   G4double newA = newZ + newN;
00161   G4ReactionProductVector* fragmentationProducts = stopDeexcitation.DoBreakUp(newA,newZ,excitation,stopAbsorption.RecoilMomentum());
00162 
00163   unsigned int nAbsorptionProducts = 0;
00164   if (absorptionProducts != 0)     
00165     { nAbsorptionProducts  =  absorptionProducts->size(); }
00166 
00167   unsigned int nFragmentationProducts = 0;
00168   if (fragmentationProducts != 0) 
00169     { nFragmentationProducts = fragmentationProducts->size(); }
00170   
00171   if (verboseLevel>0) 
00172     {
00173       G4cout << "nAbsorptionProducts = " << nAbsorptionProducts
00174              << "  nFragmentationProducts = " << nFragmentationProducts
00175              << G4endl;
00176     }
00177 
00178   // Deal with ParticleChange final state
00179 
00180   aParticleChange.Initialize(track);
00181   aParticleChange.SetNumberOfSecondaries(G4int(nAbsorptionProducts + nFragmentationProducts)); 
00182      
00183   for (i = 0; i<nAbsorptionProducts; i++)
00184     { aParticleChange.AddSecondary((*absorptionProducts)[i]); }
00185 
00186 //  for (i = 0; i<nFragmentationProducts; i++)
00187 //    { aParticleChange.AddSecondary(fragmentationProducts->at(i)); }
00188   for(i=0; i<nFragmentationProducts; i++)
00189   {
00190     G4DynamicParticle * aNew = 
00191        new G4DynamicParticle((*fragmentationProducts)[i]->GetDefinition(),
00192                              (*fragmentationProducts)[i]->GetMomentum());
00193     G4double newTime = aParticleChange.GetGlobalTime((*fragmentationProducts)[i]->GetFormationTime());
00194     aParticleChange.AddSecondary(aNew, newTime);
00195     delete (*fragmentationProducts)[i];
00196   }
00197 
00198   if (fragmentationProducts != 0)   delete fragmentationProducts;   
00199 
00200   if (_indexDeexcitation == 1) aParticleChange.ProposeLocalEnergyDeposit(excitation);
00201 
00202   // Kill the absorbed pion
00203   aParticleChange.ProposeTrackStatus(fStopAndKill); 
00204 
00205   return &aParticleChange;
00206 
00207 }
00208 
00209 G4PiMinusStopMaterial* G4PiMinusAbsorptionAtRest::LoadAlgorithm(int Z)
00210 {  
00211   if (verboseLevel>0) 
00212     {
00213       G4cout << "Load material algorithm " << Z << G4endl; 
00214     }
00215 
00216   G4int index = 0;
00217   if (Z > 0 && Z < 4) {index = 3;}
00218   if (Z > 3 && Z < 7) {index = 6;}
00219   if (Z == 7) {index = 7;}
00220   if (Z >= 8 && Z<= 11) {index = 8;}
00221   if (Z >= 12 && Z<= 18) {index = 13;}
00222   if (Z >=19 && Z<= 27) {index = 27;}
00223   if (Z >= 28 && Z<= 51) {index = 29;}
00224   if (Z >=52 ) {index = 73;}
00225 
00226   switch (index) 
00227     {
00228     case 3:
00229       if (verboseLevel>0) 
00230         { G4cout << " =================== Load Li algorithm " << G4endl; }
00231       return new G4PiMinusStopLi();
00232     case 6:
00233       if (verboseLevel>0) 
00234         { G4cout << " =================== Load C algorithm " << G4endl; }
00235       return new G4PiMinusStopC();
00236     case 7:
00237       if (verboseLevel>0) 
00238         { G4cout << " =================== Load N algorithm " << G4endl; }
00239       return new G4PiMinusStopN();
00240     case 8:
00241       if (verboseLevel>0) 
00242         { G4cout << " =================== Load O algorithm " << G4endl; }
00243       return new G4PiMinusStopO();
00244     case 13:
00245       if (verboseLevel>0) 
00246         { G4cout << " =================== Load Al algorithm " << G4endl; }
00247       return new G4PiMinusStopAl();
00248     case 27:
00249       if (verboseLevel>0) 
00250         { G4cout << " =================== Load Cu algorithm " << G4endl; }
00251       return new G4PiMinusStopCu();
00252     case 29:
00253       if (verboseLevel>0) 
00254         { G4cout << " =================== Load Co algorithm " << G4endl; }
00255       return new G4PiMinusStopCo();
00256     case 73:
00257       if (verboseLevel>0) 
00258         { G4cout << " =================== Load Ta algorithm " << G4endl; }
00259       return new G4PiMinusStopTa();
00260     default: 
00261       if (verboseLevel>0) 
00262         { G4cout << " =================== Load default material algorithm " << G4endl; }
00263       return new G4PiMinusStopC();
00264     }
00265 }
00266 
00267 G4StopDeexcitationAlgorithm* G4PiMinusAbsorptionAtRest::LoadNucleusAlgorithm()
00268 {  
00269 
00270   switch (_indexDeexcitation) 
00271     {
00272     case 0:
00273       if (verboseLevel>0) 
00274         { G4cout << " =================== Load Theo deexcitation " << G4endl; }
00275       return new G4StopTheoDeexcitation();
00276     case 1:
00277       if (verboseLevel>0) 
00278         { G4cout << " =================== Load Dummy deexcitation " << G4endl; }
00279       return new G4StopDummyDeexcitation();
00280     default:  
00281       if (verboseLevel>0) 
00282         { G4cout << " =================== Load default deexcitation " << G4endl; }
00283       return new G4StopTheoDeexcitation();
00284     }
00285 }
00286 
00287 void G4PiMinusAbsorptionAtRest::SetDeexcitationAlgorithm(G4int index)
00288 {  
00289   _indexDeexcitation = index;
00290 }

Generated on Mon May 27 17:49:20 2013 for Geant4 by  doxygen 1.4.7