ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GenericBiasingPhysics.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4GenericBiasingPhysics.cc
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 // ClassName: G4GenericBiasingPhysics
30 //
31 // Author: M. Verderi (Sept.10.2013)
32 // Modified:
33 // 07/11/2014, M. Verderi : fix bug of PhysicsBias(...) which was not taking
34 // into account the vector of processes passed, but biasing all.
35 //
36 //----------------------------------------------------------------------------
37 //
38 //
39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
41 
43 
44 #include "G4ParticleDefinition.hh"
45 #include "G4ProcessManager.hh"
46 
47 #include "G4BiasingHelper.hh"
50 
51 
52 // factory
54 //
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 
60 : G4VPhysicsConstructor(name),
61  fPhysBiasAllCharged(false), fNonPhysBiasAllCharged(false),
62  fPhysBiasAllChargedISL(false), fNonPhysBiasAllChargedISL(false),
63  fPhysBiasAllNeutral(false), fNonPhysBiasAllNeutral(false),
64  fPhysBiasAllNeutralISL(false), fNonPhysBiasAllNeutralISL(false),
65  fVerbose(false)
66 {;}
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69 
71 {;}
72 
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75 
77 {
78  fBiasedParticles.push_back(particleName);
79  std::vector< G4String > dummy;
80  fBiasedProcesses.push_back(dummy);
81  fBiasAllProcesses.push_back(true);
82 }
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85 
86 void G4GenericBiasingPhysics::PhysicsBias(const G4String& particleName, const std::vector< G4String >& processNames)
87 {
88  fBiasedParticles.push_back(particleName);
89  fBiasedProcesses.push_back(processNames);
90  fBiasAllProcesses.push_back(false);
91 }
92 
93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
94 
96 {
97  fNonPhysBiasedParticles.push_back(particleName);
98 }
99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101 
102 void G4GenericBiasingPhysics::Bias(const G4String& particleName)
103 {
104  PhysicsBias(particleName);
105  NonPhysicsBias(particleName);
106 }
107 
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
109 
110 void G4GenericBiasingPhysics::Bias(const G4String& particleName, const std::vector< G4String >& processNames)
111 {
112  PhysicsBias(particleName, processNames);
113  NonPhysicsBias(particleName);
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 void G4GenericBiasingPhysics::PhysicsBiasAddPDGRange( G4int PDGlow, G4int PDGhigh, G4bool includeAntiParticle )
118 {
119  if ( PDGlow > PDGhigh ) G4cout << " G4GenericBiasingPhysics::PhysicsBiasAddPDGRange(...) : PDGlow > PDGhigh, call ignored." << G4endl;
120  fPhysBiasByPDGRangeLow .push_back( PDGlow );
121  fPhysBiasByPDGRangeHigh.push_back( PDGhigh );
122  if ( includeAntiParticle )
123  {
124  fPhysBiasByPDGRangeLow .push_back( -PDGhigh );
125  fPhysBiasByPDGRangeHigh.push_back( -PDGlow );
126  }
127 }
128 
129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
130 void G4GenericBiasingPhysics::NonPhysicsBiasAddPDGRange( G4int PDGlow, G4int PDGhigh, G4bool includeAntiParticle )
131 {
132  if ( PDGlow > PDGhigh ) G4cout << " G4GenericBiasingPhysics::NonPhysicsBiasAddPDGRange(...) : PDGlow > PDGhigh, call ignored." << G4endl;
133  fNonPhysBiasByPDGRangeLow .push_back( PDGlow );
134  fNonPhysBiasByPDGRangeHigh.push_back( PDGhigh );
135  if ( includeAntiParticle )
136  {
137  fNonPhysBiasByPDGRangeLow .push_back( -PDGhigh );
138  fNonPhysBiasByPDGRangeHigh.push_back( -PDGlow );
139  }
140 }
141 
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144 void G4GenericBiasingPhysics::BiasAddPDGRange( G4int PDGlow, G4int PDGhigh, G4bool includeAntiParticle )
145 {
146  if ( PDGlow > PDGhigh ) G4cout << " G4GenericBiasingPhysics::BiasAddPDGRange(...) : PDGlow > PDGhigh, call ignored." << G4endl;
147  PhysicsBiasAddPDGRange ( PDGlow, PDGhigh, includeAntiParticle );
148  NonPhysicsBiasAddPDGRange( PDGlow, PDGhigh, includeAntiParticle );
149 }
150 
152 {
153  fPhysBiasAllCharged = true;
154  fPhysBiasAllChargedISL = includeShortLived;
155 }
157 {
158  fNonPhysBiasAllCharged = true;
159  fNonPhysBiasAllChargedISL = includeShortLived;
160 }
162 {
163  fPhysBiasAllCharged = true;
164  fNonPhysBiasAllCharged = true;
165  fPhysBiasAllChargedISL = includeShortLived;
166  fNonPhysBiasAllChargedISL = includeShortLived;
167 }
169 {
170  fPhysBiasAllNeutral = true;
171  fPhysBiasAllNeutralISL = includeShortLived;
172 }
174 {
175  fNonPhysBiasAllNeutral = true;
176  fNonPhysBiasAllNeutralISL = includeShortLived;
177 }
179 {
180  fPhysBiasAllNeutral = true;
181  fNonPhysBiasAllNeutral = true;
182  fPhysBiasAllNeutralISL = includeShortLived;
183  fNonPhysBiasAllNeutralISL = includeShortLived;
184 }
185 
186 
187 void G4GenericBiasingPhysics::AddParallelGeometry( const G4String& particleName, const G4String& parallelGeometryName )
188 {
189  // -- add particle, caring of possible duplication:
190  G4bool isKnown = false;
191  for ( G4String knownParticle : fParticlesWithParallelGeometries )
192  {
193  if ( knownParticle == particleName )
194  {
195  isKnown = true;
196  break;
197  }
198  }
199 
200  // -- add the geometry, caring for possible duplication of this geometry, for this particle:
201  if ( !isKnown ) fParticlesWithParallelGeometries.push_back( particleName );
202  std::vector< G4String >& geometries = fParallelGeometriesForParticle[particleName];
203 
204  isKnown = false;
205  for ( G4String knownGeometry : geometries )
206  {
207  if ( knownGeometry == parallelGeometryName )
208  {
209  isKnown = true;
210  break;
211  }
212  }
213  if ( !isKnown ) geometries.push_back( parallelGeometryName );
214 
215 }
216 
217 void G4GenericBiasingPhysics::AddParallelGeometry( const G4String& particleName, const std::vector< G4String >& parallelGeometryNames )
218 {
219  for ( G4String geometry : parallelGeometryNames ) AddParallelGeometry( particleName, geometry );
220 }
221 
222 void G4GenericBiasingPhysics::AddParallelGeometry( G4int PDGlow, G4int PDGhigh, const G4String& parallelGeometryName , G4bool includeAntiParticle )
223 {
224  if ( PDGlow > PDGhigh )
225  {
226  G4cout << "G4GenericBiasingPhysics::AddParallelGeometry( G4int PDGlow, G4int PDGhigh, const G4String& parallelGeometryName , G4bool includeAntiParticle = true ), PDGlow > PDGhigh : call ignored" << G4endl;
227  return;
228  }
229 
230  fPDGlowParallelGeometries .push_back( PDGlow );
231  fPDGhighParallelGeometries.push_back( PDGhigh );
232  G4int rangeIndex = fPDGlowParallelGeometries.size() - 1;
233  fPDGrangeParallelGeometries[rangeIndex].push_back( parallelGeometryName );
234 
235  if ( includeAntiParticle )
236  {
237  fPDGlowParallelGeometries .push_back( -PDGhigh );
238  fPDGhighParallelGeometries.push_back( -PDGlow );
239  rangeIndex = fPDGlowParallelGeometries.size() - 1;
240  fPDGrangeParallelGeometries[rangeIndex].push_back( parallelGeometryName );
241  }
242 
243 }
244 
245 void G4GenericBiasingPhysics::AddParallelGeometry( G4int PDGlow, G4int PDGhigh, const std::vector< G4String >& parallelGeometryNames, G4bool includeAntiParticle )
246 {
247  if ( PDGlow > PDGhigh )
248  {
249  G4cout << "G4GenericBiasingPhysics::AddParallelGeometry( G4int PDGlow, G4int PDGhigh, const std::vector< G4String >& parallelGeometryNames, G4bool includeAntiParticle = true ), PDGlow > PDGhigh : call ignored" << G4endl;
250  return;
251  }
252 
253  for ( G4String geometry : parallelGeometryNames ) AddParallelGeometry( PDGlow, PDGhigh, geometry, includeAntiParticle );
254 }
255 
256 void G4GenericBiasingPhysics::AddParallelGeometryAllCharged( const G4String& parallelGeometryName , G4bool includeShortLived )
257 {
258  G4bool isKnown = false;
259  for ( G4String geometry : fParallelGeometriesForCharged )
260  {
261  if ( geometry == parallelGeometryName )
262  {
263  isKnown = true;
264  break;
265  }
266  }
267  if ( !isKnown )
268  {
269  fParallelGeometriesForCharged .push_back( parallelGeometryName );
270  fAllChargedParallelGeometriesISL.push_back( includeShortLived );
271  }
272 }
273 
274 void G4GenericBiasingPhysics::AddParallelGeometryAllCharged( const std::vector< G4String >& parallelGeometryNames, G4bool includeShortLived )
275 {
276  for ( G4String geometry : parallelGeometryNames ) AddParallelGeometryAllCharged( geometry, includeShortLived );
277 }
278 
279 void G4GenericBiasingPhysics::AddParallelGeometryAllNeutral( const G4String& parallelGeometryName , G4bool includeShortLived )
280 {
281  G4bool isKnown = false;
282  for ( G4String geometry : fParallelGeometriesForNeutral )
283  {
284  if ( geometry == parallelGeometryName )
285  {
286  isKnown = true;
287  break;
288  }
289  }
290  if ( !isKnown )
291  {
292  fParallelGeometriesForNeutral .push_back( parallelGeometryName );
293  fAllNeutralParallelGeometriesISL.push_back( includeShortLived );
294  }
295 }
296 
297 void G4GenericBiasingPhysics::AddParallelGeometryAllNeutral( const std::vector< G4String >& parallelGeometryNames, G4bool includeShortLived )
298 {
299  for ( G4String geometry : parallelGeometryNames ) AddParallelGeometryAllNeutral( geometry, includeShortLived );
300 }
301 
302 
303 
304 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
305 
307 {;}
308 
309 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
310 
312 {
313 
314  // -- bias setup per individual particle name:
316  particleIterator->reset();
317 
318  while( (*particleIterator)() )
319  {
321  G4String particleName = particle->GetParticleName();
322  G4ProcessManager* pmanager = particle->GetProcessManager();
323 
324  // -- include non physics process interface for biasing:
325  if ( std::find(fNonPhysBiasedParticles.begin(),
327  particleName ) != fNonPhysBiasedParticles.end() )
328  {
330  }
331 
332  // -- wrap biased physics processes, all processes or only user selected:
333  std::vector< G4String >::const_iterator particleIt =
334  std::find(fBiasedParticles.begin(),
335  fBiasedParticles.end(),
336  particleName );
337  if ( particleIt == fBiasedParticles.end() ) continue;
338 
339  std::vector < G4String >& biasedProcesses = fBiasedProcesses [ particleIt - fBiasedParticles.begin() ];
340  G4bool biasAll = fBiasAllProcesses[ particleIt - fBiasedParticles.begin() ];
341 
342  if ( biasAll )
343  {
344  G4ProcessVector* vprocess = pmanager->GetProcessList();
345  for (std::size_t ip = 0 ; ip < vprocess->size() ; ++ip)
346  {
347  G4VProcess* process = (*vprocess)[ip];
348  biasedProcesses.push_back( process->GetProcessName() );
349  }
350  }
351 
352  G4bool restartLoop(true);
353  while ( restartLoop )
354  {
355  for (std::size_t ip = 0 ; ip < biasedProcesses.size() ; ++ip)
356  {
357  G4bool activ = G4BiasingHelper::ActivatePhysicsBiasing(pmanager, biasedProcesses[ip] );
358  restartLoop = activ;
359  if ( restartLoop ) break;
360  }
361  }
362 
363  }
364 
365 
366  // -- bias setup per group:
367  particleIterator->reset();
368 
369  while( (*particleIterator)() )
370  {
372  G4String particleName = particle->GetParticleName();
373  G4ProcessManager* pmanager = particle->GetProcessManager();
374 
375  // -- exclude particles invidually specified by name:
376  if ( std::find( fNonPhysBiasedParticles.begin(),
378  particleName ) != fNonPhysBiasedParticles.end() ) continue;
379 
380  if ( std::find( fBiasedParticles.begin(),
381  fBiasedParticles.end(),
382  particleName ) != fBiasedParticles.end() ) continue;
383 
384 
385  G4bool physBias(false), nonPhysBias(false);
386 
387  auto PDG = particle->GetPDGEncoding();
388 
389  // -- include particle if in right PDG range:
390  for ( size_t i = 0 ; i < fPhysBiasByPDGRangeLow.size() ; i++ )
391  if ( ( PDG >= fPhysBiasByPDGRangeLow[i] ) && ( PDG <= fPhysBiasByPDGRangeHigh[i] ) )
392  {
393  physBias = true;
394  break;
395  }
396  for ( size_t i = 0 ; i < fNonPhysBiasByPDGRangeLow.size() ; i++ )
397  if ( ( PDG >= fNonPhysBiasByPDGRangeLow[i] ) && ( PDG <= fNonPhysBiasByPDGRangeHigh[i] ) )
398  {
399  nonPhysBias = true;
400  break;
401  }
402 
403  // -- if particle has not yet any biasing, include it on charge criteria:
404  if ( ( physBias == false ) && ( nonPhysBias == false ) )
405  {
406  if ( std::abs( particle->GetPDGCharge() ) > DBL_MIN )
407  {
408  if ( fPhysBiasAllCharged ) if ( fPhysBiasAllChargedISL || !particle->IsShortLived() ) physBias = true;
409  if ( fNonPhysBiasAllCharged ) if ( fNonPhysBiasAllChargedISL || !particle->IsShortLived() ) nonPhysBias = true;
410  }
411  else
412  {
413  if ( fPhysBiasAllNeutral ) if ( fPhysBiasAllNeutralISL || !particle->IsShortLived() ) physBias = true;
414  if ( fNonPhysBiasAllNeutral ) if ( fNonPhysBiasAllNeutralISL || !particle->IsShortLived() ) nonPhysBias = true;
415  }
416  }
417 
418 
419  if ( nonPhysBias ) G4BiasingHelper::ActivateNonPhysicsBiasing(pmanager);
420 
421  if ( physBias )
422  {
423  std::vector < G4String > biasedProcesses;
424  G4ProcessVector* vprocess = pmanager->GetProcessList();
425  for (std::size_t ip = 0 ; ip < vprocess->size() ; ++ip)
426  {
427  G4VProcess* process = (*vprocess)[ip];
428  biasedProcesses.push_back( process->GetProcessName() );
429  }
430 
431  G4bool restartLoop(true);
432  while ( restartLoop )
433  {
434  for (std::size_t ip = 0 ; ip < biasedProcesses.size() ; ++ip)
435  {
436  G4bool activ = G4BiasingHelper::ActivatePhysicsBiasing(pmanager, biasedProcesses[ip] );
437  restartLoop = activ;
438  if ( restartLoop ) break;
439  }
440  }
441  }
442 
443  }
444 
445 
446 
447  // -- Associate parallel geometries:
449 
450 
451  // -- tells what is done:
452  if ( fVerbose )
453  {
454  // -- print:
455  particleIterator->reset();
456 
457  while( (*particleIterator)() )
458  {
460  G4String particleName = particle->GetParticleName();
461  G4ProcessManager* pmanager = particle->GetProcessManager();
462 
463  G4bool isBiased(false);
464  G4String processNames;
465  G4int icount(0);
466 
467  G4ProcessVector* vprocess = pmanager->GetProcessList();
468  for (std::size_t ip = 0 ; ip < vprocess->size() ; ++ip)
469  {
470  G4VProcess* process = (*vprocess)[ip];
471  G4BiasingProcessInterface* pb = dynamic_cast< G4BiasingProcessInterface* >(process);
472  if ( pb != nullptr )
473  {
474  isBiased = true;
475  if ( icount < 3 )
476  {
477  processNames += pb->GetProcessName();
478  processNames += " ";
479  }
480  else
481  {
482  processNames += "\n ";
483  processNames += pb->GetProcessName();
484  processNames += " ";
485  icount = 0;
486  }
487  icount++;
488  }
489  }
490  if ( isBiased )
491  {
492  if ( particle->IsShortLived() )
493  G4cout << std::setw(14) << particleName << " **** : " << processNames << G4endl;
494  else
495  G4cout << std::setw(18) << particleName << " : " << processNames << G4endl;
496  }
497  }
498  }
499 }
500 
501 
502 
504 {
505 
506  // -- parallel geometries for individual particles:
508  particleIterator->reset();
509 
510  while( (*particleIterator)() )
511  {
513  G4String particleName = particle->GetParticleName();
514  G4ProcessManager* pmanager = particle->GetProcessManager();
515 
516  G4bool requested = false;
517  for ( G4String requestedParticles : fParticlesWithParallelGeometries )
518  {
519  if ( requestedParticles == particleName )
520  {
521  requested = true;
522  break;
523  }
524  }
525  if ( requested )
526  {
527  // -- insert biasing process for handling parallel geometries:
529 
530  // -- attach the requested worlds to this process:
531  std::vector< G4String >& parallelWorlds = fParallelGeometriesForParticle[ particleName ];
532  for ( G4String world : parallelWorlds ) limiter->AddParallelWorld( world );
533  }
534 
535  }
536 
537 
538  // -- parallel geometries for particles in PDG ranges:
539  G4int i = 0; // -- index for PDG range
540  for ( G4int PDGlow : fPDGlowParallelGeometries )
541  {
542  G4int PDGhigh = fPDGhighParallelGeometries[i];
543  auto & geometries = fPDGrangeParallelGeometries[i];
544 
545  particleIterator->reset();
546 
547  while( (*particleIterator)() )
548  {
550  G4int particlePDG = particle->GetPDGEncoding();
551  G4ProcessManager* pmanager = particle->GetProcessManager();
552 
553  if ( ( particlePDG >= PDGlow ) && ( particlePDG <= PDGhigh ) )
554  {
555  // -- ยงยง exclude particles from individual list ?
556  // -- insert biasing process for handling parallel geometries:
558 
559  // -- attached the requested worlds to this process:
560  for ( auto geometry : geometries ) limiter->AddParallelWorld( geometry );
561  }
562  }
563  // -- increment index for next PDG range:
564  i++;
565  }
566 
567 
568  // -- parallel geometries for all neutral / charged particles:
569  particleIterator->reset();
570  G4bool islAllNeutral = false;
571  for(auto isln : fAllNeutralParallelGeometriesISL)
572  { islAllNeutral |= isln; }
573  G4bool islAllCharged = false;
574  for(auto islc : fAllChargedParallelGeometriesISL)
575  { islAllCharged |= islc; }
576 
577  while((*particleIterator)())
578  {
580  G4ProcessManager* pmanager = particle->GetProcessManager();
581  if(particle->GetPDGCharge() == 0.)
582  {
583  // Neutral particle
584  if(particle->IsShortLived() && !islAllNeutral) continue;
586  G4int j = 0;
588  {
589  if(!(particle->IsShortLived()) || fAllNeutralParallelGeometriesISL[j])
590  { limiter->AddParallelWorld(wNameN); }
591  j++;
592  }
593  }
594  else
595  {
596  // charged
597  if(particle->IsShortLived() && !islAllCharged) continue;
599  G4int j = 0;
601  {
602  if(!(particle->IsShortLived()) || fAllChargedParallelGeometriesISL[j])
603  { limiter->AddParallelWorld(wNameC); }
604  j++;
605  }
606  }
607  }
608 
609 }