Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/utils/src/G4EmBiasingManager.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 //
 29 // GEANT4 Class file
 30 //
 31 //
 32 // File name:     G4EmBiasingManager
 33 //
 34 // Author:        Vladimir Ivanchenko 
 35 //
 36 // Creation date: 28.07.2011
 37 //
 38 // Modifications:
 39 //
 40 // 31-05-12 D. Sawkey put back in high energy limit for brem, russian roulette 
 41 // 30-05-12 D. Sawkey  brem split gammas are unique; do weight tests for 
 42 //          brem, russian roulette
 43 // -------------------------------------------------------------------
 44 //
 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 47 
 48 #include "G4EmBiasingManager.hh"
 49 #include "G4SystemOfUnits.hh"
 50 #include "G4PhysicalConstants.hh"
 51 #include "G4MaterialCutsCouple.hh"
 52 #include "G4ProductionCutsTable.hh"
 53 #include "G4ProductionCuts.hh"
 54 #include "G4Region.hh"
 55 #include "G4RegionStore.hh"
 56 #include "G4Track.hh"
 57 #include "G4Electron.hh"
 58 #include "G4Gamma.hh"
 59 #include "G4VEmModel.hh"
 60 #include "G4LossTableManager.hh"
 61 #include "G4ParticleChangeForLoss.hh"
 62 #include "G4ParticleChangeForGamma.hh"
 63 #include "G4EmParameters.hh"
 64 
 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 66 
 67 G4EmBiasingManager::G4EmBiasingManager()
 68   : fDirectionalSplittingTarget(0.0,0.0,0.0)
 69 {
 70   fSafetyMin = 1.e-6*mm;
 71   theElectron = G4Electron::Electron();
 72   theGamma    = G4Gamma::Gamma();
 73 }
 74 
 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 76 
 77 G4EmBiasingManager::~G4EmBiasingManager() = default;
 78 
 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 80 
 81 void G4EmBiasingManager::Initialise(const G4ParticleDefinition& part,
 82                                     const G4String& procName, G4int verbose)
 83 {
 84   //G4cout << "G4EmBiasingManager::Initialise for "
 85   //         << part.GetParticleName()
 86   //         << " and " << procName << G4endl;
 87   const G4ProductionCutsTable* theCoupleTable=
 88     G4ProductionCutsTable::GetProductionCutsTable();
 89   G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
 90 
 91   if(0 < nForcedRegions) { idxForcedCouple.resize(numOfCouples, -1); }
 92   if(0 < nSecBiasedRegions) { idxSecBiasedCouple.resize(numOfCouples, -1); }
 93 
 94   // Deexcitation
 95   for (G4int j=0; j<numOfCouples; ++j) {
 96     const G4MaterialCutsCouple* couple =
 97       theCoupleTable->GetMaterialCutsCouple(j);
 98     const G4ProductionCuts* pcuts = couple->GetProductionCuts();
 99     if(0 <  nForcedRegions) {
100       for(G4int i=0; i<nForcedRegions; ++i) {
101         if(forcedRegions[i]) {
102           if(pcuts == forcedRegions[i]->GetProductionCuts()) { 
103             idxForcedCouple[j] = i;
104             break; 
105           }
106         }
107       }
108     }
109     if(0 < nSecBiasedRegions) { 
110       for(G4int i=0; i<nSecBiasedRegions; ++i) {
111         if(secBiasedRegions[i]) {
112           if(pcuts == secBiasedRegions[i]->GetProductionCuts()) { 
113             idxSecBiasedCouple[j] = i;
114             break; 
115           }
116         }
117       }
118     }
119   }
120 
121   G4EmParameters* param = G4EmParameters::Instance();
122   SetDirectionalSplitting(param->GetDirectionalSplitting());
123   if (fDirectionalSplitting) {
124     SetDirectionalSplittingTarget(param->GetDirectionalSplittingTarget());
125     SetDirectionalSplittingRadius(param->GetDirectionalSplittingRadius());
126   }
127 
128   if (nForcedRegions > 0 && 0 < verbose) {
129     G4cout << " Forced Interaction is activated for "
130            << part.GetParticleName() << " and " 
131            << procName 
132            << " inside G4Regions: " << G4endl;
133     for (G4int i=0; i<nForcedRegions; ++i) {
134       const G4Region* r = forcedRegions[i];
135       if(r) { G4cout << "           " << r->GetName() << G4endl; }
136     }
137   }
138   if (nSecBiasedRegions > 0 && 0 < verbose) {
139     G4cout << " Secondary biasing is activated for " 
140            << part.GetParticleName() << " and " 
141            << procName 
142            << " inside G4Regions: " << G4endl;
143     for (G4int i=0; i<nSecBiasedRegions; ++i) {
144       const G4Region* r = secBiasedRegions[i];
145       if(r) { 
146         G4cout << "           " << r->GetName() 
147                << "  BiasingWeight= " << secBiasedWeight[i] << G4endl; 
148       }
149     }
150     if (fDirectionalSplitting) {
151       G4cout << "     Directional splitting activated, with target position: "
152              << fDirectionalSplittingTarget/cm
153              << " cm; radius: "
154              << fDirectionalSplittingRadius/cm
155              << "cm." << G4endl;
156     }
157   }
158 }
159 
160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
161 
162 void G4EmBiasingManager::ActivateForcedInteraction(G4double val, 
163                                                    const G4String& rname)
164 {
165   G4RegionStore* regionStore = G4RegionStore::GetInstance();
166   G4String name = rname;
167   if(name == "" || name == "world" || name == "World") {
168     name = "DefaultRegionForTheWorld";
169   }
170   const G4Region* reg = regionStore->GetRegion(name, false);
171   if(!reg) { 
172     G4cout << "### G4EmBiasingManager::ForcedInteraction WARNING: "
173            << " G4Region <"
174            << rname << "> is unknown" << G4endl;
175     return; 
176   }
177 
178   // the region is in the list
179   if (0 < nForcedRegions) {
180     for (G4int i=0; i<nForcedRegions; ++i) {
181       if (reg == forcedRegions[i]) {
182         lengthForRegion[i] = val; 
183         return;
184       }
185     }
186   }
187   if(val < 0.0) { 
188     G4cout << "### G4EmBiasingManager::ForcedInteraction WARNING: "
189            << val << " < 0.0, so no activation for the G4Region <"
190            << rname << ">" << G4endl;
191     return; 
192   }
193 
194   // new region 
195   forcedRegions.push_back(reg);
196   lengthForRegion.push_back(val);
197   ++nForcedRegions;
198 }
199 
200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
201 
202 void 
203 G4EmBiasingManager::ActivateSecondaryBiasing(const G4String& rname, 
204                                              G4double factor,
205                                              G4double energyLimit)
206 {
207   //G4cout << "G4EmBiasingManager::ActivateSecondaryBiasing: "
208   //         << rname << " F= " << factor << " E(MeV)= " << energyLimit/MeV
209   //         << G4endl; 
210   G4RegionStore* regionStore = G4RegionStore::GetInstance();
211   G4String name = rname;
212   if(name == "" || name == "world" || name == "World") {
213     name = "DefaultRegionForTheWorld";
214   }
215   const G4Region* reg = regionStore->GetRegion(name, false);
216   if(!reg) { 
217     G4cout << "### G4EmBiasingManager::ActivateBremsstrahlungSplitting "
218            << "WARNING: G4Region <"
219            << rname << "> is unknown" << G4endl;
220     return; 
221   }
222 
223   // Range cut
224   G4int nsplit = 0;
225   G4double w = factor;
226 
227   // splitting
228   if(factor >= 1.0) {
229     nsplit = G4lrint(factor);
230     w = 1.0/G4double(nsplit);
231 
232     // Russian roulette 
233   } else if(0.0 < factor) { 
234     nsplit = 1;
235     w = 1.0/factor;
236   }
237 
238   // the region is in the list - overwrite parameters
239   if (0 < nSecBiasedRegions) {
240     for (G4int i=0; i<nSecBiasedRegions; ++i) {
241       if (reg == secBiasedRegions[i]) {
242         secBiasedWeight[i] = w;
243         nBremSplitting[i]  = nsplit; 
244         secBiasedEnegryLimit[i] = energyLimit;
245         return;
246       }
247     }
248   }
249   /*
250     G4cout << "### G4EmBiasingManager::ActivateSecondaryBiasing: "
251            << " nsplit= " << nsplit << " for the G4Region <"
252            << rname << ">" << G4endl; 
253   */
254 
255   // new region 
256   secBiasedRegions.push_back(reg);
257   secBiasedWeight.push_back(w);
258   nBremSplitting.push_back(nsplit);
259   secBiasedEnegryLimit.push_back(energyLimit);
260   ++nSecBiasedRegions;
261   //G4cout << "nSecBiasedRegions= " << nSecBiasedRegions << G4endl;
262 }
263 
264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
265 
266 G4double G4EmBiasingManager::GetStepLimit(G4int coupleIdx, 
267                                           G4double previousStep)
268 {
269   if(startTracking) {
270     startTracking = false;
271     G4int i = idxForcedCouple[coupleIdx];
272     if(i < 0) {
273       currentStepLimit = DBL_MAX;
274     } else {
275       currentStepLimit = lengthForRegion[i];
276       if(currentStepLimit > 0.0) { currentStepLimit *= G4UniformRand(); }
277     }
278   } else {
279     currentStepLimit -= previousStep;
280   }
281   if(currentStepLimit < 0.0) { currentStepLimit = 0.0; }
282   return currentStepLimit;
283 }
284 
285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
286 
287 G4double 
288 G4EmBiasingManager::ApplySecondaryBiasing(
289                     std::vector<G4DynamicParticle*>& vd,
290                     const G4Track& track,
291                     G4VEmModel* currentModel,
292                     G4ParticleChangeForLoss* pPartChange,
293                     G4double& eloss,  
294                     G4int coupleIdx,
295                     G4double tcut, 
296                     G4double safety)
297 {
298   G4int index = idxSecBiasedCouple[coupleIdx];
299   G4double weight = 1.;
300   if(0 <= index) {
301     std::size_t n = vd.size();
302 
303     // the check cannot be applied per secondary particle
304     // because weight correction is common, so the first
305     // secondary is checked
306     if((0 < n && vd[0]->GetKineticEnergy() < secBiasedEnegryLimit[index])
307           || fDirectionalSplitting) {
308 
309       G4int nsplit = nBremSplitting[index];
310 
311       // Range cut
312       if(0 == nsplit) { 
313         if(safety > fSafetyMin) { ApplyRangeCut(vd, track, eloss, safety); }
314 
315         // Russian Roulette
316       } else if(1 == nsplit) { 
317         weight = ApplyRussianRoulette(vd, index);
318 
319         // Splitting
320       } else {
321         if (fDirectionalSplitting) {
322           weight = ApplyDirectionalSplitting(vd, track, currentModel, index, tcut);
323         } else {
324           G4double tmpEnergy = pPartChange->GetProposedKineticEnergy();
325           G4ThreeVector tmpMomDir = pPartChange->GetProposedMomentumDirection();
326 
327           weight = ApplySplitting(vd, track, currentModel, index, tcut);
328 
329           pPartChange->SetProposedKineticEnergy(tmpEnergy);
330           pPartChange->ProposeMomentumDirection(tmpMomDir);
331         }
332       }
333     }
334   }
335   return weight;
336 }
337 
338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
339 
340 G4double 
341 G4EmBiasingManager::ApplySecondaryBiasing(
342                   std::vector<G4DynamicParticle*>& vd,
343                   const G4Track& track,
344                   G4VEmModel* currentModel, 
345                   G4ParticleChangeForGamma* pPartChange,
346                   G4double& eloss,  
347                   G4int coupleIdx,
348                   G4double tcut, 
349                   G4double safety)
350 {
351   G4int index = idxSecBiasedCouple[coupleIdx];
352   G4double weight = 1.;
353   if(0 <= index) {
354     std::size_t n = vd.size();
355 
356     // the check cannot be applied per secondary particle
357     // because weight correction is common, so the first
358     // secondary is checked
359     if((0 < n && vd[0]->GetKineticEnergy() < secBiasedEnegryLimit[index])
360           || fDirectionalSplitting) {
361 
362       G4int nsplit = nBremSplitting[index];
363 
364       // Range cut
365       if(0 == nsplit) { 
366         if(safety > fSafetyMin) { ApplyRangeCut(vd, track, eloss, safety); }
367 
368         // Russian Roulette
369       } else if(1 == nsplit) { 
370         weight = ApplyRussianRoulette(vd, index);
371 
372         // Splitting
373       } else {
374         if (fDirectionalSplitting) {
375           weight = ApplyDirectionalSplitting(vd, track, currentModel,
376                                     index, tcut, pPartChange);
377         } else {
378           G4double tmpEnergy = pPartChange->GetProposedKineticEnergy();
379           G4ThreeVector tmpMomDir = pPartChange->GetProposedMomentumDirection();
380 
381           weight = ApplySplitting(vd, track, currentModel, index, tcut);
382 
383           pPartChange->SetProposedKineticEnergy(tmpEnergy);
384           pPartChange->ProposeMomentumDirection(tmpMomDir);
385         }
386       }
387     }
388   }
389   return weight;
390 }
391 
392 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
393 
394 G4double 
395 G4EmBiasingManager::ApplySecondaryBiasing(std::vector<G4Track*>& track,
396                                           G4int coupleIdx)
397 {
398   G4int index = idxSecBiasedCouple[coupleIdx];
399   G4double weight = 1.;
400   if(0 <= index) {
401     std::size_t n = track.size();
402 
403     // the check cannot be applied per secondary particle
404     // because weight correction is common, so the first
405     // secondary is checked
406     if(0 < n && track[0]->GetKineticEnergy() < secBiasedEnegryLimit[index]) {
407 
408       G4int nsplit = nBremSplitting[index];
409 
410         // Russian Roulette only
411       if(1 == nsplit) { 
412         weight = secBiasedWeight[index];
413         for(std::size_t k=0; k<n; ++k) {
414           if(G4UniformRand()*weight > 1.0) {
415             const G4Track* t = track[k];
416             delete t;
417             track[k] = nullptr;
418           }
419         }
420       }
421     }
422   }
423   return weight;
424 }
425 
426 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
427 
428 void
429 G4EmBiasingManager::ApplyRangeCut(std::vector<G4DynamicParticle*>& vd,
430                                   const G4Track& track,
431                                   G4double& eloss, G4double safety)
432 {
433   std::size_t n = vd.size();
434   if(!eIonisation) { 
435     eIonisation = 
436       G4LossTableManager::Instance()->GetEnergyLossProcess(theElectron);
437   }
438   if(eIonisation) { 
439     for(std::size_t k=0; k<n; ++k) {
440       const G4DynamicParticle* dp = vd[k];
441       if(dp->GetDefinition() == theElectron) {
442         G4double e = dp->GetKineticEnergy();
443         if(eIonisation->GetRange(e, track.GetMaterialCutsCouple()) < safety) {
444           eloss += e;
445           delete dp;
446           vd[k] = nullptr;
447         }
448       }
449     }
450   }
451 }
452 
453 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
454 
455 G4bool G4EmBiasingManager::CheckDirection(G4ThreeVector pos,
456                                           G4ThreeVector momdir) const
457 {
458   G4ThreeVector delta = fDirectionalSplittingTarget - pos;
459   G4double angle = momdir.angle(delta);
460   G4double dist = delta.cross(momdir).mag();
461   if (dist <= fDirectionalSplittingRadius && angle < halfpi) {
462     return true;
463   }
464   return false;
465 }
466 
467 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
468 
469 G4double
470 G4EmBiasingManager::ApplySplitting(std::vector<G4DynamicParticle*>& vd,
471                                    const G4Track& track,
472                                    G4VEmModel* currentModel, 
473                                    G4int index,
474                                    G4double tcut)
475 {
476   // method is applied only if 1 secondary created PostStep 
477   // in the case of many secondaries there is a contradiction
478   G4double weight = 1.;
479   std::size_t n = vd.size();
480   G4double w = secBiasedWeight[index];
481 
482   if(1 != n || 1.0 <= w) { return weight; }
483 
484   G4double trackWeight = track.GetWeight();
485   const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
486 
487   G4int nsplit = nBremSplitting[index];
488 
489   // double splitting is suppressed 
490   if(1 < nsplit && trackWeight>w) {
491 
492     weight = w;
493     if(nsplit > (G4int)tmpSecondaries.size()) { 
494       tmpSecondaries.reserve(nsplit);
495     }
496     const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
497     // start from 1, because already one secondary created
498     for(G4int k=1; k<nsplit; ++k) {
499       tmpSecondaries.clear();
500       currentModel->SampleSecondaries(&tmpSecondaries, couple, dynParticle, 
501                                       tcut);
502       for (std::size_t kk=0; kk<tmpSecondaries.size(); ++kk) {
503         vd.push_back(tmpSecondaries[kk]);
504       }
505     }
506   }
507   return weight;
508 }
509 
510 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
511 
512 G4double
513 G4EmBiasingManager::ApplyDirectionalSplitting(
514                                    std::vector<G4DynamicParticle*>& vd,
515                                    const G4Track& track,
516                                    G4VEmModel* currentModel,
517                                    G4int index,
518                                    G4double tcut,
519                                    G4ParticleChangeForGamma* partChange)
520 {
521   // primary is gamma. do splitting/RR as appropriate
522   // method applied for any number of secondaries
523 
524   G4double weight = 1.0;
525   G4double w = secBiasedWeight[index];
526 
527   fDirectionalSplittingWeights.clear();
528   if(1.0 <= w) {
529     fDirectionalSplittingWeights.push_back(weight);
530     return weight;
531   }
532 
533   G4double trackWeight = track.GetWeight();
534   G4int nsplit = nBremSplitting[index];
535 
536   // double splitting is suppressed
537   if(1 < nsplit && trackWeight>w) {
538 
539     weight = w;
540     const G4ThreeVector pos = track.GetPosition();
541 
542     G4bool foundPrimaryParticle = false;
543     G4double primaryEnergy = 0.;
544     G4ThreeVector primaryMomdir(0.,0.,0.);
545     G4double primaryWeight = trackWeight;
546 
547     tmpSecondaries = vd;
548     vd.clear();
549     vd.reserve(nsplit);
550     for (G4int k=0; k<nsplit; ++k) {
551       if (k>0) {  // for k==0, SampleSecondaries has already been called
552         tmpSecondaries.clear();
553         // SampleSecondaries modifies primary info stored in partChange
554         currentModel->SampleSecondaries(&tmpSecondaries,
555                                         track.GetMaterialCutsCouple(),
556                                         track.GetDynamicParticle(), tcut);
557       }
558       for (std::size_t kk=0; kk<tmpSecondaries.size(); ++kk) {
559         if (tmpSecondaries[kk]->GetParticleDefinition() == theGamma) {
560           if (CheckDirection(pos, tmpSecondaries[kk]->GetMomentumDirection())){
561             vd.push_back(tmpSecondaries[kk]);
562             fDirectionalSplittingWeights.push_back(1.);
563           } else if (G4UniformRand() < w) {
564             vd.push_back(tmpSecondaries[kk]);
565             fDirectionalSplittingWeights.push_back(1./weight);
566           } else {
567             delete tmpSecondaries[kk];
568             tmpSecondaries[kk] = nullptr;
569           }
570         } else if (k==0) { // keep charged 2ry from first splitting
571           vd.push_back(tmpSecondaries[kk]);
572           fDirectionalSplittingWeights.push_back(1./weight);
573         } else {
574           delete tmpSecondaries[kk];
575           tmpSecondaries[kk] = nullptr;
576         }
577       }
578 
579       // primary
580       G4double en = partChange->GetProposedKineticEnergy();
581       if (en>0.) { // don't add if kinetic energy = 0
582         G4ThreeVector momdir = partChange->GetProposedMomentumDirection();
583         if (CheckDirection(pos,momdir)) {
584           // keep only one primary; others are secondaries
585           if (!foundPrimaryParticle) {
586             primaryEnergy = en;
587             primaryMomdir = momdir;
588             foundPrimaryParticle = true;
589             primaryWeight = weight;
590           } else {
591             auto dp = new G4DynamicParticle(theGamma,
592                           partChange->GetProposedMomentumDirection(),
593                           partChange->GetProposedKineticEnergy());
594             vd.push_back(dp);
595             fDirectionalSplittingWeights.push_back(1.);
596           }
597         } else if (G4UniformRand()<w) { // not going to target. play RR.
598           if (!foundPrimaryParticle) {
599             foundPrimaryParticle = true;
600             primaryEnergy = en;
601             primaryMomdir = momdir;
602             primaryWeight = 1.;
603           } else {
604             auto dp = new G4DynamicParticle(theGamma,
605                           partChange->GetProposedMomentumDirection(),
606                           partChange->GetProposedKineticEnergy());
607             vd.push_back(dp);
608             fDirectionalSplittingWeights.push_back(1./weight);
609           }
610         }
611       }
612     }  // end of loop over nsplit
613 
614     partChange->ProposeWeight(primaryWeight);
615     partChange->SetProposedKineticEnergy(primaryEnergy);
616     partChange->ProposeMomentumDirection(primaryMomdir);
617   } else {
618     for (std::size_t i = 0; i < vd.size(); ++i) {
619       fDirectionalSplittingWeights.push_back(1.);
620     }
621   }
622 
623   return weight;
624 }
625 
626 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
627 
628 G4double G4EmBiasingManager::GetWeight(G4int i)
629 {
630   // normally return 1. If a directionally split particle survives RR,
631   //  return 1./(splitting factor)
632   if (fDirectionalSplittingWeights.size() >= (unsigned int)(i+1) ) {
633     G4double w = fDirectionalSplittingWeights[i];
634     fDirectionalSplittingWeights[i] = 1.; // ensure it's not used again
635     return w;
636   } else {
637     return 1.;
638   }
639 }
640 
641 G4double
642 G4EmBiasingManager::ApplyDirectionalSplitting(
643                                   std::vector<G4DynamicParticle*>& vd,
644                                   const G4Track& track,
645                                   G4VEmModel* currentModel,
646                                   G4int index,
647                                   G4double tcut)
648 {
649   // primary is not a gamma. Do nothing with primary
650 
651   G4double weight = 1.0;
652   G4double w = secBiasedWeight[index];
653 
654   fDirectionalSplittingWeights.clear();
655   if(1.0 <= w) {
656     fDirectionalSplittingWeights.push_back(weight);
657     return weight;
658   }
659 
660   G4double trackWeight = track.GetWeight();
661   G4int nsplit = nBremSplitting[index];
662 
663   // double splitting is suppressed
664   if(1 < nsplit && trackWeight>w) {
665 
666     weight = w;
667     const G4ThreeVector pos = track.GetPosition();
668 
669     tmpSecondaries = vd;
670     vd.clear();
671     vd.reserve(nsplit);
672     for (G4int k=0; k<nsplit; ++k) {
673       if (k>0) {
674         tmpSecondaries.clear();
675         currentModel->SampleSecondaries(&tmpSecondaries,
676                                         track.GetMaterialCutsCouple(),
677                                         track.GetDynamicParticle(), tcut);
678       }
679       //for (auto sec : tmpSecondaries) {
680       for (std::size_t kk=0; kk < tmpSecondaries.size(); ++kk) {
681         if (CheckDirection(pos, tmpSecondaries[kk]->GetMomentumDirection())) {
682           vd.push_back(tmpSecondaries[kk]);
683           fDirectionalSplittingWeights.push_back(1.);
684         } else if (G4UniformRand()<w) {
685           vd.push_back(tmpSecondaries[kk]);
686           fDirectionalSplittingWeights.push_back(1./weight);
687         } else {
688           delete tmpSecondaries[kk];
689           tmpSecondaries[kk] = nullptr;
690         }
691       }
692     }  // end of loop over nsplit
693   } else { // no splitting was done; still need weights
694     for (std::size_t i = 0; i < vd.size(); ++i) {
695       fDirectionalSplittingWeights.push_back(1.0);
696     }
697   }
698   return weight;
699 }
700