source: trunk/source/processes/hadronic/management/src/G4HadronicProcess.cc @ 1315

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 27.9 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26//
27
28#include "G4Types.hh"
29
30//#include <fstream>
31//#include <sstream>
32//#include <stdlib.h>
33#include "G4HadronicProcess.hh"
34
35#include "G4HadProjectile.hh"
36#include "G4ElementVector.hh"
37#include "G4Track.hh"
38#include "G4Step.hh"
39#include "G4Element.hh"
40#include "G4ParticleChange.hh"
41#include "G4TransportationManager.hh"
42#include "G4Navigator.hh"
43#include "G4ProcessVector.hh"
44#include "G4ProcessManager.hh"
45#include "G4StableIsotopes.hh"
46#include "G4HadTmpUtil.hh"
47#include "G4NucleiProperties.hh"
48
49//VI #include "G4HadLeadBias.hh"
50#include "G4HadronicException.hh"
51//VI #include "G4HadReentrentException.hh"
52//VI #include "G4HadronicWhiteBoard.hh"
53//VI #include "G4HadSignalHandler.hh"
54#include "G4HadronicProcessStore.hh"
55
56
57/*VI
58
59namespace G4HadronicProcess_local
60{
61  extern "C" void G4HadronicProcessHandler_1(int)
62  {
63    G4HadronicWhiteBoard::Instance().Dump();
64  }
65}
66VI*/
67
68G4IsoParticleChange * G4HadronicProcess::theIsoResult = 0;
69G4IsoParticleChange * G4HadronicProcess::theOldIsoResult = 0;
70G4bool G4HadronicProcess::isoIsEnabled = true;
71
72void G4HadronicProcess::
73EnableIsotopeProductionGlobally()  {isoIsEnabled = true;}
74
75void G4HadronicProcess::
76DisableIsotopeProductionGlobally() {isoIsEnabled = false;}
77
78//////////////////////////////////////////////////////////////////
79
80G4HadronicProcess::G4HadronicProcess(const G4String& processName,
81                                     G4ProcessType aType)
82 :G4VDiscreteProcess(processName, aType)
83{
84  ModelingState = 0;
85  isoIsOnAnyway = -1;
86  theTotalResult = new G4ParticleChange();
87  theInteraction = 0;
88  theCrossSectionDataStore = new G4CrossSectionDataStore();
89  G4HadronicProcessStore::Instance()->Register(this);
90  aScaleFactor = 1;
91  xBiasOn = false;
92  G4HadronicProcess_debug_flag = false;
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());
99}
100
101G4HadronicProcess::~G4HadronicProcess()
102{ 
103  G4HadronicProcessStore::Instance()->DeRegister(this);
104  delete theTotalResult;
105
106  std::for_each(theProductionModels.begin(),
107                theProductionModels.end(), G4Delete());
108  //VI  std::for_each(theBias.begin(), theBias.end(), G4Delete());
109 
110  delete theOldIsoResult; 
111  delete theIsoResult;
112  delete theCrossSectionDataStore;
113}
114
115void G4HadronicProcess::RegisterMe( G4HadronicInteraction *a )
116{ 
117  if(!a) { return; }
118  try{GetManagerPointer()->RegisterMe( a );}   
119  catch(G4HadronicException & aE)
120  {
121    G4cout << "Unrecoverable error in " << GetProcessName() 
122           << " to register " << a->GetModelName() << G4endl;
123    G4Exception("G4HadronicProcess", "007", FatalException,
124    "Could not register G4HadronicInteraction");
125  }
126  G4HadronicProcessStore::Instance()->RegisterInteraction(this, a); 
127}
128
129void G4HadronicProcess::PreparePhysicsTable(const G4ParticleDefinition& p)
130{
131  if(getenv("G4HadronicProcess_debug")) G4HadronicProcess_debug_flag = true;
132  G4HadronicProcessStore::Instance()->RegisterParticle(this, &p);
133}
134
135void G4HadronicProcess::BuildPhysicsTable(const G4ParticleDefinition& p)
136{
137  theCrossSectionDataStore->BuildPhysicsTable(p);
138  G4HadronicProcessStore::Instance()->PrintInfo(&p);
139}
140
141G4double G4HadronicProcess::
142GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *)
143{ 
144  try
145  {
146    //VI ModelingState = 1;   
147    theLastCrossSection = aScaleFactor* 
148      theCrossSectionDataStore->GetCrossSection(aTrack.GetDynamicParticle(), 
149                                                aTrack.GetMaterial());
150  }
151  catch(G4HadronicException aR)
152  { 
153    DumpState(aTrack,"GetMeanFreePath");
154    G4Exception("G4HadronicProcess", "007", FatalException,
155    "G4HadronicProcess::GetMeanFreePath failed");
156  } 
157  G4double res = DBL_MAX;
158  if( theLastCrossSection > 0.0 ) { res = 1.0/theLastCrossSection; }
159  return res;
160}
161
162G4double G4HadronicProcess::
163GetMicroscopicCrossSection(const G4DynamicParticle *aParticle, 
164                           const G4Element *anElement, 
165                           G4double aTemp )
166{
167  return
168    theCrossSectionDataStore->GetCrossSection(aParticle, anElement, aTemp);
169}
170
171G4VParticleChange *G4HadronicProcess::PostStepDoIt(
172                   const G4Track &aTrack, const G4Step &)
173{
174  // Debugging stuff
175  //VI if(G4HadronicProcess_debug_flag)
176  //VI       std::cout << "@@@@ hadronic process start "<< std::endl;
177  // G4cout << theNumberOfInteractionLengthLeft<<G4endl;
178  //VI #ifndef G4HadSignalHandler_off
179  //VI G4HadSignalHandler aHandler(G4HadronicProcess_local::G4HadronicProcessHandler_1);
180  //VI #endif
181
182  if(aTrack.GetTrackStatus() != fAlive && aTrack.GetTrackStatus() != fSuspend) {
183    if (aTrack.GetTrackStatus() == fStopAndKill ||
184        aTrack.GetTrackStatus() == fKillTrackAndSecondaries ||
185        aTrack.GetTrackStatus() == fPostponeToNextEvent) {
186      G4cout << "G4HadronicProcess: track in unusable state - "
187             << aTrack.GetTrackStatus() << G4endl;
188      G4cout << "G4HadronicProcess: returning unchanged track " << G4endl;
189      DumpState(aTrack,"PostStepDoIt");
190      G4Exception("G4HadronicProcess", "001", JustWarning, "bailing out");
191    }
192    // No warning for fStopButAlive which is a legal status here
193    theTotalResult->Clear();
194    theTotalResult->Initialize(aTrack);
195    return theTotalResult;
196  }
197
198  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
199  G4Material* aMaterial = aTrack.GetMaterial();
200  G4double originalEnergy = aParticle->GetKineticEnergy();
201  G4double kineticEnergy = originalEnergy;
202
203  /*
204  // It is not needed with standard NaN check
205  // More debugging 
206  G4Nancheck go_wild;
207  if(go_wild(originalEnergy) ||
208    go_wild(aParticle->Get4Momentum().x()) ||
209  go_wild(aParticle->Get4Momentum().y()) ||
210  go_wild(aParticle->Get4Momentum().z()) ||
211  go_wild(aParticle->Get4Momentum().t())
212  )
213  {
214    G4Exception("G4HadronicProcess", "001", JustWarning, "NaN in input energy or momentum - bailing out.");
215    theTotalResult->Clear();
216    theTotalResult->Initialize(aTrack);
217    return theTotalResult;
218  }
219  */
220
221  // Get kinetic energy per nucleon for ions
222  if(aParticle->GetDefinition()->GetBaryonNumber() > 1.5)
223          kineticEnergy/=aParticle->GetDefinition()->GetBaryonNumber();
224
225  G4Element* anElement = 0;
226  try
227  {
228    //    anElement = ChooseAandZ( aParticle, aMaterial );
229    anElement = theCrossSectionDataStore->SampleZandA(aParticle, 
230                                                      aMaterial, 
231                                                      targetNucleus);
232  }
233  catch(G4HadronicException & aR)
234  {
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;
241    G4Exception("G4HadronicProcess", "007", FatalException,
242    "PostStepDoIt failed on element selection.");
243  }
244
245  try
246  {
247    theInteraction = 
248      ChooseHadronicInteraction( kineticEnergy, aMaterial, anElement );
249  }
250  catch(G4HadronicException & aE)
251  {
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;
259    G4Exception("G4HadronicProcess", "007", FatalException,
260    "ChooseHadronicInteraction failed.");
261  }
262
263  // Initialize the hadronic projectile from the track
264
265  G4HadProjectile thePro(aTrack);
266 
267  G4HadFinalState* result = 0;
268  G4int reentryCount = 0;
269
270  do
271  {
272    try
273    {
274      // Call the interaction
275      result = theInteraction->ApplyYourself( thePro, targetNucleus);
276      ++reentryCount;
277    }
278    /*VI
279    catch(G4HadReentrentException aR)
280    {
281      G4HadronicWhiteBoard & theBoard = G4HadronicWhiteBoard::Instance();
282      theBoard.SetProjectile(thePro);
283      theBoard.SetTargetNucleus(targetNucleus);
284      theBoard.SetProcessName(GetProcessName());
285      theBoard.SetModelName(theInteraction->GetModelName());
286
287      aR.Report(G4cout);
288      G4cout << " G4HadronicProcess re-entering the ApplyYourself call for "
289             <<G4endl;
290      G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
291      G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
292      G4cout << " - Particle type = "
293             << aParticle->GetDefinition()->GetParticleName() << G4endl;
294      result = 0; // here would still be leaking...
295      G4Exception("G4HadronicProcess", "007", FatalException,
296      "GetHadronicProcess: PostStepDoIt failed (Reentering ApplyYourself not yet supported.)"); 
297    }
298    VI*/
299    catch(G4HadronicException aR)
300    {
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;
315      G4Exception("G4HadronicProcess", "007", FatalException,
316      "PostStepDoIt failed.");
317    }
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    }
326  }
327  while(!result);
328
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}
334
335  // NOT USED ?? Projectile particle has changed character during interaction
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 }
343
344  result->SetTrafoToLab(thePro.GetTrafoToLab());
345
346  /*VI
347  if(getenv("HadronicDoitLogging") )
348  {
349    G4cout << "HadronicDoitLogging "
350           << GetProcessName() <<" "
351           << aParticle->GetDefinition()->GetPDGEncoding()<<" "
352           << originalEnergy<<" "
353           << aParticle->GetMomentum()<<" "
354           << targetNucleus.GetN()<<" "
355           << targetNucleus.GetZ()<<" "
356           << G4endl;
357  }
358  VI*/
359
360  ClearNumberOfInteractionLengthLeft();
361  if(isoIsOnAnyway!=-1)
362  {
363    if(isoIsEnabled||isoIsOnAnyway)
364    {
365      result = DoIsotopeCounting(result, aTrack, targetNucleus);
366    }
367  }
368 
369  /*VI
370  G4double e=aTrack.GetKineticEnergy();
371  ModelingState = 0;
372  if(e<5*GeV)
373  {
374    for(size_t i=0; i<theBias.size(); i++)
375    {
376      result = theBias[i]->Bias(result);
377    }
378  }
379  VI*/
380  // Put hadronic final state particles into G4ParticleChange
381
382  FillTotalResult(result, aTrack);
383  //VI if(G4HadronicProcess_debug_flag)
384  //VI  std::cout << "@@@@ hadronic process end "<< std::endl;
385
386  if (epReportLevel > 0) CheckEnergyMomentumConservation(aTrack, targetNucleus);
387 
388  return theTotalResult;
389}
390
391#include <typeinfo>
392
393G4HadFinalState* 
394G4HadronicProcess::DoIsotopeCounting(G4HadFinalState * aResult,
395                                     const G4Track & aTrack,
396                                     const G4Nucleus & aNucleus)
397{
398  // get the PC from iso-production
399  delete theOldIsoResult;
400  theOldIsoResult = 0;
401  delete theIsoResult;
402  theIsoResult = new G4IsoParticleChange;
403  G4bool done = false;
404  G4IsoResult * anIsoResult = 0;
405  for(unsigned int i=0; i<theProductionModels.size(); i++)
406  {
407    anIsoResult = theProductionModels[i]->GetIsotope(aTrack, aNucleus);
408    if(anIsoResult!=0)
409    {
410      done = true;
411      break;
412    }
413  }
414
415  // If no production models active, use default iso production
416  if(!done) anIsoResult = ExtractResidualNucleus(aTrack, aNucleus, aResult); 
417
418  // Add all info explicitely and add typename from model called.
419  theIsoResult->SetIsotope(anIsoResult->GetIsotope());
420  theIsoResult->SetProductionPosition(aTrack.GetPosition());
421  theIsoResult->SetProductionTime(aTrack.GetGlobalTime());
422  theIsoResult->SetParentParticle(*aTrack.GetDynamicParticle());
423  theIsoResult->SetMotherNucleus(anIsoResult->GetMotherNucleus());
424  theIsoResult->SetProducer(typeid(*theInteraction).name());
425 
426  delete anIsoResult;
427
428  // If isotope production is enabled the GetIsotopeProductionInfo()
429  // method must be called or else a memory leak will result
430  //
431  // The following code will fix the memory leak, but remove the
432  // isotope information:
433  //
434  //  if(theIsoResult) {
435  //    delete theIsoResult;
436  //    theIsoResult = 0;
437  //  }
438 
439  return aResult;
440}
441
442G4IsoResult* 
443G4HadronicProcess::ExtractResidualNucleus(const G4Track&,
444                                          const G4Nucleus& aNucleus,
445                                          G4HadFinalState* aResult)
446{
447  G4double A = aNucleus.GetN();
448  G4double Z = aNucleus.GetZ();
449  G4double bufferA = 0;
450  G4double bufferZ = 0;
451 
452  // loop over aResult, and decrement A, Z accordingly
453  // cash the max
454  for(G4int i=0; i<aResult->GetNumberOfSecondaries(); i++)
455  {
456    G4HadSecondary* aSecTrack = aResult->GetSecondary(i);
457    if(bufferA<aSecTrack->GetParticle()->GetDefinition()->GetBaryonNumber())
458    {
459      bufferA = aSecTrack->GetParticle()->GetDefinition()->GetBaryonNumber();
460      bufferZ = aSecTrack->GetParticle()->GetDefinition()->GetPDGCharge();
461    }
462    Z-=aSecTrack->GetParticle()->GetDefinition()->GetPDGCharge();
463    A-=aSecTrack->GetParticle()->GetDefinition()->GetBaryonNumber();
464  }
465 
466  // if the fragment was part of the final state, it is
467  // assumed to be the heaviest secondary.
468  if(A<0.1)
469  {
470    A = bufferA;
471    Z = bufferZ;
472  }
473 
474  // prepare the IsoResult.
475
476  std::ostringstream ost1;
477  ost1 <<Z<<"_"<<A;
478  G4String biff = ost1.str();
479  G4IsoResult * theResult = new G4IsoResult(biff, aNucleus);
480
481  return theResult;
482}
483
484G4double G4HadronicProcess::XBiasSurvivalProbability()
485{
486  G4double result = 0;
487  G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
488  G4double biasedProbability = 1.-std::exp(-nLTraversed);
489  G4double realProbability = 1-std::exp(-nLTraversed/aScaleFactor);
490  result = (biasedProbability-realProbability)/biasedProbability;
491  return result;
492}
493
494G4double G4HadronicProcess::XBiasSecondaryWeight()
495{
496  G4double result = 0;
497  G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
498  result = 
499     1./aScaleFactor*std::exp(-nLTraversed/aScaleFactor*(1-1./aScaleFactor));
500  return result;
501}
502
503void 
504G4HadronicProcess::FillTotalResult(G4HadFinalState * aR, const G4Track & aT)
505{
506  //  G4Nancheck go_wild;
507  theTotalResult->Clear();
508  theTotalResult->ProposeLocalEnergyDeposit(0.);
509  theTotalResult->Initialize(aT);
510  theTotalResult->SetSecondaryWeightByProcess(true);
511  theTotalResult->ProposeTrackStatus(fAlive);
512  G4double rotation = CLHEP::twopi*G4UniformRand();
513  G4ThreeVector it(0., 0., 1.);
514
515  /*
516  if(xBiasOn)
517  {
518    G4cout << "BiasDebug "<<GetProcessName()<<" "
519    <<aScaleFactor<<" "
520    <<XBiasSurvivalProbability()<<" "
521    <<XBiasSecondaryWeight()<<" "
522    <<G4endl;
523  }
524  */
525  // if(GetProcessName() != "LElastic") std::cout << "Debug -1 "<<aR->GetStatusChange()<<std::endl;
526  if(aR->GetStatusChange()==stopAndKill)
527  {
528    if( xBiasOn && G4UniformRand()<XBiasSurvivalProbability() )
529    {
530      theTotalResult->ProposeParentWeight( XBiasSurvivalProbability()*aT.GetWeight() );
531    }
532    else
533    {
534      theTotalResult->ProposeTrackStatus(fStopAndKill);
535      theTotalResult->ProposeEnergy( 0.0 );
536    }
537  }
538  else if(aR->GetStatusChange()!=stopAndKill )
539  {
540    if(aR->GetStatusChange()==suspend)
541    {
542      theTotalResult->ProposeTrackStatus(fSuspend);
543      if(xBiasOn)
544      {
545        DumpState(aT,"FillTotalResult");
546        G4Exception("G4HadronicProcess", "007", FatalException,
547        "Cannot cross-section bias a process that suspends tracks.");
548      }
549    } else if (aT.GetKineticEnergy() == 0) {
550      theTotalResult->ProposeTrackStatus(fStopButAlive);
551    }
552
553    if(xBiasOn && G4UniformRand()<XBiasSurvivalProbability())
554    {
555      theTotalResult->ProposeParentWeight( XBiasSurvivalProbability()*aT.GetWeight() );
556      G4double newWeight = aR->GetWeightChange()*aT.GetWeight();
557      /*
558      if(go_wild(aR->GetEnergyChange()))
559      {
560        G4Exception("G4HadronicProcess", "007", FatalException,
561        "surviving track received NaN energy."); 
562      }
563      if(go_wild(aR->GetMomentumChange().x()) ||
564        go_wild(aR->GetMomentumChange().y()) ||
565      go_wild(aR->GetMomentumChange().z()))
566      {
567        G4Exception("G4HadronicProcess", "007", FatalException,
568        "surviving track received NaN momentum."); 
569      }
570      */
571      G4double newM=aT.GetDefinition()->GetPDGMass();
572      G4double newE=aR->GetEnergyChange() + newM;
573      G4double newP=std::sqrt(newE*newE - newM*newM);
574      G4DynamicParticle * aNew = 
575      new G4DynamicParticle(aT.GetDefinition(), newE, newP*aR->GetMomentumChange());
576      G4HadSecondary * theSec = new G4HadSecondary(aNew, newWeight);
577      aR->AddSecondary(theSec);
578    }
579    else
580    {
581      G4double newWeight = aR->GetWeightChange()*aT.GetWeight();
582      theTotalResult->ProposeParentWeight(newWeight); // This is multiplicative
583      if(aR->GetEnergyChange()>-.5) 
584      {
585        /*
586        if(go_wild(aR->GetEnergyChange()))
587        {
588          G4Exception("G4HadronicProcess", "007", FatalException,
589          "track received NaN energy."); 
590        }
591        */
592        theTotalResult->ProposeEnergy(aR->GetEnergyChange());
593      }
594      G4LorentzVector newDirection(aR->GetMomentumChange().unit(), 1.);
595      newDirection*=aR->GetTrafoToLab();
596      theTotalResult->ProposeMomentumDirection(newDirection.vect());
597    }
598  }
599  else
600  {
601    G4cout << "Call for " << theInteraction->GetModelName() << G4endl;
602    G4cout << "Target Z= " 
603           << targetNucleus.GetZ() << "  A= " << targetNucleus.GetN() << G4endl;
604    DumpState(aT,"FillTotalResult");
605    G4Exception("G4HadronicProcess", "007", FatalException,
606    "use of unsupported track-status.");
607  }
608
609  if(GetProcessName() != "hElastic" && GetProcessName() != "HadronElastic"
610     &&  theTotalResult->GetTrackStatus()==fAlive
611     && aR->GetStatusChange()==isAlive)
612    {
613    // Use for debugging:   G4double newWeight = theTotalResult->GetParentWeight();
614
615    G4double newKE = std::max(DBL_MIN, aR->GetEnergyChange());
616    G4DynamicParticle* aNew = new G4DynamicParticle(aT.GetDefinition(), 
617                                                    aR->GetMomentumChange(), 
618                                                    newKE);
619    G4HadSecondary* theSec = new G4HadSecondary(aNew, 1.0);
620    aR->AddSecondary(theSec);
621    aR->SetStatusChange(stopAndKill);
622
623    theTotalResult->ProposeTrackStatus(fStopAndKill);
624    theTotalResult->ProposeEnergy( 0.0 );
625
626  }
627  theTotalResult->ProposeLocalEnergyDeposit(aR->GetLocalEnergyDeposit());
628  theTotalResult->SetNumberOfSecondaries(aR->GetNumberOfSecondaries());
629
630  if(aR->GetStatusChange() != stopAndKill)
631  {
632    G4double newM=aT.GetDefinition()->GetPDGMass();
633    G4double newE=aR->GetEnergyChange() + newM;
634    G4double newP=std::sqrt(newE*newE - newM*newM);
635    G4ThreeVector newPV = newP*aR->GetMomentumChange();
636    G4LorentzVector newP4(newE, newPV);
637    newP4.rotate(rotation, it);
638    newP4*=aR->GetTrafoToLab();
639    theTotalResult->ProposeMomentumDirection(newP4.vect().unit());
640  }
641
642  for(G4int i=0; i<aR->GetNumberOfSecondaries(); i++)
643  {
644    G4LorentzVector theM = aR->GetSecondary(i)->GetParticle()->Get4Momentum();
645    theM.rotate(rotation, it);
646    theM*=aR->GetTrafoToLab();
647    /*
648    if(go_wild(theM.e()))
649    {
650      G4Exception("G4HadronicProcess", "007", FatalException,
651      "secondary track received NaN energy."); 
652    }
653    if(go_wild(theM.x()) ||
654      go_wild(theM.y()) ||
655    go_wild(theM.z()))
656    {
657      G4Exception("G4HadronicProcess", "007", FatalException,
658      "secondary track received NaN momentum."); 
659    }
660    */
661    aR->GetSecondary(i)->GetParticle()->Set4Momentum(theM);
662    G4double time = aR->GetSecondary(i)->GetTime();
663    if(time<0) time = aT.GetGlobalTime();
664
665    G4Track* track = new G4Track(aR->GetSecondary(i)->GetParticle(),
666                                 time,
667                                 aT.GetPosition());
668
669    G4double newWeight = aT.GetWeight()*aR->GetSecondary(i)->GetWeight();
670    //static G4double pinelcount=0;
671    if(xBiasOn) newWeight *= XBiasSecondaryWeight();
672    // G4cout << "#### ParticleDebug "
673    // <<GetProcessName()<<" "
674    // <<aR->GetSecondary(i)->GetParticle()->GetDefinition()->GetParticleName()<<" "
675    // <<aScaleFactor<<" "
676    // <<XBiasSurvivalProbability()<<" "
677    // <<XBiasSecondaryWeight()<<" "
678    // <<aT.GetWeight()<<" "
679    // <<aR->GetSecondary(i)->GetWeight()<<" "
680    // <<aR->GetSecondary(i)->GetParticle()->Get4Momentum()<<" "
681    // <<G4endl;
682    track->SetWeight(newWeight);
683    /*
684    G4double trackDeb = track->GetKineticEnergy();
685    if( (  trackDeb<0
686      || (trackDeb>aT.GetKineticEnergy()+1*GeV) ) && getenv("GHADEnergyBalanceDebug") )
687      {
688        G4cout << "Debugging hadronic processes: "<<track->GetKineticEnergy()
689        <<" "<<aT.GetKineticEnergy()
690        <<" "<<GetProcessName()
691        <<" "<<aT.GetDefinition()->GetParticleName()
692        <<G4endl;
693      }
694    */
695    track->SetTouchableHandle(aT.GetTouchableHandle());
696    theTotalResult->AddSecondary(track);
697  }
698
699  aR->Clear();
700  return;
701}
702
703G4IsoParticleChange* G4HadronicProcess::GetIsotopeProductionInfo() 
704{ 
705  G4IsoParticleChange * anIsoResult = theIsoResult;
706  if(theIsoResult) theOldIsoResult = theIsoResult;
707  theIsoResult = 0;
708  return anIsoResult;
709}
710
711void G4HadronicProcess::BiasCrossSectionByFactor(G4double aScale) 
712{
713  xBiasOn = true;
714  aScaleFactor = aScale;
715  G4String it = GetProcessName(); 
716  if( (it != "PhotonInelastic") && 
717      (it != "ElectroNuclear") && 
718      (it != "PositronNuclear") )
719    {
720      G4Exception("G4HadronicProcess::BiasCrossSectionByFactor", "007", FatalException,
721                  "Cross-section biasing available only for gamma and electro nuclear reactions.");
722    }
723  if(aScale<100)
724    {
725      G4Exception("G4HadronicProcess::BiasCrossSectionByFactor", "001", JustWarning,
726                  "Cross-section bias readjusted to be above safe limit. New value is 100");
727      aScaleFactor = 100.;
728    }
729}
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
825/* end of file */
Note: See TracBrowser for help on using the repository browser.