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

Last change on this file since 1331 was 1315, checked in by garnier, 15 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.