Ignore:
Timestamp:
Jun 18, 2010, 11:42:07 AM (14 years ago)
Author:
garnier
Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r962 r1315  
    4545#include "G4StableIsotopes.hh"
    4646#include "G4HadTmpUtil.hh"
    47 
    48 #include "G4HadLeadBias.hh"
     47#include "G4NucleiProperties.hh"
     48
     49//VI #include "G4HadLeadBias.hh"
    4950#include "G4HadronicException.hh"
    50 #include "G4HadReentrentException.hh"
    51 
    52 #include "G4HadronicWhiteBoard.hh"
    53 
    54 #include "G4HadSignalHandler.hh"
    55 
     51//VI #include "G4HadReentrentException.hh"
     52//VI #include "G4HadronicWhiteBoard.hh"
     53//VI #include "G4HadSignalHandler.hh"
    5654#include "G4HadronicProcessStore.hh"
    5755
    58 #include <typeinfo>
     56
     57/*VI
    5958
    6059namespace G4HadronicProcess_local
    … …  
    6564  }
    6665}
     66VI*/
    6767
    6868G4IsoParticleChange * G4HadronicProcess::theIsoResult = 0;
    … …  
    7878//////////////////////////////////////////////////////////////////
    7979
    80 G4HadronicProcess::G4HadronicProcess( const G4String &processName,
    81                                       G4ProcessType   aType ) :
    82 G4VDiscreteProcess( processName, aType)
    83 { 
     80G4HadronicProcess::G4HadronicProcess(const G4String& processName,
     81                                     G4ProcessType aType)
     82 :G4VDiscreteProcess(processName, aType)
     83{
    8484  ModelingState = 0;
    8585  isoIsOnAnyway = -1;
    8686  theTotalResult = new G4ParticleChange();
     87  theInteraction = 0;
    8788  theCrossSectionDataStore = new G4CrossSectionDataStore();
    8889  G4HadronicProcessStore::Instance()->Register(this);
    … …  
    9091  xBiasOn = false;
    9192  G4HadronicProcess_debug_flag = false;
    92   if(getenv("SwitchLeadBiasOn")) theBias.push_back(new G4HadLeadBias());
     93  epReportLevel = 0;
     94  epCheckLevels.first = DBL_MAX;
     95  epCheckLevels.second = DBL_MAX;
     96  levelsSetByProcess = false;
     97
     98  //VI   if(getenv("SwitchLeadBiasOn")) theBias.push_back(new G4HadLeadBias());
    9399}
    94100
    … …  
    100106  std::for_each(theProductionModels.begin(),
    101107                theProductionModels.end(), G4Delete());
    102   std::for_each(theBias.begin(), theBias.end(), G4Delete());
     108  //VI  std::for_each(theBias.begin(), theBias.end(), G4Delete());
    103109 
    104110  delete theOldIsoResult;
    … …  
    109115void G4HadronicProcess::RegisterMe( G4HadronicInteraction *a )
    110116{
     117  if(!a) { return; }
    111118  try{GetManagerPointer()->RegisterMe( a );}   
    112119  catch(G4HadronicException & aE)
    113120  {
    114     aE.Report(std::cout);
     121    G4cout << "Unrecoverable error in " << GetProcessName()
     122           << " to register " << a->GetModelName() << G4endl;
    115123    G4Exception("G4HadronicProcess", "007", FatalException,
    116124    "Could not register G4HadronicInteraction");
    … …  
    136144  try
    137145  {
    138     ModelingState = 1;   
     146    //VI ModelingState = 1;   
    139147    theLastCrossSection = aScaleFactor*
    140148      theCrossSectionDataStore->GetCrossSection(aTrack.GetDynamicParticle(),
    … …  
    143151  catch(G4HadronicException aR)
    144152  {
    145     aR.Report(G4cout);
     153    DumpState(aTrack,"GetMeanFreePath");
    146154    G4Exception("G4HadronicProcess", "007", FatalException,
    147155    "G4HadronicProcess::GetMeanFreePath failed");
    148156  }
    149157  G4double res = DBL_MAX;
    150   if( theLastCrossSection > 0.0 ) res = 1.0/theLastCrossSection;
     158  if( theLastCrossSection > 0.0 ) { res = 1.0/theLastCrossSection; }
    151159  return res;
    152160}
    … …  
    165173{
    166174  // Debugging stuff
    167   if(G4HadronicProcess_debug_flag)
    168          std::cout << "@@@@ hadronic process start "<< std::endl;
     175  //VI if(G4HadronicProcess_debug_flag)
     176  //VI       std::cout << "@@@@ hadronic process start "<< std::endl;
    169177  // G4cout << theNumberOfInteractionLengthLeft<<G4endl;
    170 #ifndef G4HadSignalHandler_off
    171   G4HadSignalHandler aHandler(G4HadronicProcess_local::G4HadronicProcessHandler_1);
    172 #endif
     178  //VI #ifndef G4HadSignalHandler_off
     179  //VI G4HadSignalHandler aHandler(G4HadronicProcess_local::G4HadronicProcessHandler_1);
     180  //VI #endif
    173181
    174182  if(aTrack.GetTrackStatus() != fAlive && aTrack.GetTrackStatus() != fSuspend) {
    … …  
    176184        aTrack.GetTrackStatus() == fKillTrackAndSecondaries ||
    177185        aTrack.GetTrackStatus() == fPostponeToNextEvent) {
    178       G4cerr << "G4HadronicProcess: track in unusable state - "
     186      G4cout << "G4HadronicProcess: track in unusable state - "
    179187             << aTrack.GetTrackStatus() << G4endl;
    180       G4cerr << "G4HadronicProcess: returning unchanged track " << G4endl;
     188      G4cout << "G4HadronicProcess: returning unchanged track " << G4endl;
     189      DumpState(aTrack,"PostStepDoIt");
    181190      G4Exception("G4HadronicProcess", "001", JustWarning, "bailing out");
    182191    }
    … …  
    224233  catch(G4HadronicException & aR)
    225234  {
    226     aR.Report(G4cout);
    227     G4cout << "Unrecoverable error for:"<<G4endl;
    228     G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
    229     G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
    230     G4cout << " - Particle type = "
    231     <<aParticle->GetDefinition()->GetParticleName()<<G4endl;
     235    DumpState(aTrack,"SampleZandA");
     236    //VI    G4cout << "Unrecoverable error for:"<<G4endl;
     237    //VI G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
     238    //VI G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
     239    //VI G4cout << " - Particle type = "
     240    //VI <<aParticle->GetDefinition()->GetParticleName()<<G4endl;
    232241    G4Exception("G4HadronicProcess", "007", FatalException,
    233242    "PostStepDoIt failed on element selection.");
    … …  
    236245  try
    237246  {
    238     theInteraction = ChooseHadronicInteraction( kineticEnergy,
    239     aMaterial, anElement );
     247    theInteraction =
     248      ChooseHadronicInteraction( kineticEnergy, aMaterial, anElement );
    240249  }
    241250  catch(G4HadronicException & aE)
    242251  {
    243     aE.Report(std::cout);
    244     G4cout << "Unrecoverable error for:"<<G4endl;
    245     G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
    246     G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
    247     G4cout << " - Particle type = "
    248            << aParticle->GetDefinition()->GetParticleName()<<G4endl;
     252    G4cout << "Target element "<<anElement->GetName()<<"  Z= "
     253           << targetNucleus.GetZ() << "  A= " << targetNucleus.GetN() << G4endl;
     254    DumpState(aTrack,"ChooseHadronicInteraction");
     255    //VI G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
     256    //VI G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
     257    //VI G4cout << " - Particle type = "
     258    //VI       << aParticle->GetDefinition()->GetParticleName()<<G4endl;
    249259    G4Exception("G4HadronicProcess", "007", FatalException,
    250260    "ChooseHadronicInteraction failed.");
    … …  
    264274      // Call the interaction
    265275      result = theInteraction->ApplyYourself( thePro, targetNucleus);
    266     }
     276      ++reentryCount;
     277    }
     278    /*VI
    267279    catch(G4HadReentrentException aR)
    268280    {
    … …  
    281293             << aParticle->GetDefinition()->GetParticleName() << G4endl;
    282294      result = 0; // here would still be leaking...
    283       if(reentryCount>100)
    284       {
    285         G4Exception("G4HadronicProcess", "007", FatalException,
    286         "GetHadronicProcess: Reentering ApplyYourself too often - PostStepDoIt failed."); 
    287       }
    288295      G4Exception("G4HadronicProcess", "007", FatalException,
    289296      "GetHadronicProcess: PostStepDoIt failed (Reentering ApplyYourself not yet supported.)"); 
    290297    }
     298    VI*/
    291299    catch(G4HadronicException aR)
    292300    {
    293       G4HadronicWhiteBoard & theBoard = G4HadronicWhiteBoard::Instance();
    294       theBoard.SetProjectile(thePro);
    295       theBoard.SetTargetNucleus(targetNucleus);
    296       theBoard.SetProcessName(GetProcessName());
    297       theBoard.SetModelName(theInteraction->GetModelName());
    298 
    299       aR.Report(G4cout);
    300       G4cout << " G4HadronicProcess failed in ApplyYourself call for"
    301              << G4endl;
    302       G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
    303       G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
    304       G4cout << " - Particle type = "
    305              << aParticle->GetDefinition()->GetParticleName() << G4endl;
     301      //VI G4HadronicWhiteBoard & theBoard = G4HadronicWhiteBoard::Instance();
     302      //VI theBoard.SetProjectile(thePro);
     303      //VI theBoard.SetTargetNucleus(targetNucleus);
     304      //VI theBoard.SetProcessName(GetProcessName());
     305      //VI theBoard.SetModelName(theInteraction->GetModelName());
     306
     307      G4cout << "Call for " << theInteraction->GetModelName() << G4endl;
     308      G4cout << "Target element "<<anElement->GetName()<<"  Z= "
     309             << targetNucleus.GetZ() << "  A= " << targetNucleus.GetN() << G4endl;
     310      DumpState(aTrack,"ApplyYourself");
     311      //VI G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
     312      //VI G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
     313      //VI G4cout << " - Particle type = "
     314      //VI        << aParticle->GetDefinition()->GetParticleName() << G4endl;
    306315      G4Exception("G4HadronicProcess", "007", FatalException,
    307316      "PostStepDoIt failed.");
    308317    }
     318    if(reentryCount>100) {
     319      G4cout << "Call for " << theInteraction->GetModelName() << G4endl;
     320      G4cout << "Target element "<<anElement->GetName()<<"  Z= "
     321             << targetNucleus.GetZ() << "  A= " << targetNucleus.GetN() << G4endl;
     322      DumpState(aTrack,"ApplyYourself");
     323      G4Exception("G4HadronicProcess", "007", FatalException,
     324                  "Reentering ApplyYourself too often - PostStepDoIt failed."); 
     325    }
    309326  }
    310327  while(!result);
    311328
    312   if(!ModelingState && !getenv("BypassAllSafetyChecks") )
    313   {
    314     G4cout << "ERROR IN EXECUTION -- HADRONIC PROCESS STATE NOT VALID"<<G4endl;
    315     G4cout << "Result will be of undefined quality."<<G4endl;
    316   }
     329  //VI if(!ModelingState && !getenv("BypassAllSafetyChecks") )
     330  //VI{
     331  //VI  G4cout << "ERROR IN EXECUTION -- HADRONIC PROCESS STATE NOT VALID"<<G4endl;
     332  //VI  G4cout << "Result will be of undefined quality."<<G4endl;
     333  //VI}
    317334
    318335  // NOT USED ?? Projectile particle has changed character during interaction
    319   if(result->GetStatusChange() == isAlive &&
    320      thePro.GetDefinition() != aTrack.GetDefinition())
    321   {
    322     G4DynamicParticle * aP =
    323              const_cast<G4DynamicParticle *>(aTrack.GetDynamicParticle());
    324     aP->SetDefinition(const_cast<G4ParticleDefinition *>(thePro.GetDefinition()));
    325   }
     336  //VI if(result->GetStatusChange() == isAlive &&
     337  //VI    thePro.GetDefinition() != aTrack.GetDefinition())
     338  //VI {
     339  //VI  G4DynamicParticle * aP =
     340  //VI           const_cast<G4DynamicParticle *>(aTrack.GetDynamicParticle());
     341  //VI   aP->SetDefinition(const_cast<G4ParticleDefinition *>(thePro.GetDefinition()));
     342  //VI }
    326343
    327344  result->SetTrafoToLab(thePro.GetTrafoToLab());
    328345
     346  /*VI
    329347  if(getenv("HadronicDoitLogging") )
    330348  {
    … …  
    338356           << G4endl;
    339357  }
     358  VI*/
    340359
    341360  ClearNumberOfInteractionLengthLeft();
    … …  
    348367  }
    349368 
     369  /*VI
    350370  G4double e=aTrack.GetKineticEnergy();
    351371  ModelingState = 0;
    … …  
    357377    }
    358378  }
    359 
     379  VI*/
    360380  // Put hadronic final state particles into G4ParticleChange
    361381
    362382  FillTotalResult(result, aTrack);
    363   if(G4HadronicProcess_debug_flag)
    364     std::cout << "@@@@ hadronic process end "<< std::endl;
    365    
     383  //VI if(G4HadronicProcess_debug_flag)
     384  //VI  std::cout << "@@@@ hadronic process end "<< std::endl;
     385
     386  if (epReportLevel > 0) CheckEnergyMomentumConservation(aTrack, targetNucleus);
     387 
    366388  return theTotalResult;
    367389}
    368390
     391#include <typeinfo>
    369392
    370393G4HadFinalState*
    … …  
    487510  theTotalResult->SetSecondaryWeightByProcess(true);
    488511  theTotalResult->ProposeTrackStatus(fAlive);
    489   G4double rotation = 2.*pi*G4UniformRand();
     512  G4double rotation = CLHEP::twopi*G4UniformRand();
    490513  G4ThreeVector it(0., 0., 1.);
    491514
    … …  
    520543      if(xBiasOn)
    521544      {
     545        DumpState(aT,"FillTotalResult");
    522546        G4Exception("G4HadronicProcess", "007", FatalException,
    523547        "Cannot cross-section bias a process that suspends tracks.");
    … …  
    575599  else
    576600  {
    577     G4cerr << "Track status is "<< aR->GetStatusChange()<<G4endl;
     601    G4cout << "Call for " << theInteraction->GetModelName() << G4endl;
     602    G4cout << "Target Z= "
     603           << targetNucleus.GetZ() << "  A= " << targetNucleus.GetN() << G4endl;
     604    DumpState(aT,"FillTotalResult");
    578605    G4Exception("G4HadronicProcess", "007", FatalException,
    579606    "use of unsupported track-status.");
    … …  
    691718      (it != "PositronNuclear") )
    692719    {
    693       G4Exception("G4HadronicProcess", "007", FatalException,
     720      G4Exception("G4HadronicProcess::BiasCrossSectionByFactor", "007", FatalException,
    694721                  "Cross-section biasing available only for gamma and electro nuclear reactions.");
    695722    }
    696723  if(aScale<100)
    697724    {
    698       G4Exception("G4HadronicProcess", "001", JustWarning,
     725      G4Exception("G4HadronicProcess::BiasCrossSectionByFactor", "001", JustWarning,
    699726                  "Cross-section bias readjusted to be above safe limit. New value is 100");
    700727      aScaleFactor = 100.;
    701728    }
    702729}
     730
     731void
     732G4HadronicProcess::CheckEnergyMomentumConservation(const G4Track& aTrack,
     733                                                   const G4Nucleus& aNucleus)
     734{
     735  G4double targetMass = G4NucleiProperties::GetNuclearMass(aNucleus.GetN(), aNucleus.GetZ());
     736  G4LorentzVector projectile4mom = aTrack.GetDynamicParticle()->Get4Momentum();
     737  G4LorentzVector target4mom(0, 0, 0, targetMass);
     738  G4LorentzVector initial4mom = projectile4mom + target4mom;
     739
     740  G4Track* sec;
     741  G4LorentzVector final4mom;
     742  G4int nSec = theTotalResult->GetNumberOfSecondaries();
     743  for (G4int i = 0; i < nSec; i++) {
     744    sec = theTotalResult->GetSecondary(i);
     745    final4mom += sec->GetDynamicParticle()->Get4Momentum();
     746  }
     747
     748  G4LorentzVector diff = initial4mom - final4mom;
     749  G4double absolute = diff.e();
     750  G4double relative = absolute/aTrack.GetKineticEnergy();
     751
     752  G4String processName = GetProcessName();
     753  G4HadronicInteraction* theModel = GetHadronicInteraction();
     754  G4String modelName("none");
     755  if (theModel) modelName = theModel->GetModelName();
     756
     757  std::pair<G4double, G4double> checkLevels = epCheckLevels;;
     758  if (!levelsSetByProcess) {
     759    if (theModel) checkLevels = theModel->GetEnergyMomentumCheckLevels();
     760  }
     761
     762  G4bool relPass = false;
     763  G4String relResult = "fail";
     764  if (std::abs(relative) < checkLevels.first) {
     765    relPass = true;
     766    relResult = "pass";
     767  }
     768
     769  G4bool absPass = false;
     770  G4String absResult = "fail";
     771  if (std::abs(absolute) < checkLevels.second) {
     772    absPass = true;
     773    absResult = "pass";
     774  }
     775
     776  // Options for level of reporting detail:
     777  //  0. off
     778  //  1. report only when E/p not conserved
     779  //  2. report regardless of E/p conservation
     780  //  3. report only when E/p not conserved, with model names, process names, and limits
     781  //  4. report regardless of E/p conservation, with model names, process names, and limits
     782
     783  if(epReportLevel == 4) {
     784    G4cout << " Process: " << processName << " , Model: " <<  modelName << G4endl;
     785    G4cout << " relative limit " << checkLevels.first << " relative value = "
     786           << relative << " " << relResult << G4endl;
     787    G4cout << " absolute limit " << checkLevels.second << " absolute value = "
     788           << absolute << " " << absResult << G4endl;
     789
     790  } else if(epReportLevel == 3) {
     791    if (!absPass || !relPass) {
     792      G4cout << " Process: " << processName << " , Model: " <<  modelName << G4endl;
     793      G4cout << " relative limit " << checkLevels.first << " relative value = "
     794             << relative << " " << relResult << G4endl; 
     795      G4cout << " absolute limit " << checkLevels.second << " absolute value = "
     796             << absolute << " " << absResult << G4endl;
     797    }
     798
     799  } else if(epReportLevel == 2) {
     800    G4cout << " relative value = " << relative << " " << relPass
     801           << " absolute value = " << absolute << " " << absPass << G4endl;
     802
     803  } else if(epReportLevel == 1) {
     804    if (!absPass || !relPass) {
     805      G4cout << " relative value = " << relative << " " << relPass
     806             << " absolute value = " << absolute << " " << absPass << G4endl;
     807    }
     808  }
     809}
     810
     811
     812void G4HadronicProcess::DumpState(const G4Track& aTrack, const G4String& method)
     813{
     814  G4cout << "Unrecoverable error in the method " << method << " of "
     815         << GetProcessName() << G4endl;
     816  G4cout << "TrackID= "<< aTrack.GetTrackID() << "  ParentID= " << aTrack.GetParentID()
     817         << "  " << aTrack.GetDefinition()->GetParticleName() << G4endl;
     818  G4cout << "Ekin(GeV)= " << aTrack.GetKineticEnergy()/CLHEP::GeV
     819         << ";  direction= " << aTrack.GetMomentumDirection() << G4endl;
     820  G4cout << "Position(mm)= " << aTrack.GetPosition()/CLHEP::mm
     821         << ";  material " << aTrack.GetMaterial()->GetName() << G4endl;
     822  G4cout << "PhysicalVolume  <" << aTrack.GetVolume()->GetName() << ">" << G4endl;
     823}
     824
    703825/* end of file */
Note: See TracChangeset for help on using the changeset viewer.