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

Last change on this file since 900 was 819, checked in by garnier, 17 years ago

import all except CVS
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// #include "G4EffectiveCharge.hh"
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#include "G4HadronicInteractionWrapper.hh"
52
53#include "G4HadSignalHandler.hh"
54
55#include <typeinfo>
56
57namespace G4HadronicProcess_local
58{
59 extern "C" void G4HadronicProcessHandler_1(int)
60 {
61 G4HadronicWhiteBoard::Instance().Dump();
62 }
63}
64
65G4IsoParticleChange * G4HadronicProcess::theIsoResult = 0;
66G4IsoParticleChange * G4HadronicProcess::theOldIsoResult = 0;
67G4bool G4HadronicProcess::isoIsEnabled = true;
68
69void G4HadronicProcess::
70EnableIsotopeProductionGlobally() {isoIsEnabled = true;}
71
72void G4HadronicProcess::
73DisableIsotopeProductionGlobally() {isoIsEnabled = false;}
74
75G4HadronicProcess::G4HadronicProcess( const G4String &processName,
76 G4ProcessType aType ) :
77G4VDiscreteProcess( processName, aType)
78{
79 ModelingState = 0;
80 isoIsOnAnyway = -1;
81 theTotalResult = new G4ParticleChange();
82 theCrossSectionDataStore = new G4CrossSectionDataStore();
83 aScaleFactor = 1;
84 xBiasOn = false;
85 if(getenv("SwitchLeadBiasOn")) theBias.push_back(new G4HadLeadBias());
86}
87
88G4HadronicProcess::~G4HadronicProcess()
89{
90 delete theTotalResult;
91
92 std::for_each(theProductionModels.begin(),
93 theProductionModels.end(), G4Delete());
94 std::for_each(theBias.begin(), theBias.end(), G4Delete());
95
96 delete theOldIsoResult; delete theIsoResult;
97 delete theCrossSectionDataStore;
98}
99
100void G4HadronicProcess::RegisterMe( G4HadronicInteraction *a )
101{
102 try{GetManagerPointer()->RegisterMe( a );}
103 catch(G4HadronicException & aE)
104 {
105 aE.Report(std::cout);
106 G4Exception("G4HadronicProcess", "007", FatalException,
107 "Could not register G4HadronicInteraction");
108 }
109}
110
111G4double G4HadronicProcess::
112GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *)
113{
114 G4double sigma = 0.0;
115 try
116 {
117 const G4DynamicParticle *aParticle = aTrack.GetDynamicParticle();
118 if( !IsApplicable(*aParticle->GetDefinition()))
119 {
120 G4cout << "Unrecoverable error: "<<G4endl;
121 G4ProcessManager * it = aParticle->GetDefinition()->GetProcessManager();
122 G4ProcessVector * itv = it->GetProcessList();
123 G4cout <<aParticle->GetDefinition()->GetParticleName()<<
124 " has the following processes:"<<G4endl;
125 for(G4int i=0; i<itv->size(); i++)
126 {
127 G4cout <<" "<<(*itv)[i]->GetProcessName()<<G4endl;
128 }
129 G4cout << "for kinetic energy "<<aParticle->GetKineticEnergy()<<G4endl;
130 G4cout << "and material "<<aTrack.GetMaterial()->GetName()<<G4endl;
131 G4Exception("G4HadronicProcess", "007", FatalException,
132 std::string(this->GetProcessName()+
133 " was called for "+
134 aParticle->GetDefinition()->GetParticleName()).c_str() );
135 }
136 G4Material *aMaterial = aTrack.GetMaterial();
137 ModelingState = 1;
138
139 sigma = theCrossSectionDataStore->GetCrossSection(aParticle, aMaterial);
140
141 sigma *= aScaleFactor;
142 theLastCrossSection = sigma;
143 }
144 catch(G4HadronicException aR)
145 {
146 aR.Report(G4cout);
147 G4Exception("G4HadronicProcess", "007", FatalException,
148 "G4HadronicProcess::GetMeanFreePath failed");
149 }
150 if( sigma > 0.0 )
151 return 1.0/sigma;
152 else
153 return DBL_MAX;
154}
155
156
157G4Element* G4HadronicProcess::ChooseAandZ(
158const G4DynamicParticle *aParticle, const G4Material *aMaterial )
159{
160 std::pair<G4double, G4double> ZA =
161 theCrossSectionDataStore->SelectRandomIsotope(aParticle, aMaterial);
162 G4double ZZ = ZA.first;
163 G4double AA = ZA.second;
164
165 targetNucleus.SetParameters(AA, ZZ);
166
167 const G4int numberOfElements = aMaterial->GetNumberOfElements();
168 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
169 G4Element* chosen = 0;
170 for (G4int i = 0; i < numberOfElements; i++) {
171 chosen = (*theElementVector)[i];
172 if (chosen->GetZ() == ZZ) break;
173 }
174 return chosen;
175}
176
177
178struct G4Nancheck{ bool operator()(G4double aV){return (!(aV<1))&&(!(aV>-1));}};
179
180G4VParticleChange *G4HadronicProcess::GeneralPostStepDoIt(
181const G4Track &aTrack, const G4Step &)
182{
183 // Debugging stuff
184
185 bool G4HadronicProcess_debug_flag = false;
186 if(getenv("G4HadronicProcess_debug")) G4HadronicProcess_debug_flag = true;
187 if(G4HadronicProcess_debug_flag)
188 std::cout << "@@@@ hadronic process start "<< std::endl;
189 // G4cout << theNumberOfInteractionLengthLeft<<G4endl;
190 #ifndef G4HadSignalHandler_off
191 G4HadSignalHandler aHandler(G4HadronicProcess_local::G4HadronicProcessHandler_1);
192 #endif
193
194 if(aTrack.GetTrackStatus() != fAlive && aTrack.GetTrackStatus() != fSuspend)
195 {
196 G4cerr << "G4HadronicProcess: track in unusable state - "
197 <<aTrack.GetTrackStatus()<<G4endl;
198 G4cerr << "G4HadronicProcess: returning unchanged track "<<G4endl;
199 G4Exception("G4HadronicProcess", "001", JustWarning, "bailing out");
200 theTotalResult->Clear();
201 theTotalResult->Initialize(aTrack);
202 return theTotalResult;
203 }
204
205 const G4DynamicParticle *aParticle = aTrack.GetDynamicParticle();
206 G4Material *aMaterial = aTrack.GetMaterial();
207 G4double originalEnergy = aParticle->GetKineticEnergy();
208 G4double kineticEnergy = originalEnergy;
209
210 // More debugging
211
212 G4Nancheck go_wild;
213 if(go_wild(originalEnergy) ||
214 go_wild(aParticle->Get4Momentum().x()) ||
215 go_wild(aParticle->Get4Momentum().y()) ||
216 go_wild(aParticle->Get4Momentum().z()) ||
217 go_wild(aParticle->Get4Momentum().t())
218 )
219 {
220 G4Exception("G4HadronicProcess", "001", JustWarning, "NaN in input energy or momentum - bailing out.");
221 theTotalResult->Clear();
222 theTotalResult->Initialize(aTrack);
223 return theTotalResult;
224 }
225
226 // Get kinetic energy per nucleon for ions
227
228 if(aParticle->GetDefinition()->GetBaryonNumber() > 1.5)
229 kineticEnergy/=aParticle->GetDefinition()->GetBaryonNumber();
230
231 G4Element* anElement = 0;
232 try
233 {
234 anElement = ChooseAandZ( aParticle, aMaterial );
235 }
236 catch(G4HadronicException & aR)
237 {
238 aR.Report(G4cout);
239 G4cout << "Unrecoverable error for:"<<G4endl;
240 G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
241 G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
242 G4cout << " - Particle type = "
243 <<aParticle->GetDefinition()->GetParticleName()<<G4endl;
244 G4Exception("G4HadronicProcess", "007", FatalException,
245 "GeneralPostStepDoIt failed on element selection.");
246 }
247
248 try
249 {
250 theInteraction = ChooseHadronicInteraction( kineticEnergy,
251 aMaterial, anElement );
252 }
253 catch(G4HadronicException & aE)
254 {
255 aE.Report(std::cout);
256 G4cout << "Unrecoverable error for:"<<G4endl;
257 G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
258 G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
259 G4cout << " - Particle type = "
260 << aParticle->GetDefinition()->GetParticleName()<<G4endl;
261 G4Exception("G4HadronicProcess", "007", FatalException,
262 "ChooseHadronicInteraction failed.");
263 }
264
265 // Initialize the hadronic projectile from the track
266
267 G4HadProjectile thePro(aTrack);
268
269 G4HadFinalState* result = 0;
270 G4int reentryCount = 0;
271
272 do
273 {
274 try
275 {
276 // Call the interaction
277
278 G4HadronicInteractionWrapper aW;
279 result = aW.ApplyInteraction(thePro, targetNucleus, theInteraction,
280 GetProcessName(),
281 theInteraction->GetModelName());
282 }
283 catch(G4HadReentrentException aR)
284 {
285 aR.Report(G4cout);
286 G4cout << " G4HadronicProcess re-entering the ApplyYourself call for "
287 <<G4endl;
288 G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
289 G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
290 G4cout << " - Particle type = "
291 << aParticle->GetDefinition()->GetParticleName() << G4endl;
292 result = 0; // here would still be leaking...
293 if(reentryCount>100)
294 {
295 G4Exception("G4HadronicProcess", "007", FatalException,
296 "GetHadronicProcess: Reentering ApplyYourself too often - GeneralPostStepDoIt failed.");
297 }
298 G4Exception("G4HadronicProcess", "007", FatalException,
299 "GetHadronicProcess: GeneralPostStepDoIt failed (Reentering ApplyYourself not yet supported.)");
300 }
301 catch(G4HadronicException aR)
302 {
303 aR.Report(G4cout);
304 G4cout << " G4HadronicProcess failed in ApplyYourself call for"
305 << G4endl;
306 G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
307 G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
308 G4cout << " - Particle type = "
309 << aParticle->GetDefinition()->GetParticleName() << G4endl;
310 G4Exception("G4HadronicProcess", "007", FatalException,
311 "GeneralPostStepDoIt failed.");
312 }
313 }
314 while(!result);
315
316 if(!ModelingState && !getenv("BypassAllSafetyChecks") )
317 {
318 G4cout << "ERROR IN EXECUTION -- HADRONIC PROCESS STATE NOT VALID"<<G4endl;
319 G4cout << "Result will be of undefined quality."<<G4endl;
320 }
321
322 // NOT USED ?? Projectile particle has changed character during interaction
323 if(result->GetStatusChange() == isAlive &&
324 thePro.GetDefinition() != aTrack.GetDefinition())
325 {
326 G4DynamicParticle * aP =
327 const_cast<G4DynamicParticle *>(aTrack.GetDynamicParticle());
328 aP->SetDefinition(const_cast<G4ParticleDefinition *>(thePro.GetDefinition()));
329 }
330
331 result->SetTrafoToLab(thePro.GetTrafoToLab());
332
333 /*
334 // Loop over charged ion secondaries
335
336 for(G4int i=0; i<result->GetNumberOfSecondaries(); i++)
337 {
338 G4DynamicParticle* aSecTrack = result->GetSecondary(i)->GetParticle();
339 if(aSecTrack->GetDefinition()->GetPDGCharge()>1.5)
340 {
341 G4EffectiveCharge aCalculator;
342 G4double charge =
343 aCalculator.GetCharge(aMaterial, aSecTrack->GetKineticEnergy(),
344 aSecTrack->GetDefinition()->GetPDGMass(),
345 aSecTrack->GetDefinition()->GetPDGCharge());
346 if(getenv("GHADChargeDebug"))
347 {
348 std::cout << "Recoil fractional charge is "
349 << charge/aSecTrack->GetDefinition()->GetPDGCharge()<<" "
350 << charge <<" "<<aSecTrack->GetDefinition()->GetPDGCharge()<<std::endl;
351 }
352 aSecTrack->SetCharge(charge);
353 }
354 }
355 */
356
357 if(getenv("HadronicDoitLogging") )
358 {
359 G4cout << "HadronicDoitLogging "
360 << GetProcessName() <<" "
361 << aParticle->GetDefinition()->GetPDGEncoding()<<" "
362 << originalEnergy<<" "
363 << aParticle->GetMomentum()<<" "
364 << targetNucleus.GetN()<<" "
365 << targetNucleus.GetZ()<<" "
366 << G4endl;
367 }
368
369 ClearNumberOfInteractionLengthLeft();
370 if(isoIsOnAnyway!=-1)
371 {
372 if(isoIsEnabled||isoIsOnAnyway)
373 {
374 result = DoIsotopeCounting(result, aTrack, targetNucleus);
375 }
376 }
377
378 G4double e=aTrack.GetKineticEnergy();
379 ModelingState = 0;
380 if(e<5*GeV)
381 {
382 for(size_t i=0; i<theBias.size(); i++)
383 {
384 result = theBias[i]->Bias(result);
385 }
386 }
387
388 // Put hadronic final state particles into G4ParticleChange
389
390 FillTotalResult(result, aTrack);
391 if(G4HadronicProcess_debug_flag)
392 std::cout << "@@@@ hadronic process end "<< std::endl;
393
394 return theTotalResult;
395}
396
397
398G4HadFinalState*
399G4HadronicProcess::DoIsotopeCounting(G4HadFinalState * aResult,
400 const G4Track & aTrack,
401 const G4Nucleus & aNucleus)
402{
403 // get the PC from iso-production
404 delete theOldIsoResult;
405 theOldIsoResult = 0;
406 delete theIsoResult;
407 theIsoResult = new G4IsoParticleChange;
408 G4bool done = false;
409 G4IsoResult * anIsoResult = 0;
410 for(unsigned int i=0; i<theProductionModels.size(); i++)
411 {
412 anIsoResult = theProductionModels[i]->GetIsotope(aTrack, aNucleus);
413 if(anIsoResult!=0)
414 {
415 done = true;
416 break;
417 }
418 }
419
420 // If no production models active, use default iso production
421 if(!done) anIsoResult = ExtractResidualNucleus(aTrack, aNucleus, aResult);
422
423 // Add all info explicitely and add typename from model called.
424 theIsoResult->SetIsotope(anIsoResult->GetIsotope());
425 theIsoResult->SetProductionPosition(aTrack.GetPosition());
426 theIsoResult->SetProductionTime(aTrack.GetGlobalTime());
427 theIsoResult->SetParentParticle(*aTrack.GetDynamicParticle());
428 theIsoResult->SetMotherNucleus(anIsoResult->GetMotherNucleus());
429 theIsoResult->SetProducer(typeid(*theInteraction).name());
430
431 delete anIsoResult;
432
433 // If isotope production is enabled the GetIsotopeProductionInfo()
434 // method must be called or else a memory leak will result
435 //
436 // The following code will fix the memory leak, but remove the
437 // isotope information:
438 //
439 // if(theIsoResult) {
440 // delete theIsoResult;
441 // theIsoResult = 0;
442 // }
443
444 return aResult;
445}
446
447G4IsoResult*
448G4HadronicProcess::ExtractResidualNucleus(const G4Track&,
449 const G4Nucleus& aNucleus,
450 G4HadFinalState* aResult)
451{
452 G4double A = aNucleus.GetN();
453 G4double Z = aNucleus.GetZ();
454 G4double bufferA = 0;
455 G4double bufferZ = 0;
456
457 // loop over aResult, and decrement A, Z accordingly
458 // cash the max
459 for(G4int i=0; i<aResult->GetNumberOfSecondaries(); i++)
460 {
461 G4HadSecondary* aSecTrack = aResult->GetSecondary(i);
462 if(bufferA<aSecTrack->GetParticle()->GetDefinition()->GetBaryonNumber())
463 {
464 bufferA = aSecTrack->GetParticle()->GetDefinition()->GetBaryonNumber();
465 bufferZ = aSecTrack->GetParticle()->GetDefinition()->GetPDGCharge();
466 }
467 Z-=aSecTrack->GetParticle()->GetDefinition()->GetPDGCharge();
468 A-=aSecTrack->GetParticle()->GetDefinition()->GetBaryonNumber();
469 }
470
471 // if the fragment was part of the final state, it is
472 // assumed to be the heaviest secondary.
473 if(A<0.1)
474 {
475 A = bufferA;
476 Z = bufferZ;
477 }
478
479 // prepare the IsoResult.
480
481 std::ostringstream ost1;
482 ost1 <<Z<<"_"<<A;
483 G4String biff = ost1.str();
484 G4IsoResult * theResult = new G4IsoResult(biff, aNucleus);
485
486 return theResult;
487}
488
489G4double G4HadronicProcess::XBiasSurvivalProbability()
490{
491 G4double result = 0;
492 G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
493 G4double biasedProbability = 1.-std::exp(-nLTraversed);
494 G4double realProbability = 1-std::exp(-nLTraversed/aScaleFactor);
495 result = (biasedProbability-realProbability)/biasedProbability;
496 return result;
497}
498
499G4double G4HadronicProcess::XBiasSecondaryWeight()
500{
501 G4double result = 0;
502 G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
503 result =
504 1./aScaleFactor*std::exp(-nLTraversed/aScaleFactor*(1-1./aScaleFactor));
505 return result;
506}
507
508void
509G4HadronicProcess::FillTotalResult(G4HadFinalState * aR, const G4Track & aT)
510{
511 G4Nancheck go_wild;
512 theTotalResult->Clear();
513 theTotalResult->ProposeLocalEnergyDeposit(0.);
514 theTotalResult->Initialize(aT);
515 theTotalResult->SetSecondaryWeightByProcess(true);
516 theTotalResult->ProposeTrackStatus(fAlive);
517 G4double rotation = 2.*pi*G4UniformRand();
518 G4ThreeVector it(0., 0., 1.);
519
520 /*
521 if(xBiasOn)
522 {
523 G4cout << "BiasDebug "<<GetProcessName()<<" "
524 <<aScaleFactor<<" "
525 <<XBiasSurvivalProbability()<<" "
526 <<XBiasSecondaryWeight()<<" "
527 <<G4endl;
528 }
529 */
530 // if(GetProcessName() != "LElastic") std::cout << "Debug -1 "<<aR->GetStatusChange()<<std::endl;
531 if(aR->GetStatusChange()==stopAndKill)
532 {
533 if( xBiasOn && G4UniformRand()<XBiasSurvivalProbability() )
534 {
535 theTotalResult->ProposeParentWeight( XBiasSurvivalProbability()*aT.GetWeight() );
536 }
537 else
538 {
539 theTotalResult->ProposeTrackStatus(fStopAndKill);
540 theTotalResult->ProposeEnergy( 0.0 );
541 }
542 }
543 else if(aR->GetStatusChange()!=stopAndKill )
544 {
545 if(aR->GetStatusChange()==suspend)
546 {
547 theTotalResult->ProposeTrackStatus(fSuspend);
548 if(xBiasOn)
549 {
550 G4Exception("G4HadronicProcess", "007", FatalException,
551 "Cannot cross-section bias a process that suspends tracks.");
552 }
553 } else if (aT.GetKineticEnergy() == 0) {
554 theTotalResult->ProposeTrackStatus(fStopButAlive);
555 }
556
557 if(xBiasOn && G4UniformRand()<XBiasSurvivalProbability())
558 {
559 theTotalResult->ProposeParentWeight( XBiasSurvivalProbability()*aT.GetWeight() );
560 G4double newWeight = aR->GetWeightChange()*aT.GetWeight();
561 if(go_wild(aR->GetEnergyChange()))
562 {
563 G4Exception("G4HadronicProcess", "007", FatalException,
564 "surviving track received NaN energy.");
565 }
566 if(go_wild(aR->GetMomentumChange().x()) ||
567 go_wild(aR->GetMomentumChange().y()) ||
568 go_wild(aR->GetMomentumChange().z()))
569 {
570 G4Exception("G4HadronicProcess", "007", FatalException,
571 "surviving track received NaN momentum.");
572 }
573 G4double newM=aT.GetDefinition()->GetPDGMass();
574 G4double newE=aR->GetEnergyChange() + newM;
575 G4double newP=std::sqrt(newE*newE - newM*newM);
576 G4DynamicParticle * aNew =
577 new G4DynamicParticle(aT.GetDefinition(), newE, newP*aR->GetMomentumChange());
578 G4HadSecondary * theSec = new G4HadSecondary(aNew, newWeight);
579 aR->AddSecondary(theSec);
580 }
581 else
582 {
583 G4double newWeight = aR->GetWeightChange()*aT.GetWeight();
584 theTotalResult->ProposeParentWeight(newWeight); // This is multiplicative
585 if(aR->GetEnergyChange()>-.5)
586 {
587 if(go_wild(aR->GetEnergyChange()))
588 {
589 G4Exception("G4HadronicProcess", "007", FatalException,
590 "track received NaN energy.");
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 G4cerr << "Track status is "<< aR->GetStatusChange()<<G4endl;
602 G4Exception("G4HadronicProcess", "007", FatalException,
603 "use of unsupported track-status.");
604 }
605
606 if(GetProcessName() != "hElastic" && GetProcessName() != "HadronElastic"
607 && theTotalResult->GetTrackStatus()==fAlive
608 && aR->GetStatusChange()==isAlive
609 )
610 {
611 // Use for debugging: G4double newWeight = theTotalResult->GetParentWeight();
612
613 G4double newKE = std::max(DBL_MIN, aR->GetEnergyChange());
614 G4DynamicParticle* aNew = new G4DynamicParticle(aT.GetDefinition(),
615 aR->GetMomentumChange(),
616 newKE);
617 G4HadSecondary* theSec = new G4HadSecondary(aNew, 1.0);
618 aR->AddSecondary(theSec);
619 aR->SetStatusChange(stopAndKill);
620
621 theTotalResult->ProposeTrackStatus(fStopAndKill);
622 theTotalResult->ProposeEnergy( 0.0 );
623
624 }
625 theTotalResult->ProposeLocalEnergyDeposit(aR->GetLocalEnergyDeposit());
626 theTotalResult->SetNumberOfSecondaries(aR->GetNumberOfSecondaries());
627
628 if(aR->GetStatusChange() != stopAndKill)
629 {
630 G4double newM=aT.GetDefinition()->GetPDGMass();
631 G4double newE=aR->GetEnergyChange() + newM;
632 G4double newP=std::sqrt(newE*newE - newM*newM);
633 G4ThreeVector newPV = newP*aR->GetMomentumChange();
634 G4LorentzVector newP4(newE, newPV);
635 newP4.rotate(rotation, it);
636 newP4*=aR->GetTrafoToLab();
637 theTotalResult->ProposeMomentumDirection(newP4.vect().unit());
638 }
639
640 for(G4int i=0; i<aR->GetNumberOfSecondaries(); i++)
641 {
642 G4LorentzVector theM = aR->GetSecondary(i)->GetParticle()->Get4Momentum();
643 theM.rotate(rotation, it);
644 theM*=aR->GetTrafoToLab();
645
646 if(go_wild(theM.e()))
647 {
648 G4Exception("G4HadronicProcess", "007", FatalException,
649 "secondary track received NaN energy.");
650 }
651 if(go_wild(theM.x()) ||
652 go_wild(theM.y()) ||
653 go_wild(theM.z()))
654 {
655 G4Exception("G4HadronicProcess", "007", FatalException,
656 "secondary track received NaN momentum.");
657 }
658
659 aR->GetSecondary(i)->GetParticle()->Set4Momentum(theM);
660 G4double time = aR->GetSecondary(i)->GetTime();
661 if(time<0) time = aT.GetGlobalTime();
662
663 G4Track* track = new G4Track(aR->GetSecondary(i)->GetParticle(),
664 time,
665 aT.GetPosition());
666
667 G4double newWeight = aT.GetWeight()*aR->GetSecondary(i)->GetWeight();
668 //static G4double pinelcount=0;
669 if(xBiasOn) newWeight *= XBiasSecondaryWeight();
670 // G4cout << "#### ParticleDebug "
671 // <<GetProcessName()<<" "
672 // <<aR->GetSecondary(i)->GetParticle()->GetDefinition()->GetParticleName()<<" "
673 // <<aScaleFactor<<" "
674 // <<XBiasSurvivalProbability()<<" "
675 // <<XBiasSecondaryWeight()<<" "
676 // <<aT.GetWeight()<<" "
677 // <<aR->GetSecondary(i)->GetWeight()<<" "
678 // <<aR->GetSecondary(i)->GetParticle()->Get4Momentum()<<" "
679 // <<G4endl;
680 track->SetWeight(newWeight);
681 G4double trackDeb = track->GetKineticEnergy();
682 if( ( trackDeb<0
683 || (trackDeb>aT.GetKineticEnergy()+1*GeV) ) && getenv("GHADEnergyBalanceDebug") )
684 {
685 G4cout << "Debugging hadronic processes: "<<track->GetKineticEnergy()
686 <<" "<<aT.GetKineticEnergy()
687 <<" "<<GetProcessName()
688 <<" "<<aT.GetDefinition()->GetParticleName()
689 <<G4endl;
690 }
691 track->SetTouchableHandle(aT.GetTouchableHandle());
692 theTotalResult->AddSecondary(track);
693 }
694
695 aR->Clear();
696 return;
697}
698/* end of file */
Note: See TracBrowser for help on using the repository browser.