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

Last change on this file since 1253 was 962, checked in by garnier, 15 years ago

update processes
File size: 22.8 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
48#include "G4HadLeadBias.hh"
49#include "G4HadronicException.hh"
50#include "G4HadReentrentException.hh"
51
52#include "G4HadronicWhiteBoard.hh"
53
54#include "G4HadSignalHandler.hh"
55
56#include "G4HadronicProcessStore.hh"
57
58#include <typeinfo>
59
60namespace G4HadronicProcess_local
61{
62  extern "C" void G4HadronicProcessHandler_1(int)
63  {
64    G4HadronicWhiteBoard::Instance().Dump();
65  }
66} 
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 ) :
82G4VDiscreteProcess( processName, aType)
83{ 
84  ModelingState = 0;
85  isoIsOnAnyway = -1;
86  theTotalResult = new G4ParticleChange();
87  theCrossSectionDataStore = new G4CrossSectionDataStore();
88  G4HadronicProcessStore::Instance()->Register(this);
89  aScaleFactor = 1;
90  xBiasOn = false;
91  G4HadronicProcess_debug_flag = false;
92  if(getenv("SwitchLeadBiasOn")) theBias.push_back(new G4HadLeadBias());
93}
94
95G4HadronicProcess::~G4HadronicProcess()
96{ 
97  G4HadronicProcessStore::Instance()->DeRegister(this);
98  delete theTotalResult;
99
100  std::for_each(theProductionModels.begin(),
101                theProductionModels.end(), G4Delete());
102  std::for_each(theBias.begin(), theBias.end(), G4Delete());
103 
104  delete theOldIsoResult; 
105  delete theIsoResult;
106  delete theCrossSectionDataStore;
107}
108
109void G4HadronicProcess::RegisterMe( G4HadronicInteraction *a )
110{ 
111  try{GetManagerPointer()->RegisterMe( a );}   
112  catch(G4HadronicException & aE)
113  {
114    aE.Report(std::cout);
115    G4Exception("G4HadronicProcess", "007", FatalException,
116    "Could not register G4HadronicInteraction");
117  }
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);
131}
132
133G4double G4HadronicProcess::
134GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *)
135{ 
136  try
137  {
138    ModelingState = 1;   
139    theLastCrossSection = aScaleFactor* 
140      theCrossSectionDataStore->GetCrossSection(aTrack.GetDynamicParticle(), 
141                                                aTrack.GetMaterial());
142  }
143  catch(G4HadronicException aR)
144  { 
145    aR.Report(G4cout); 
146    G4Exception("G4HadronicProcess", "007", FatalException,
147    "G4HadronicProcess::GetMeanFreePath failed");
148  } 
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 &)
165{
166  // Debugging stuff
167  if(G4HadronicProcess_debug_flag) 
168         std::cout << "@@@@ hadronic process start "<< std::endl;
169  // G4cout << theNumberOfInteractionLengthLeft<<G4endl;
170#ifndef G4HadSignalHandler_off
171  G4HadSignalHandler aHandler(G4HadronicProcess_local::G4HadronicProcessHandler_1);
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
184    theTotalResult->Clear();
185    theTotalResult->Initialize(aTrack);
186    return theTotalResult;
187  }
188
189  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
190  G4Material* aMaterial = aTrack.GetMaterial();
191  G4double originalEnergy = aParticle->GetKineticEnergy();
192  G4double kineticEnergy = originalEnergy;
193
194  /*
195  // It is not needed with standard NaN check
196  // More debugging 
197  G4Nancheck go_wild;
198  if(go_wild(originalEnergy) ||
199    go_wild(aParticle->Get4Momentum().x()) ||
200  go_wild(aParticle->Get4Momentum().y()) ||
201  go_wild(aParticle->Get4Momentum().z()) ||
202  go_wild(aParticle->Get4Momentum().t())
203  )
204  {
205    G4Exception("G4HadronicProcess", "001", JustWarning, "NaN in input energy or momentum - bailing out.");
206    theTotalResult->Clear();
207    theTotalResult->Initialize(aTrack);
208    return theTotalResult;
209  }
210  */
211
212  // Get kinetic energy per nucleon for ions
213  if(aParticle->GetDefinition()->GetBaryonNumber() > 1.5)
214          kineticEnergy/=aParticle->GetDefinition()->GetBaryonNumber();
215
216  G4Element* anElement = 0;
217  try
218  {
219    //    anElement = ChooseAandZ( aParticle, aMaterial );
220    anElement = theCrossSectionDataStore->SampleZandA(aParticle, 
221                                                      aMaterial, 
222                                                      targetNucleus);
223  }
224  catch(G4HadronicException & aR)
225  {
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;
232    G4Exception("G4HadronicProcess", "007", FatalException,
233    "PostStepDoIt failed on element selection.");
234  }
235
236  try
237  {
238    theInteraction = ChooseHadronicInteraction( kineticEnergy,
239    aMaterial, anElement );
240  }
241  catch(G4HadronicException & aE)
242  {
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;
249    G4Exception("G4HadronicProcess", "007", FatalException,
250    "ChooseHadronicInteraction failed.");
251  }
252
253  // Initialize the hadronic projectile from the track
254
255  G4HadProjectile thePro(aTrack);
256 
257  G4HadFinalState* result = 0;
258  G4int reentryCount = 0;
259
260  do
261  {
262    try
263    {
264      // Call the interaction
265      result = theInteraction->ApplyYourself( thePro, targetNucleus);
266    }
267    catch(G4HadReentrentException aR)
268    {
269      G4HadronicWhiteBoard & theBoard = G4HadronicWhiteBoard::Instance();
270      theBoard.SetProjectile(thePro);
271      theBoard.SetTargetNucleus(targetNucleus);
272      theBoard.SetProcessName(GetProcessName());
273      theBoard.SetModelName(theInteraction->GetModelName());
274
275      aR.Report(G4cout);
276      G4cout << " G4HadronicProcess re-entering the ApplyYourself call for "
277             <<G4endl;
278      G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
279      G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
280      G4cout << " - Particle type = " 
281             << aParticle->GetDefinition()->GetParticleName() << G4endl;
282      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      }
288      G4Exception("G4HadronicProcess", "007", FatalException,
289      "GetHadronicProcess: PostStepDoIt failed (Reentering ApplyYourself not yet supported.)"); 
290    }
291    catch(G4HadronicException aR)
292    {
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;
306      G4Exception("G4HadronicProcess", "007", FatalException,
307      "PostStepDoIt failed.");
308    }
309  }
310  while(!result);
311
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  }
317
318  // 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  }
326
327  result->SetTrafoToLab(thePro.GetTrafoToLab());
328
329  if(getenv("HadronicDoitLogging") )
330  {
331    G4cout << "HadronicDoitLogging "
332           << GetProcessName() <<" "
333           << aParticle->GetDefinition()->GetPDGEncoding()<<" "
334           << originalEnergy<<" "
335           << aParticle->GetMomentum()<<" "
336           << targetNucleus.GetN()<<" "
337           << targetNucleus.GetZ()<<" "
338           << G4endl;
339  }
340
341  ClearNumberOfInteractionLengthLeft();
342  if(isoIsOnAnyway!=-1)
343  {
344    if(isoIsEnabled||isoIsOnAnyway)
345    {
346      result = DoIsotopeCounting(result, aTrack, targetNucleus);
347    }
348  }
349 
350  G4double e=aTrack.GetKineticEnergy();
351  ModelingState = 0;
352  if(e<5*GeV)
353  {
354    for(size_t i=0; i<theBias.size(); i++)
355    {
356      result = theBias[i]->Bias(result);
357    }
358  }
359
360  // Put hadronic final state particles into G4ParticleChange
361
362  FillTotalResult(result, aTrack);
363  if(G4HadronicProcess_debug_flag) 
364    std::cout << "@@@@ hadronic process end "<< std::endl;
365   
366  return theTotalResult;
367}
368
369
370G4HadFinalState* 
371G4HadronicProcess::DoIsotopeCounting(G4HadFinalState * aResult,
372                                     const G4Track & aTrack,
373                                     const G4Nucleus & aNucleus)
374{
375  // get the PC from iso-production
376  delete theOldIsoResult;
377  theOldIsoResult = 0;
378  delete theIsoResult;
379  theIsoResult = new G4IsoParticleChange;
380  G4bool done = false;
381  G4IsoResult * anIsoResult = 0;
382  for(unsigned int i=0; i<theProductionModels.size(); i++)
383  {
384    anIsoResult = theProductionModels[i]->GetIsotope(aTrack, aNucleus);
385    if(anIsoResult!=0)
386    {
387      done = true;
388      break;
389    }
390  }
391
392  // If no production models active, use default iso production
393  if(!done) anIsoResult = ExtractResidualNucleus(aTrack, aNucleus, aResult); 
394
395  // Add all info explicitely and add typename from model called.
396  theIsoResult->SetIsotope(anIsoResult->GetIsotope());
397  theIsoResult->SetProductionPosition(aTrack.GetPosition());
398  theIsoResult->SetProductionTime(aTrack.GetGlobalTime());
399  theIsoResult->SetParentParticle(*aTrack.GetDynamicParticle());
400  theIsoResult->SetMotherNucleus(anIsoResult->GetMotherNucleus());
401  theIsoResult->SetProducer(typeid(*theInteraction).name());
402 
403  delete anIsoResult;
404
405  // If isotope production is enabled the GetIsotopeProductionInfo()
406  // method must be called or else a memory leak will result
407  //
408  // The following code will fix the memory leak, but remove the
409  // isotope information:
410  //
411  //  if(theIsoResult) {
412  //    delete theIsoResult;
413  //    theIsoResult = 0;
414  //  }
415 
416  return aResult;
417}
418
419G4IsoResult* 
420G4HadronicProcess::ExtractResidualNucleus(const G4Track&,
421                                          const G4Nucleus& aNucleus,
422                                          G4HadFinalState* aResult)
423{
424  G4double A = aNucleus.GetN();
425  G4double Z = aNucleus.GetZ();
426  G4double bufferA = 0;
427  G4double bufferZ = 0;
428 
429  // loop over aResult, and decrement A, Z accordingly
430  // cash the max
431  for(G4int i=0; i<aResult->GetNumberOfSecondaries(); i++)
432  {
433    G4HadSecondary* aSecTrack = aResult->GetSecondary(i);
434    if(bufferA<aSecTrack->GetParticle()->GetDefinition()->GetBaryonNumber())
435    {
436      bufferA = aSecTrack->GetParticle()->GetDefinition()->GetBaryonNumber();
437      bufferZ = aSecTrack->GetParticle()->GetDefinition()->GetPDGCharge();
438    }
439    Z-=aSecTrack->GetParticle()->GetDefinition()->GetPDGCharge();
440    A-=aSecTrack->GetParticle()->GetDefinition()->GetBaryonNumber();
441  }
442 
443  // if the fragment was part of the final state, it is
444  // assumed to be the heaviest secondary.
445  if(A<0.1)
446  {
447    A = bufferA;
448    Z = bufferZ;
449  }
450 
451  // prepare the IsoResult.
452
453  std::ostringstream ost1;
454  ost1 <<Z<<"_"<<A;
455  G4String biff = ost1.str();
456  G4IsoResult * theResult = new G4IsoResult(biff, aNucleus);
457
458  return theResult;
459}
460
461G4double G4HadronicProcess::XBiasSurvivalProbability()
462{
463  G4double result = 0;
464  G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
465  G4double biasedProbability = 1.-std::exp(-nLTraversed);
466  G4double realProbability = 1-std::exp(-nLTraversed/aScaleFactor);
467  result = (biasedProbability-realProbability)/biasedProbability;
468  return result;
469}
470
471G4double G4HadronicProcess::XBiasSecondaryWeight()
472{
473  G4double result = 0;
474  G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
475  result = 
476     1./aScaleFactor*std::exp(-nLTraversed/aScaleFactor*(1-1./aScaleFactor));
477  return result;
478}
479
480void 
481G4HadronicProcess::FillTotalResult(G4HadFinalState * aR, const G4Track & aT)
482{
483  //  G4Nancheck go_wild;
484  theTotalResult->Clear();
485  theTotalResult->ProposeLocalEnergyDeposit(0.);
486  theTotalResult->Initialize(aT);
487  theTotalResult->SetSecondaryWeightByProcess(true);
488  theTotalResult->ProposeTrackStatus(fAlive);
489  G4double rotation = 2.*pi*G4UniformRand();
490  G4ThreeVector it(0., 0., 1.);
491
492  /*
493  if(xBiasOn)
494  {
495    G4cout << "BiasDebug "<<GetProcessName()<<" "
496    <<aScaleFactor<<" "
497    <<XBiasSurvivalProbability()<<" "
498    <<XBiasSecondaryWeight()<<" "
499    <<G4endl;
500  }
501  */
502  // if(GetProcessName() != "LElastic") std::cout << "Debug -1 "<<aR->GetStatusChange()<<std::endl;
503  if(aR->GetStatusChange()==stopAndKill)
504  {
505    if( xBiasOn && G4UniformRand()<XBiasSurvivalProbability() )
506    {
507      theTotalResult->ProposeParentWeight( XBiasSurvivalProbability()*aT.GetWeight() );
508    }
509    else
510    {
511      theTotalResult->ProposeTrackStatus(fStopAndKill);
512      theTotalResult->ProposeEnergy( 0.0 );
513    }
514  }
515  else if(aR->GetStatusChange()!=stopAndKill )
516  {
517    if(aR->GetStatusChange()==suspend)
518    {
519      theTotalResult->ProposeTrackStatus(fSuspend);
520      if(xBiasOn)
521      {
522        G4Exception("G4HadronicProcess", "007", FatalException,
523        "Cannot cross-section bias a process that suspends tracks.");
524      }
525    } else if (aT.GetKineticEnergy() == 0) {
526      theTotalResult->ProposeTrackStatus(fStopButAlive);
527    }
528
529    if(xBiasOn && G4UniformRand()<XBiasSurvivalProbability())
530    {
531      theTotalResult->ProposeParentWeight( XBiasSurvivalProbability()*aT.GetWeight() );
532      G4double newWeight = aR->GetWeightChange()*aT.GetWeight();
533      /*
534      if(go_wild(aR->GetEnergyChange()))
535      {
536        G4Exception("G4HadronicProcess", "007", FatalException,
537        "surviving track received NaN energy."); 
538      }
539      if(go_wild(aR->GetMomentumChange().x()) ||
540        go_wild(aR->GetMomentumChange().y()) ||
541      go_wild(aR->GetMomentumChange().z()))
542      {
543        G4Exception("G4HadronicProcess", "007", FatalException,
544        "surviving track received NaN momentum."); 
545      }
546      */
547      G4double newM=aT.GetDefinition()->GetPDGMass();
548      G4double newE=aR->GetEnergyChange() + newM;
549      G4double newP=std::sqrt(newE*newE - newM*newM);
550      G4DynamicParticle * aNew = 
551      new G4DynamicParticle(aT.GetDefinition(), newE, newP*aR->GetMomentumChange());
552      G4HadSecondary * theSec = new G4HadSecondary(aNew, newWeight);
553      aR->AddSecondary(theSec);
554    }
555    else
556    {
557      G4double newWeight = aR->GetWeightChange()*aT.GetWeight();
558      theTotalResult->ProposeParentWeight(newWeight); // This is multiplicative
559      if(aR->GetEnergyChange()>-.5) 
560      {
561        /*
562        if(go_wild(aR->GetEnergyChange()))
563        {
564          G4Exception("G4HadronicProcess", "007", FatalException,
565          "track received NaN energy."); 
566        }
567        */
568        theTotalResult->ProposeEnergy(aR->GetEnergyChange());
569      }
570      G4LorentzVector newDirection(aR->GetMomentumChange().unit(), 1.);
571      newDirection*=aR->GetTrafoToLab();
572      theTotalResult->ProposeMomentumDirection(newDirection.vect());
573    }
574  }
575  else
576  {
577    G4cerr << "Track status is "<< aR->GetStatusChange()<<G4endl;
578    G4Exception("G4HadronicProcess", "007", FatalException,
579    "use of unsupported track-status.");
580  }
581
582  if(GetProcessName() != "hElastic" && GetProcessName() != "HadronElastic"
583     &&  theTotalResult->GetTrackStatus()==fAlive
584     && aR->GetStatusChange()==isAlive)
585    {
586    // Use for debugging:   G4double newWeight = theTotalResult->GetParentWeight();
587
588    G4double newKE = std::max(DBL_MIN, aR->GetEnergyChange());
589    G4DynamicParticle* aNew = new G4DynamicParticle(aT.GetDefinition(), 
590                                                    aR->GetMomentumChange(), 
591                                                    newKE);
592    G4HadSecondary* theSec = new G4HadSecondary(aNew, 1.0);
593    aR->AddSecondary(theSec);
594    aR->SetStatusChange(stopAndKill);
595
596    theTotalResult->ProposeTrackStatus(fStopAndKill);
597    theTotalResult->ProposeEnergy( 0.0 );
598
599  }
600  theTotalResult->ProposeLocalEnergyDeposit(aR->GetLocalEnergyDeposit());
601  theTotalResult->SetNumberOfSecondaries(aR->GetNumberOfSecondaries());
602
603  if(aR->GetStatusChange() != stopAndKill)
604  {
605    G4double newM=aT.GetDefinition()->GetPDGMass();
606    G4double newE=aR->GetEnergyChange() + newM;
607    G4double newP=std::sqrt(newE*newE - newM*newM);
608    G4ThreeVector newPV = newP*aR->GetMomentumChange();
609    G4LorentzVector newP4(newE, newPV);
610    newP4.rotate(rotation, it);
611    newP4*=aR->GetTrafoToLab();
612    theTotalResult->ProposeMomentumDirection(newP4.vect().unit());
613  }
614
615  for(G4int i=0; i<aR->GetNumberOfSecondaries(); i++)
616  {
617    G4LorentzVector theM = aR->GetSecondary(i)->GetParticle()->Get4Momentum();
618    theM.rotate(rotation, it);
619    theM*=aR->GetTrafoToLab();
620    /*
621    if(go_wild(theM.e()))
622    {
623      G4Exception("G4HadronicProcess", "007", FatalException,
624      "secondary track received NaN energy."); 
625    }
626    if(go_wild(theM.x()) ||
627      go_wild(theM.y()) ||
628    go_wild(theM.z()))
629    {
630      G4Exception("G4HadronicProcess", "007", FatalException,
631      "secondary track received NaN momentum."); 
632    }
633    */
634    aR->GetSecondary(i)->GetParticle()->Set4Momentum(theM);
635    G4double time = aR->GetSecondary(i)->GetTime();
636    if(time<0) time = aT.GetGlobalTime();
637
638    G4Track* track = new G4Track(aR->GetSecondary(i)->GetParticle(),
639                                 time,
640                                 aT.GetPosition());
641
642    G4double newWeight = aT.GetWeight()*aR->GetSecondary(i)->GetWeight();
643    //static G4double pinelcount=0;
644    if(xBiasOn) newWeight *= XBiasSecondaryWeight();
645    // G4cout << "#### ParticleDebug "
646    // <<GetProcessName()<<" "
647    // <<aR->GetSecondary(i)->GetParticle()->GetDefinition()->GetParticleName()<<" "
648    // <<aScaleFactor<<" "
649    // <<XBiasSurvivalProbability()<<" "
650    // <<XBiasSecondaryWeight()<<" "
651    // <<aT.GetWeight()<<" "
652    // <<aR->GetSecondary(i)->GetWeight()<<" "
653    // <<aR->GetSecondary(i)->GetParticle()->Get4Momentum()<<" "
654    // <<G4endl;
655    track->SetWeight(newWeight);
656    /*
657    G4double trackDeb = track->GetKineticEnergy();
658    if( (  trackDeb<0
659      || (trackDeb>aT.GetKineticEnergy()+1*GeV) ) && getenv("GHADEnergyBalanceDebug") )
660      {
661        G4cout << "Debugging hadronic processes: "<<track->GetKineticEnergy()
662        <<" "<<aT.GetKineticEnergy()
663        <<" "<<GetProcessName()
664        <<" "<<aT.GetDefinition()->GetParticleName()
665        <<G4endl;
666      }
667    */
668    track->SetTouchableHandle(aT.GetTouchableHandle());
669    theTotalResult->AddSecondary(track);
670  }
671
672  aR->Clear();
673  return;
674}
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}
703/* end of file */
Note: See TracBrowser for help on using the repository browser.