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

Last change on this file since 968 was 962, checked in by garnier, 17 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.