Ignore:
Timestamp:
Apr 6, 2009, 12:30:29 PM (17 years ago)
Author:
garnier
Message:

update processes

Location:
trunk/source/processes/hadronic/management/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/management/src/G4EnergyRangeManager.cc

    r819 r962  
    2626//
    2727// $Id: G4EnergyRangeManager.cc,v 1.15 2006/06/29 19:58:21 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-02-ref-02 $
    2929//
    3030 // Hadronic Process: Energy Range Manager

  • trunk/source/processes/hadronic/management/src/G4HadronInelasticProcess.cc

    r819 r962  
    2525//
    2626//
    27 //
    28  // Hadronic Inelastic Process Class
    29  // J.L. Chuma, TRIUMF, 24-Mar-1997
    30  // Last modified: 27-Mar-1997
    31  // J.P. Wellisch: Bug hunting, 23-Apr-97
    32  // Modified by J.L.Chuma 8-Jul-97 to eliminate possible division by zero for sigma
     27// Hadronic Inelastic Process Class
     28// J.L. Chuma, TRIUMF, 24-Mar-1997
     29// Last modified: 27-Mar-1997
     30// J.P. Wellisch: Bug hunting, 23-Apr-97
     31// Modified by J.L.Chuma 8-Jul-97 to eliminate possible division by zero for sigma
    3332//
    3433// 14-APR-98 F.W.Jones: variant G4HadronInelastic process for
    … …  
    3635//
    3736// 17-JUN-98 F.W.Jones: removed extraneous code causing core dump.
     37// 01-SEP-2008 V.Ivanchenko: use methods from the base class
    3838//
    3939 
    4040#include "G4HadronInelasticProcess.hh"
     41#include "G4HadronInelasticDataSet.hh"
    4142#include "G4GenericIon.hh"
    42 #include "G4ProcessManager.hh"
    43 #include "G4ProcessVector.hh"
    44 #include "G4HadronicException.hh"
     43#include "G4ParticleDefinition.hh"
    4544 
    46  void G4HadronInelasticProcess::BuildThePhysicsTable()
    47   {
    48     if (!G4HadronicProcess::GetCrossSectionDataStore()) {
    49       return;
    50     }
    51     G4HadronicProcess::GetCrossSectionDataStore()->BuildPhysicsTable(*theParticle);
    52   }
    53  
    54  G4HadronInelasticProcess::G4HadronInelasticProcess(
    55   const G4String &processName,
    56   G4ParticleDefinition *aParticle ) :
    57    G4HadronicProcess( processName )
    58  {
    59    G4HadronicProcess::AddDataSet(new G4HadronInelasticDataSet);
    60    theParticle = aParticle;
    61  }
     45G4HadronInelasticProcess::G4HadronInelasticProcess(const G4String& processName,
     46                                                   G4ParticleDefinition* aParticle):
     47  G4HadronicProcess(processName)
     48{
     49  SetProcessSubType(fHadronInelastic);
     50  AddDataSet(new G4HadronInelasticDataSet());
     51  theParticle = aParticle;
     52}
    6253
    63  G4HadronInelasticProcess::~G4HadronInelasticProcess() { }
     54G4HadronInelasticProcess::~G4HadronInelasticProcess()
     55{}
    6456
    65  G4VParticleChange *G4HadronInelasticProcess::
    66  PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
    67  {
    68    if(0==GetLastCrossSection()&&!getenv("DebugNeutronHP"))
    69    {
    70      G4cerr << "G4HadronInelasticProcess: called for final state, while cross-section was zero"<<G4endl;
    71      G4cerr << "                          Returning empty particle change...."<<G4endl;
    72      G4double dummy=0;
    73      G4ForceCondition condition;
    74      G4double it = GetMeanFreePath(aTrack, dummy, &condition);
    75      G4cerr << "                          current MeanFreePath is "<<it<<G4endl;
    76      theParticleChange.Initialize(aTrack);
    77      return &theParticleChange;
    78    }
    79    SetDispatch( this );
    80    return G4HadronicProcess::GeneralPostStepDoIt( aTrack, aStep );
    81  }
    82 
    83  G4bool G4HadronInelasticProcess::
    84  IsApplicable(const G4ParticleDefinition& aP)
    85  {
    86     return  theParticle == &aP || theParticle == G4GenericIon::GenericIon();
    87  }
    88 
    89  G4double G4HadronInelasticProcess::GetMicroscopicCrossSection(
    90   const G4DynamicParticle *aParticle,
    91   const G4Element *anElement,
    92   G4double aTemp)
    93   {
    94     // returns the microscopic cross section in GEANT4 internal units
    95    
    96    if (!G4HadronicProcess::GetCrossSectionDataStore())
    97    {
    98       throw G4HadronicException(__FILE__, __LINE__,
    99       "G4HadronInelasticProcess::GetMicroscopicCrossSection: "
    100                   "no CrossSectionDataStore");
    101       return DBL_MIN;
    102    }
    103    return G4HadronicProcess::GetCrossSectionDataStore()->GetCrossSection(aParticle, anElement, aTemp);
    104 
    105   }
    106  
    107  /* end of file */
     57G4bool G4HadronInelasticProcess::IsApplicable(const G4ParticleDefinition& aP)
     58{
     59  return  theParticle == &aP || theParticle == G4GenericIon::GenericIon();
     60}

  • trunk/source/processes/hadronic/management/src/G4HadronicProcess.cc

    r819 r962  
    2828#include "G4Types.hh"
    2929
    30 #include <fstream>
    31 #include <sstream>
    32 #include <stdlib.h>
     30//#include <fstream>
     31//#include <sstream>
     32//#include <stdlib.h>
    3333#include "G4HadronicProcess.hh"
    34 // #include "G4EffectiveCharge.hh"
     34
    3535#include "G4HadProjectile.hh"
    3636#include "G4ElementVector.hh"
    … …  
    4949#include "G4HadronicException.hh"
    5050#include "G4HadReentrentException.hh"
    51 #include "G4HadronicInteractionWrapper.hh"
     51
     52#include "G4HadronicWhiteBoard.hh"
    5253
    5354#include "G4HadSignalHandler.hh"
     55
     56#include "G4HadronicProcessStore.hh"
    5457
    5558#include <typeinfo>
    … …  
    7275void G4HadronicProcess::
    7376DisableIsotopeProductionGlobally() {isoIsEnabled = false;}
     77
     78//////////////////////////////////////////////////////////////////
    7479
    7580G4HadronicProcess::G4HadronicProcess( const G4String &processName,
    … …  
    8186  theTotalResult = new G4ParticleChange();
    8287  theCrossSectionDataStore = new G4CrossSectionDataStore();
     88  G4HadronicProcessStore::Instance()->Register(this);
    8389  aScaleFactor = 1;
    8490  xBiasOn = false;
     91  G4HadronicProcess_debug_flag = false;
    8592  if(getenv("SwitchLeadBiasOn")) theBias.push_back(new G4HadLeadBias());
    8693}
    … …  
    8895G4HadronicProcess::~G4HadronicProcess()
    8996{
     97  G4HadronicProcessStore::Instance()->DeRegister(this);
    9098  delete theTotalResult;
    9199
    … …  
    94102  std::for_each(theBias.begin(), theBias.end(), G4Delete());
    95103 
    96   delete theOldIsoResult; delete theIsoResult;
     104  delete theOldIsoResult;
     105  delete theIsoResult;
    97106  delete theCrossSectionDataStore;
    98107}
    … …  
    107116    "Could not register G4HadronicInteraction");
    108117  }
     118  G4HadronicProcessStore::Instance()->RegisterInteraction(this, a); 
     119}
     120
     121void G4HadronicProcess::PreparePhysicsTable(const G4ParticleDefinition& p)
     122{
     123  if(getenv("G4HadronicProcess_debug")) G4HadronicProcess_debug_flag = true;
     124  G4HadronicProcessStore::Instance()->RegisterParticle(this, &p);
     125}
     126
     127void G4HadronicProcess::BuildPhysicsTable(const G4ParticleDefinition& p)
     128{
     129  theCrossSectionDataStore->BuildPhysicsTable(p);
     130  G4HadronicProcessStore::Instance()->PrintInfo(&p);
    109131}
    110132
    … …  
    112134GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *)
    113135{
    114   G4double sigma = 0.0;
    115136  try
    116137  {
    117     const G4DynamicParticle *aParticle = aTrack.GetDynamicParticle();
    118     if( !IsApplicable(*aParticle->GetDefinition()))
    119     {
    120       G4cout << "Unrecoverable error: "<<G4endl;
    121       G4ProcessManager * it = aParticle->GetDefinition()->GetProcessManager();
    122       G4ProcessVector * itv = it->GetProcessList();
    123       G4cout <<aParticle->GetDefinition()->GetParticleName()<<
    124       " has the following processes:"<<G4endl;
    125       for(G4int i=0; i<itv->size(); i++)
    126       {
    127         G4cout <<"  "<<(*itv)[i]->GetProcessName()<<G4endl;             
    128       }
    129       G4cout << "for kinetic energy "<<aParticle->GetKineticEnergy()<<G4endl;
    130       G4cout << "and material "<<aTrack.GetMaterial()->GetName()<<G4endl;
    131       G4Exception("G4HadronicProcess", "007", FatalException,
    132       std::string(this->GetProcessName()+
    133       " was called for "+
    134       aParticle->GetDefinition()->GetParticleName()).c_str() );
    135     }
    136     G4Material *aMaterial = aTrack.GetMaterial();
    137     ModelingState = 1;
    138    
    139     sigma = theCrossSectionDataStore->GetCrossSection(aParticle, aMaterial);
    140 
    141     sigma *= aScaleFactor;
    142     theLastCrossSection = sigma;
     138    ModelingState = 1;   
     139    theLastCrossSection = aScaleFactor*
     140      theCrossSectionDataStore->GetCrossSection(aTrack.GetDynamicParticle(),
     141                                                aTrack.GetMaterial());
    143142  }
    144143  catch(G4HadronicException aR)
    … …  
    148147    "G4HadronicProcess::GetMeanFreePath failed");
    149148  }
    150   if( sigma > 0.0 )
    151     return 1.0/sigma;
    152   else
    153     return DBL_MAX;
    154 }
    155 
    156 
    157 G4Element* G4HadronicProcess::ChooseAandZ(
    158 const G4DynamicParticle *aParticle, const G4Material *aMaterial )
    159 {
    160   std::pair<G4double, G4double> ZA =
    161       theCrossSectionDataStore->SelectRandomIsotope(aParticle, aMaterial);
    162   G4double ZZ = ZA.first;
    163   G4double AA = ZA.second;
    164 
    165   targetNucleus.SetParameters(AA, ZZ);
    166 
    167   const G4int numberOfElements = aMaterial->GetNumberOfElements();
    168   const G4ElementVector* theElementVector = aMaterial->GetElementVector();
    169   G4Element* chosen = 0;
    170   for (G4int i = 0; i < numberOfElements; i++) {
    171     chosen = (*theElementVector)[i];
    172     if (chosen->GetZ() == ZZ) break;
    173   }
    174   return chosen;
    175 }
    176 
    177 
    178 struct G4Nancheck{ bool operator()(G4double aV){return (!(aV<1))&&(!(aV>-1));}};
    179 
    180 G4VParticleChange *G4HadronicProcess::GeneralPostStepDoIt(
    181 const G4Track &aTrack, const G4Step &)
     149  G4double res = DBL_MAX;
     150  if( theLastCrossSection > 0.0 ) res = 1.0/theLastCrossSection;
     151  return res;
     152}
     153
     154G4double G4HadronicProcess::
     155GetMicroscopicCrossSection(const G4DynamicParticle *aParticle,
     156                           const G4Element *anElement,
     157                           G4double aTemp )
     158{
     159  return
     160    theCrossSectionDataStore->GetCrossSection(aParticle, anElement, aTemp);
     161}
     162
     163G4VParticleChange *G4HadronicProcess::PostStepDoIt(
     164                   const G4Track &aTrack, const G4Step &)
    182165{
    183166  // Debugging stuff
    184 
    185   bool G4HadronicProcess_debug_flag = false;
    186   if(getenv("G4HadronicProcess_debug")) G4HadronicProcess_debug_flag = true;
    187167  if(G4HadronicProcess_debug_flag)
    188168         std::cout << "@@@@ hadronic process start "<< std::endl;
    189169  // G4cout << theNumberOfInteractionLengthLeft<<G4endl;
    190   #ifndef G4HadSignalHandler_off
     170#ifndef G4HadSignalHandler_off
    191171  G4HadSignalHandler aHandler(G4HadronicProcess_local::G4HadronicProcessHandler_1);
    192   #endif
    193 
    194   if(aTrack.GetTrackStatus() != fAlive && aTrack.GetTrackStatus() != fSuspend)
    195   {
    196     G4cerr << "G4HadronicProcess: track in unusable state - "
    197     <<aTrack.GetTrackStatus()<<G4endl;
    198     G4cerr << "G4HadronicProcess: returning unchanged track "<<G4endl;
    199     G4Exception("G4HadronicProcess", "001", JustWarning, "bailing out");
     172#endif
     173
     174  if(aTrack.GetTrackStatus() != fAlive && aTrack.GetTrackStatus() != fSuspend) {
     175    if (aTrack.GetTrackStatus() == fStopAndKill ||
     176        aTrack.GetTrackStatus() == fKillTrackAndSecondaries ||
     177        aTrack.GetTrackStatus() == fPostponeToNextEvent) {
     178      G4cerr << "G4HadronicProcess: track in unusable state - "
     179             << aTrack.GetTrackStatus() << G4endl;
     180      G4cerr << "G4HadronicProcess: returning unchanged track " << G4endl;
     181      G4Exception("G4HadronicProcess", "001", JustWarning, "bailing out");
     182    }
     183    // No warning for fStopButAlive which is a legal status here
    200184    theTotalResult->Clear();
    201185    theTotalResult->Initialize(aTrack);
    … …  
    203187  }
    204188
    205   const G4DynamicParticle *aParticle = aTrack.GetDynamicParticle();
    206   G4Material *aMaterial = aTrack.GetMaterial();
     189  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
     190  G4Material* aMaterial = aTrack.GetMaterial();
    207191  G4double originalEnergy = aParticle->GetKineticEnergy();
    208192  G4double kineticEnergy = originalEnergy;
    209193
    210   // More debugging
    211 
     194  /*
     195  // It is not needed with standard NaN check
     196  // More debugging 
    212197  G4Nancheck go_wild;
    213198  if(go_wild(originalEnergy) ||
    … …  
    223208    return theTotalResult;
    224209  }
     210  */
    225211
    226212  // Get kinetic energy per nucleon for ions
    227 
    228213  if(aParticle->GetDefinition()->GetBaryonNumber() > 1.5)
    229214          kineticEnergy/=aParticle->GetDefinition()->GetBaryonNumber();
    … …  
    232217  try
    233218  {
    234     anElement = ChooseAandZ( aParticle, aMaterial );
     219    //    anElement = ChooseAandZ( aParticle, aMaterial );
     220    anElement = theCrossSectionDataStore->SampleZandA(aParticle,
     221                                                      aMaterial,
     222                                                      targetNucleus);
    235223  }
    236224  catch(G4HadronicException & aR)
    … …  
    243231    <<aParticle->GetDefinition()->GetParticleName()<<G4endl;
    244232    G4Exception("G4HadronicProcess", "007", FatalException,
    245     "GeneralPostStepDoIt failed on element selection.");
     233    "PostStepDoIt failed on element selection.");
    246234  }
    247235
    … …  
    275263    {
    276264      // Call the interaction
    277 
    278       G4HadronicInteractionWrapper aW;
    279       result = aW.ApplyInteraction(thePro, targetNucleus, theInteraction,
    280                                    GetProcessName(),
    281                                    theInteraction->GetModelName());
     265      result = theInteraction->ApplyYourself( thePro, targetNucleus);
    282266    }
    283267    catch(G4HadReentrentException aR)
    284268    {
     269      G4HadronicWhiteBoard & theBoard = G4HadronicWhiteBoard::Instance();
     270      theBoard.SetProjectile(thePro);
     271      theBoard.SetTargetNucleus(targetNucleus);
     272      theBoard.SetProcessName(GetProcessName());
     273      theBoard.SetModelName(theInteraction->GetModelName());
     274
    285275      aR.Report(G4cout);
    286276      G4cout << " G4HadronicProcess re-entering the ApplyYourself call for "
    … …  
    294284      {
    295285        G4Exception("G4HadronicProcess", "007", FatalException,
    296         "GetHadronicProcess: Reentering ApplyYourself too often - GeneralPostStepDoIt failed."); 
     286        "GetHadronicProcess: Reentering ApplyYourself too often - PostStepDoIt failed."); 
    297287      }
    298288      G4Exception("G4HadronicProcess", "007", FatalException,
    299       "GetHadronicProcess: GeneralPostStepDoIt failed (Reentering ApplyYourself not yet supported.)"); 
     289      "GetHadronicProcess: PostStepDoIt failed (Reentering ApplyYourself not yet supported.)"); 
    300290    }
    301291    catch(G4HadronicException aR)
    302292    {
     293      G4HadronicWhiteBoard & theBoard = G4HadronicWhiteBoard::Instance();
     294      theBoard.SetProjectile(thePro);
     295      theBoard.SetTargetNucleus(targetNucleus);
     296      theBoard.SetProcessName(GetProcessName());
     297      theBoard.SetModelName(theInteraction->GetModelName());
     298
    303299      aR.Report(G4cout);
    304300      G4cout << " G4HadronicProcess failed in ApplyYourself call for"
    … …  
    309305             << aParticle->GetDefinition()->GetParticleName() << G4endl;
    310306      G4Exception("G4HadronicProcess", "007", FatalException,
    311       "GeneralPostStepDoIt failed.");
     307      "PostStepDoIt failed.");
    312308    }
    313309  }
    … …  
    331327  result->SetTrafoToLab(thePro.GetTrafoToLab());
    332328
    333   /*
    334   // Loop over charged ion secondaries
    335 
    336   for(G4int i=0; i<result->GetNumberOfSecondaries(); i++)
    337   {
    338     G4DynamicParticle* aSecTrack = result->GetSecondary(i)->GetParticle();
    339     if(aSecTrack->GetDefinition()->GetPDGCharge()>1.5)
    340     {
    341       G4EffectiveCharge aCalculator;
    342       G4double charge =
    343            aCalculator.GetCharge(aMaterial, aSecTrack->GetKineticEnergy(),
    344                                  aSecTrack->GetDefinition()->GetPDGMass(),
    345                                  aSecTrack->GetDefinition()->GetPDGCharge());
    346       if(getenv("GHADChargeDebug"))
    347       {
    348         std::cout << "Recoil fractional charge is "
    349         << charge/aSecTrack->GetDefinition()->GetPDGCharge()<<" "
    350         << charge <<" "<<aSecTrack->GetDefinition()->GetPDGCharge()<<std::endl;
    351       }
    352       aSecTrack->SetCharge(charge);
    353     }
    354   }
    355   */
    356 
    357329  if(getenv("HadronicDoitLogging") )
    358330  {
    359331    G4cout << "HadronicDoitLogging "
    360     << GetProcessName() <<" "
    361     << aParticle->GetDefinition()->GetPDGEncoding()<<" "
    362     << originalEnergy<<" "
    363     << aParticle->GetMomentum()<<" "
    364     << targetNucleus.GetN()<<" "
    365     << targetNucleus.GetZ()<<" "
    366     << G4endl;
     332           << GetProcessName() <<" "
     333           << aParticle->GetDefinition()->GetPDGEncoding()<<" "
     334           << originalEnergy<<" "
     335           << aParticle->GetMomentum()<<" "
     336           << targetNucleus.GetN()<<" "
     337           << targetNucleus.GetZ()<<" "
     338           << G4endl;
    367339  }
    368340
    … …  
    509481G4HadronicProcess::FillTotalResult(G4HadFinalState * aR, const G4Track & aT)
    510482{
    511   G4Nancheck go_wild;
     483  //  G4Nancheck go_wild;
    512484  theTotalResult->Clear();
    513485  theTotalResult->ProposeLocalEnergyDeposit(0.);
    … …  
    559531      theTotalResult->ProposeParentWeight( XBiasSurvivalProbability()*aT.GetWeight() );
    560532      G4double newWeight = aR->GetWeightChange()*aT.GetWeight();
     533      /*
    561534      if(go_wild(aR->GetEnergyChange()))
    562535      {
    … …  
    571544        "surviving track received NaN momentum."); 
    572545      }
     546      */
    573547      G4double newM=aT.GetDefinition()->GetPDGMass();
    574548      G4double newE=aR->GetEnergyChange() + newM;
    … …  
    585559      if(aR->GetEnergyChange()>-.5)
    586560      {
     561        /*
    587562        if(go_wild(aR->GetEnergyChange()))
    588563        {
    … …  
    590565          "track received NaN energy."); 
    591566        }
     567        */
    592568        theTotalResult->ProposeEnergy(aR->GetEnergyChange());
    593569      }
    … …  
    606582  if(GetProcessName() != "hElastic" && GetProcessName() != "HadronElastic"
    607583     &&  theTotalResult->GetTrackStatus()==fAlive
    608      && aR->GetStatusChange()==isAlive
    609     )
    610   {
     584     && aR->GetStatusChange()==isAlive)
     585    {
    611586    // Use for debugging:   G4double newWeight = theTotalResult->GetParentWeight();
    612587
    … …  
    643618    theM.rotate(rotation, it);
    644619    theM*=aR->GetTrafoToLab();
    645 
     620    /*
    646621    if(go_wild(theM.e()))
    647622    {
    … …  
    656631      "secondary track received NaN momentum."); 
    657632    }
    658 
     633    */
    659634    aR->GetSecondary(i)->GetParticle()->Set4Momentum(theM);
    660635    G4double time = aR->GetSecondary(i)->GetTime();
    … …  
    662637
    663638    G4Track* track = new G4Track(aR->GetSecondary(i)->GetParticle(),
    664     time,
    665     aT.GetPosition());
     639                                time,
     640                                aT.GetPosition());
    666641
    667642    G4double newWeight = aT.GetWeight()*aR->GetSecondary(i)->GetWeight();
    … …  
    679654    // <<G4endl;
    680655    track->SetWeight(newWeight);
     656    /*
    681657    G4double trackDeb = track->GetKineticEnergy();
    682658    if( (  trackDeb<0
    … …  
    689665        <<G4endl;
    690666      }
    691       track->SetTouchableHandle(aT.GetTouchableHandle());
    692       theTotalResult->AddSecondary(track);
     667    */
     668    track->SetTouchableHandle(aT.GetTouchableHandle());
     669    theTotalResult->AddSecondary(track);
    693670  }
    694671
    … …  
    696673  return;
    697674}
     675
     676G4IsoParticleChange* G4HadronicProcess::GetIsotopeProductionInfo()
     677{
     678  G4IsoParticleChange * anIsoResult = theIsoResult;
     679  if(theIsoResult) theOldIsoResult = theIsoResult;
     680  theIsoResult = 0;
     681  return anIsoResult;
     682}
     683
     684void G4HadronicProcess::BiasCrossSectionByFactor(G4double aScale)
     685{
     686  xBiasOn = true;
     687  aScaleFactor = aScale;
     688  G4String it = GetProcessName();
     689  if( (it != "PhotonInelastic") &&
     690      (it != "ElectroNuclear") &&
     691      (it != "PositronNuclear") )
     692    {
     693      G4Exception("G4HadronicProcess", "007", FatalException,
     694                  "Cross-section biasing available only for gamma and electro nuclear reactions.");
     695    }
     696  if(aScale<100)
     697    {
     698      G4Exception("G4HadronicProcess", "001", JustWarning,
     699                  "Cross-section bias readjusted to be above safe limit. New value is 100");
     700      aScaleFactor = 100.;
     701    }
     702}
    698703/* end of file */
Note: See TracChangeset for help on using the changeset viewer.