Ignore:
Timestamp:
Nov 5, 2010, 3:45:55 PM (14 years ago)
Author:
garnier
Message:

update ti head

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/binary_cascade/src/G4GeneratorPrecompoundInterface.cc

    r1196 r1340  
    2424// ********************************************************************
    2525//
     26// $Id: G4GeneratorPrecompoundInterface.cc,v 1.10 2010/08/31 16:16:51 vnivanch Exp $
     27// GEANT4 tag $Name: had-binary-V09-03-03 $
     28//
     29// -----------------------------------------------------------------------------
     30//      GEANT 4 class file
     31//
     32//      History: first implementation
     33//      HPW, 10DEC 98, the decay part originally written by Gunter Folger
     34//                in his FTF-test-program.
     35//
     36//
     37// -----------------------------------------------------------------------------
     38
    2639#include "G4GeneratorPrecompoundInterface.hh"
    2740#include "G4DynamicParticleVector.hh"
    28 #include "G4IonTable.hh"
     41#include "G4KineticTrackVector.hh"
     42#include "G4Proton.hh"
     43#include "G4Neutron.hh"
     44#include "G4V3DNucleus.hh"
     45#include "G4Nucleon.hh"
     46#include "G4FragmentVector.hh"
     47#include "G4ReactionProduct.hh"
     48#include "G4PreCompoundModel.hh"
     49#include "G4ExcitationHandler.hh"
    2950
    30 //
    31 // HPW, 10DEC 98, the decay part originally written by Gunter Folger in his FTF-test-program.
    32 //
    33    
    34     G4GeneratorPrecompoundInterface::G4GeneratorPrecompoundInterface()
    35     : CaptureThreshold(80*MeV)
    36     {}
    37      
    38    G4HadFinalState* G4GeneratorPrecompoundInterface::
    39    ApplyYourself(const G4HadProjectile &, G4Nucleus & )
    40    {
    41      std::cout << "G4GeneratorPrecompoundInterface: ApplyYourself interface called stand-allone."<< G4endl;
    42      std::cout << "This class is only a mediator between generator and precompound"<<G4endl;
    43      std::cout << "Please remove from your physics list."<<G4endl;
    44      throw G4HadronicException(__FILE__, __LINE__, "SEVERE: G4GeneratorPrecompoundInterface model interface called stand-allone.");
    45      return new G4HadFinalState;
    46    }
    47    
    48    G4ReactionProductVector* G4GeneratorPrecompoundInterface::
    49    Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
    50    {
    51      G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
     51G4GeneratorPrecompoundInterface::G4GeneratorPrecompoundInterface(G4VPreCompoundModel* p)
     52  : CaptureThreshold(80*MeV)
     53{
     54  proton = G4Proton::Proton();
     55  neutron = G4Neutron::Neutron();
     56  if(p) { SetDeExcitation(p); }
     57  else  { SetDeExcitation(new G4PreCompoundModel(new G4ExcitationHandler())); }
     58}
     59         
     60G4GeneratorPrecompoundInterface::~G4GeneratorPrecompoundInterface()
     61{}
     62         
     63G4ReactionProductVector* G4GeneratorPrecompoundInterface::
     64Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
     65{
     66  G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
    5267
    53      // decay the strong resonances
    54      G4KineticTrackVector *result1, *secondaries, *result;
    55      result1=theSecondaries;
    56      result=new G4KineticTrackVector();
    57 
    58      for (unsigned int aResult=0; aResult < result1->size(); aResult++)
    59      {
    60        G4ParticleDefinition * pdef;
    61        pdef=result1->operator[](aResult)->GetDefinition();
    62        secondaries=NULL;
    63        if ( pdef->IsShortLived() )
    64        {
     68  // decay the strong resonances
     69  G4KineticTrackVector *result1, *secondaries, *result;
     70  result1=theSecondaries;
     71  result=new G4KineticTrackVector();
     72  //G4cout << "### G4GeneratorPrecompoundInterface::Propagate "
     73  //     << result1->size() << " tracks " << theDeExcitation << G4endl;
     74  for (unsigned int aResult=0; aResult < result1->size(); ++aResult)
     75    {
     76      G4ParticleDefinition * pdef;
     77      pdef=result1->operator[](aResult)->GetDefinition();
     78      secondaries=0;
     79      if ( pdef->IsShortLived() )
     80        {
    6581          secondaries = result1->operator[](aResult)->Decay();
    66        }
    67        if ( secondaries == NULL )
    68        {
     82        }
     83      if ( 0 == secondaries )
     84        {
    6985          result->push_back(result1->operator[](aResult));
    7086          result1->operator[](aResult)=NULL;    //protect for clearAndDestroy
    71        }
    72        else
    73        {
    74          for (unsigned int aSecondary=0; aSecondary<secondaries->size(); aSecondary++)
     87        }
     88      else
     89        {
     90          unsigned int amax = secondaries->size();
     91          for (unsigned int aSecondary=0; aSecondary<amax; ++aSecondary)
     92            {
     93              result1->push_back(secondaries->operator[](aSecondary));
     94            }
     95          delete secondaries;
     96        }
     97    }
     98  //G4cout << "Delete tracks" << G4endl;
     99  std::for_each(result1->begin(), result1->end(), DeleteKineticTrack());
     100  delete result1;
     101     
     102  // prepare the fragment
     103  G4int anA=theNucleus->GetMassNumber();
     104  G4int aZ=theNucleus->GetCharge();
     105  G4int numberOfEx = 0;
     106  G4int numberOfCh = 0;
     107  G4int numberOfHoles = 0;
     108  G4double exEnergy = 0.0;
     109  G4double R = theNucleus->GetNuclearRadius();
     110  G4ThreeVector exciton3Momentum(0.,0.,0.);
     111
     112  // loop over secondaries
     113  unsigned int amax = result->size();
     114  for(unsigned int list=0; list<amax; ++list)
     115    {
     116      G4KineticTrack *aTrack = result->operator[](list);
     117      G4ParticleDefinition* part = aTrack->GetDefinition();
     118      G4double e = aTrack->Get4Momentum().e();
     119      G4double mass = aTrack->Get4Momentum().mag();
     120      G4ThreeVector mom = aTrack->Get4Momentum().vect();
     121      if((part != proton && part != neutron) ||
     122         (e > mass + CaptureThreshold) ||
     123         (aTrack->GetPosition().mag() > R))
     124        {
     125          G4ReactionProduct * theNew = new G4ReactionProduct(part);
     126          theNew->SetMomentum(mom);
     127          theNew->SetTotalEnergy(e);
     128          theTotalResult->push_back(theNew);           
     129        }
     130      else
     131        {
     132          // within the nucleus, neutron or proton
     133          // now calculate  A, Z of the fragment, momentum, number of exciton states
     134          ++anA;
     135          ++numberOfEx;
     136          G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
     137          aZ += Z;
     138          numberOfCh += Z;
     139          exciton3Momentum += mom;
     140          exEnergy += (e - mass);
     141        }
     142    }
     143     
     144  // loop over wounded nucleus
     145  G4Nucleon * theCurrentNucleon =
     146    theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
     147  while(0 != theCurrentNucleon)
     148    {
     149      if(theCurrentNucleon->AreYouHit())
     150        {
     151          ++numberOfHoles;
     152          ++numberOfEx;
     153          --anA;
     154          aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/eplus + 0.1);
     155          exciton3Momentum -= theCurrentNucleon->Get4Momentum().vect();
     156          exEnergy += theCurrentNucleon->GetBindingEnergy();
     157        }
     158      theCurrentNucleon = theNucleus->GetNextNucleon();
     159    }   
     160 
     161  if(0!=anA && 0!=aZ)
     162    {
     163      G4double fMass =  G4NucleiProperties::GetNuclearMass(anA, aZ);
     164      fMass += exEnergy;
     165
     166      G4LorentzVector exciton4Momentum(exciton3Momentum,
     167                                       std::sqrt(exciton3Momentum.mag2() + fMass*fMass));
     168   
     169      G4Fragment anInitialState(anA, aZ, exciton4Momentum);
     170      anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
     171      anInitialState.SetNumberOfCharged(numberOfCh);
     172      anInitialState.SetNumberOfHoles(numberOfHoles);
     173      G4ReactionProductVector * aPreResult = theDeExcitation->DeExcite(anInitialState);
     174
     175      // fill pre-compound part into the result, and return
     176      unsigned int amax = aPreResult->size();
     177       for(unsigned int ll=0; ll<amax; ++ll)
    75178         {
    76              result1->push_back(secondaries->operator[](aSecondary));
     179           theTotalResult->push_back(aPreResult->operator[](ll));
    77180         }
    78          delete secondaries;
    79        }
    80      }
    81      std::for_each(result1->begin(), result1->end(), DeleteKineticTrack());
    82      delete result1;
     181       delete aPreResult;
     182    }
    83183     
    84      
    85      
    86      // prepare the fragment
    87      G4Fragment anInitialState;
    88      G4int anA=theNucleus->GetMassNumber();
    89      G4int aZ=theNucleus->GetCharge();
    90      G4int numberOfEx = 0;
    91      G4int numberOfCh = 0;
    92      G4int numberOfHoles = 0;
    93      G4double exEnergy = 0;
    94      G4ThreeVector exciton3Momentum(0,0,0);
    95      // loop over secondaries
    96      for(unsigned int list=0; list < result->size(); list++)
    97      {
    98        G4KineticTrack *aTrack = result->operator[](list);
    99        if(aTrack->GetDefinition() != G4Proton::Proton() &&
    100           aTrack->GetDefinition() != G4Neutron::Neutron())
    101        {
    102          G4ReactionProduct * theNew = new G4ReactionProduct(aTrack->GetDefinition());
    103          theNew->SetMomentum(aTrack->Get4Momentum().vect());
    104          theNew->SetTotalEnergy(aTrack->Get4Momentum().e());
    105          theTotalResult->push_back(theNew);           
    106        }
    107        else if(aTrack->Get4Momentum().t() - aTrack->Get4Momentum().mag()>CaptureThreshold)
    108        {
    109          G4ReactionProduct * theNew = new G4ReactionProduct(aTrack->GetDefinition());
    110          theNew->SetMomentum(aTrack->Get4Momentum().vect());
    111          theNew->SetTotalEnergy(aTrack->Get4Momentum().e());
    112          theTotalResult->push_back(theNew);           
    113        }
    114        else if(aTrack->GetPosition().mag() > theNucleus->GetNuclearRadius())
    115        {
    116          G4ReactionProduct * theNew = new G4ReactionProduct(aTrack->GetDefinition());
    117          theNew->SetMomentum(aTrack->Get4Momentum().vect());
    118          theNew->SetTotalEnergy(aTrack->Get4Momentum().e());
    119          theTotalResult->push_back(theNew);           
    120        }
    121        else
    122        {
    123          // within the nucleus, neutron or proton
    124          // now calculate  A, Z of the fragment, momentum, number of exciton states
    125          anA++;;
    126          numberOfEx++;
    127          aZ += G4int(aTrack->GetDefinition()->GetPDGCharge());
    128          numberOfCh += G4int(aTrack->GetDefinition()->GetPDGCharge());
    129          exciton3Momentum += aTrack->Get4Momentum().vect();
    130          exEnergy += (aTrack->Get4Momentum().t()-aTrack->Get4Momentum().m());
    131        }
    132      }
    133      
    134      // loop over wounded nucleus
    135      G4Nucleon * theCurrentNucleon = theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : NULL;
    136      while(theCurrentNucleon != NULL)
    137      {
    138        if(theCurrentNucleon->AreYouHit())
    139        {
    140          numberOfHoles++;
    141          numberOfEx++;
    142          anA--;
    143          aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge());
    144          exciton3Momentum -= theCurrentNucleon->Get4Momentum().vect();
    145          exEnergy+=theCurrentNucleon->GetBindingEnergy();
    146        }
    147        theCurrentNucleon = theNucleus->GetNextNucleon();
    148      }   
    149  
    150      if(!theDeExcitation)
    151      {
    152        // throw G4HadronicException(__FILE__, __LINE__, "Please register an evaporation phase with G4GeneratorPrecompoundInterface.");
    153      }
    154      else if(0!=anA && 0!=aZ)
    155      {
    156        G4double residualMass = 
    157                 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(aZ ,anA);
    158          residualMass += exEnergy;
    159 
    160        G4LorentzVector exciton4Momentum(exciton3Momentum,
    161                                         std::sqrt(exciton3Momentum.mag2()+residualMass*residualMass));
    162    
    163        anInitialState.SetA(anA);
    164        anInitialState.SetZ(aZ);
    165        anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
    166        anInitialState.SetNumberOfCharged(numberOfCh);
    167        anInitialState.SetNumberOfHoles(numberOfHoles);
    168        anInitialState.SetMomentum(exciton4Momentum);
    169 //      anInitialState.SetExcitationEnergy(exEnergy); // now a redundant call.
    170 
    171      // call pre-compound
    172        const G4Fragment aFragment(anInitialState);
    173        G4ReactionProductVector * aPreResult = theDeExcitation->DeExcite(aFragment);
    174 
    175        // fill pre-compound part into the result, and return
    176        for(unsigned int ll=0; ll<aPreResult->size(); ll++)
    177        {
    178          theTotalResult->push_back(aPreResult->operator[](ll));
    179        }
    180        delete aPreResult;
    181      }
    182      else
    183      {
    184        // throw G4HadronicException(__FILE__, __LINE__, "Please register an evaporation phase with G4GeneratorPrecompoundInterface.");
    185      }
    186      // now return
    187      
    188      std::for_each(result->begin(), result->end(), DeleteKineticTrack());
    189      delete result;
    190      return theTotalResult;
    191    }
     184  std::for_each(result->begin(), result->end(), DeleteKineticTrack());
     185  delete result;
     186  return theTotalResult;
     187}
    192188 
    193 G4double G4GeneratorPrecompoundInterface::SetCaptureThreshold(G4double value)
     189G4HadFinalState* G4GeneratorPrecompoundInterface::
     190ApplyYourself(const G4HadProjectile &, G4Nucleus & )
    194191{
    195     G4double old=CaptureThreshold;
    196     CaptureThreshold=value;
    197     return old;
    198 
     192  G4cout << "G4GeneratorPrecompoundInterface: ApplyYourself interface called stand-allone."
     193         << G4endl;
     194  G4cout << "This class is only a mediator between generator and precompound"<<G4endl;
     195  G4cout << "Please remove from your physics list."<<G4endl;
     196  throw G4HadronicException(__FILE__, __LINE__, "SEVERE: G4GeneratorPrecompoundInterface model interface called stand-allone.");
     197  return new G4HadFinalState;
    199198}
Note: See TracChangeset for help on using the changeset viewer.