Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/scoring/src/G4ParallelWorldProcess.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 
 30 #include "G4ParallelWorldProcess.hh"
 31 
 32 #include "G4FieldTrackUpdator.hh"
 33 #include "G4Material.hh"
 34 #include "G4Navigator.hh"
 35 #include "G4ParallelWorldProcessStore.hh"
 36 #include "G4ParticleChange.hh"
 37 #include "G4PathFinder.hh"
 38 #include "G4ProductionCuts.hh"
 39 #include "G4ProductionCutsTable.hh"
 40 #include "G4SDManager.hh"
 41 #include "G4Step.hh"
 42 #include "G4StepPoint.hh"
 43 #include "G4TransportationManager.hh"
 44 #include "G4TransportationProcessType.hh"
 45 #include "G4VPhysicalVolume.hh"
 46 #include "G4VSensitiveDetector.hh"
 47 #include "G4VTouchable.hh"
 48 #include "G4ios.hh"
 49 
 50 G4ThreadLocal G4Step* G4ParallelWorldProcess::fpHyperStep = nullptr;
 51 G4ThreadLocal G4int G4ParallelWorldProcess::nParallelWorlds = 0;
 52 G4ThreadLocal G4int G4ParallelWorldProcess::fNavIDHyp = 0;
 53 
 54 const G4Step* G4ParallelWorldProcess::GetHyperStep()
 55 {
 56   return fpHyperStep;
 57 }
 58 
 59 G4int G4ParallelWorldProcess::GetHypNavigatorID()
 60 {
 61   return fNavIDHyp;
 62 }
 63 
 64 G4ParallelWorldProcess::G4ParallelWorldProcess(const G4String& processName, G4ProcessType theType)
 65   : G4VProcess(processName, theType), fFieldTrack('0')
 66 {
 67   SetProcessSubType(PARALLEL_WORLD_PROCESS);
 68   if (fpHyperStep == nullptr) fpHyperStep = new G4Step();
 69   iParallelWorld = ++nParallelWorlds;
 70 
 71   pParticleChange = &aDummyParticleChange;
 72 
 73   fGhostStep = new G4Step();
 74   fGhostPreStepPoint = fGhostStep->GetPreStepPoint();
 75   fGhostPostStepPoint = fGhostStep->GetPostStepPoint();
 76 
 77   fTransportationManager = G4TransportationManager::GetTransportationManager();
 78   fTransportationManager->GetNavigatorForTracking()->SetPushVerbosity(false);
 79   fPathFinder = G4PathFinder::GetInstance();
 80 
 81   fGhostWorldName = "** NotDefined **";
 82   G4ParallelWorldProcessStore::GetInstance()->SetParallelWorld(this, processName);
 83 
 84   if (verboseLevel > 0) {
 85     G4cout << GetProcessName() << " is created " << G4endl;
 86   }
 87 }
 88 
 89 G4ParallelWorldProcess::~G4ParallelWorldProcess()
 90 {
 91   delete fGhostStep;
 92   nParallelWorlds--;
 93   if (nParallelWorlds == 0) {
 94     delete fpHyperStep;
 95     fpHyperStep = nullptr;
 96   }
 97 }
 98 
 99 void G4ParallelWorldProcess::SetParallelWorld(G4String parallelWorldName)
100 {
101   fGhostWorldName = parallelWorldName;
102   fGhostWorld = fTransportationManager->GetParallelWorld(fGhostWorldName);
103   fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
104   fGhostNavigator->SetPushVerbosity(false);
105 }
106 
107 void G4ParallelWorldProcess::SetParallelWorld(G4VPhysicalVolume* parallelWorld)
108 {
109   fGhostWorldName = parallelWorld->GetName();
110   fGhostWorld = parallelWorld;
111   fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
112   fGhostNavigator->SetPushVerbosity(false);
113 }
114 
115 void G4ParallelWorldProcess::StartTracking(G4Track* trk)
116 {
117   if (fGhostNavigator != nullptr) {
118     fNavigatorID = fTransportationManager->ActivateNavigator(fGhostNavigator);
119   }
120   else {
121     G4Exception(
122       "G4ParallelWorldProcess::StartTracking", "ProcParaWorld000", FatalException,
123       "G4ParallelWorldProcess is used for tracking without having a parallel world assigned");
124   }
125   fPathFinder->PrepareNewTrack(trk->GetPosition(), trk->GetMomentumDirection());
126 
127   fOldGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
128   fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
129   fNewGhostTouchable = fOldGhostTouchable;
130   fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
131 
132   fGhostSafety = -1.;
133   fOnBoundary = false;
134   fGhostPreStepPoint->SetStepStatus(fUndefined);
135   fGhostPostStepPoint->SetStepStatus(fUndefined);
136 
137   *(fpHyperStep->GetPostStepPoint()) = *(trk->GetStep()->GetPostStepPoint());
138   if (layeredMaterialFlag) {
139     G4StepPoint* realWorldPostStepPoint = trk->GetStep()->GetPostStepPoint();
140     SwitchMaterial(realWorldPostStepPoint);
141     G4StepPoint* realWorldPreStepPoint = trk->GetStep()->GetPreStepPoint();
142     SwitchMaterial(realWorldPreStepPoint);
143     G4double velocity = trk->CalculateVelocity();
144     realWorldPostStepPoint->SetVelocity(velocity);
145     realWorldPreStepPoint->SetVelocity(velocity);
146     trk->SetVelocity(velocity);
147   }
148   *(fpHyperStep->GetPreStepPoint()) = *(fpHyperStep->GetPostStepPoint());
149 }
150 
151 G4double G4ParallelWorldProcess::AtRestGetPhysicalInteractionLength(const G4Track& /*track*/,
152                                                                     G4ForceCondition* condition)
153 {
154   //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
155   // At Rest must be registered ONLY for the particle which has other At Rest
156   // process(es).
157   //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
158   *condition = Forced;
159   return DBL_MAX;
160 }
161 
162 G4VParticleChange* G4ParallelWorldProcess::AtRestDoIt(const G4Track& track, const G4Step& step)
163 {
164   //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
165   // At Rest must be registered ONLY for the particle which has other At Rest
166   // process(es).
167   //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
168   fOldGhostTouchable = fGhostPostStepPoint->GetTouchableHandle();
169   G4VSensitiveDetector* aSD = nullptr;
170   if (fOldGhostTouchable->GetVolume() != nullptr) {
171     aSD = fOldGhostTouchable->GetVolume()->GetLogicalVolume()->GetSensitiveDetector();
172   }
173   fOnBoundary = false;
174   if (aSD != nullptr) {
175     CopyStep(step);
176     fGhostPreStepPoint->SetSensitiveDetector(aSD);
177 
178     fNewGhostTouchable = fOldGhostTouchable;
179 
180     fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
181     fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
182     if (fNewGhostTouchable->GetVolume() != nullptr) {
183       fGhostPostStepPoint->SetSensitiveDetector(
184         fNewGhostTouchable->GetVolume()->GetLogicalVolume()->GetSensitiveDetector());
185     }
186     else {
187       fGhostPostStepPoint->SetSensitiveDetector(nullptr);
188     }
189 
190     aSD->Hit(fGhostStep);
191   }
192 
193   pParticleChange->Initialize(track);
194   return pParticleChange;
195 }
196 
197 G4double G4ParallelWorldProcess::PostStepGetPhysicalInteractionLength(const G4Track& /*track*/,
198                                                                       G4double /*previousStepSize*/,
199                                                                       G4ForceCondition* condition)
200 {
201   *condition = StronglyForced;
202   return DBL_MAX;
203 }
204 
205 G4VParticleChange* G4ParallelWorldProcess::PostStepDoIt(const G4Track& track, const G4Step& step)
206 {
207   fOldGhostTouchable = fGhostPostStepPoint->GetTouchableHandle();
208   G4VSensitiveDetector* aSD = nullptr;
209   if (fOldGhostTouchable->GetVolume() != nullptr) {
210     aSD = fOldGhostTouchable->GetVolume()->GetLogicalVolume()->GetSensitiveDetector();
211   }
212   CopyStep(step);
213   fGhostPreStepPoint->SetSensitiveDetector(aSD);
214 
215   if (fOnBoundary) {
216     fNewGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
217   }
218   else {
219     fNewGhostTouchable = fOldGhostTouchable;
220   }
221 
222   fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
223   fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
224 
225   if (fNewGhostTouchable->GetVolume() != nullptr) {
226     fGhostPostStepPoint->SetSensitiveDetector(
227       fNewGhostTouchable->GetVolume()->GetLogicalVolume()->GetSensitiveDetector());
228   }
229   else {
230     fGhostPostStepPoint->SetSensitiveDetector(nullptr);
231   }
232 
233   G4VSensitiveDetector* sd = fGhostPreStepPoint->GetSensitiveDetector();
234   if (sd != nullptr) {
235     sd->Hit(fGhostStep);
236   }
237 
238   pParticleChange->Initialize(track);
239   if (layeredMaterialFlag) {
240     G4StepPoint* realWorldPostStepPoint = ((G4Step*)(track.GetStep()))->GetPostStepPoint();
241     SwitchMaterial(realWorldPostStepPoint);
242   }
243   return pParticleChange;
244 }
245 
246 G4double G4ParallelWorldProcess::AlongStepGetPhysicalInteractionLength(const G4Track& track,
247                                                                        G4double previousStepSize,
248                                                                        G4double currentMinimumStep,
249                                                                        G4double& proposedSafety,
250                                                                        G4GPILSelection* selection)
251 {
252   static G4ThreadLocal G4FieldTrack* endTrack_G4MT_TLS_ = nullptr;
253   if (endTrack_G4MT_TLS_ == nullptr) endTrack_G4MT_TLS_ = new G4FieldTrack('0');
254   G4FieldTrack& endTrack = *endTrack_G4MT_TLS_;
255   // static ELimited eLimited;
256   ELimited eLimited;
257   ELimited eLim = kUndefLimited;
258 
259   *selection = NotCandidateForSelection;
260   G4double returnedStep = DBL_MAX;
261 
262   if (previousStepSize > 0.) {
263     fGhostSafety -= previousStepSize;
264   }
265   if (fGhostSafety < 0.) fGhostSafety = 0.0;
266 
267   if (currentMinimumStep <= fGhostSafety && currentMinimumStep > 0.) {
268     // I have no chance to limit
269     returnedStep = currentMinimumStep;
270     fOnBoundary = false;
271     proposedSafety = fGhostSafety - currentMinimumStep;
272     eLim = kDoNot;
273   }
274   else {
275     G4FieldTrackUpdator::Update(&fFieldTrack, &track);
276 
277 #ifdef G4DEBUG_PARALLEL_WORLD_PROCESS
278     if (verboseLevel > 0) {
279       int localVerb = verboseLevel - 1;
280 
281       if (localVerb == 1) {
282         G4cout << " Pll Wrl  proc::AlongStepGPIL " << this->GetProcessName() << G4endl;
283       }
284       else if (localVerb > 1) {
285         G4cout << "----------------------------------------------" << G4endl;
286         G4cout << " ParallelWorldProcess: field Track set to : " << G4endl;
287         G4cout << "----------------------------------------------" << G4endl;
288         G4cout << fFieldTrack << G4endl;
289         G4cout << "----------------------------------------------" << G4endl;
290       }
291     }
292 #endif
293 
294     returnedStep = fPathFinder->ComputeStep(fFieldTrack, currentMinimumStep, fNavigatorID,
295                                             track.GetCurrentStepNumber(), fGhostSafety, eLimited,
296                                             endTrack, track.GetVolume());
297     if (eLimited == kDoNot) {
298       fOnBoundary = false;
299       fGhostSafety = fGhostNavigator->ComputeSafety(endTrack.GetPosition());
300     }
301     else {
302       fOnBoundary = true;
303       // fGhostSafetyEnd = 0.0;    // At end-point of expected step only
304     }
305     proposedSafety = fGhostSafety;
306     if (eLimited == kUnique || eLimited == kSharedOther) {
307       *selection = CandidateForSelection;
308     }
309     else if (eLimited == kSharedTransport) {
310       returnedStep *= (1.0 + 1.0e-9);
311     }
312     eLim = eLimited;
313   }
314 
315   if (iParallelWorld == nParallelWorlds) fNavIDHyp = 0;
316   if (eLim == kUnique || eLim == kSharedOther) fNavIDHyp = fNavigatorID;
317   return returnedStep;
318 }
319 
320 G4VParticleChange* G4ParallelWorldProcess::AlongStepDoIt(const G4Track& track, const G4Step&)
321 {
322   pParticleChange->Initialize(track);
323   return pParticleChange;
324 }
325 
326 void G4ParallelWorldProcess::CopyStep(const G4Step& step)
327 {
328   G4StepStatus prevStat = fGhostPostStepPoint->GetStepStatus();
329 
330   fGhostStep->SetTrack(step.GetTrack());
331   fGhostStep->SetStepLength(step.GetStepLength());
332   fGhostStep->SetTotalEnergyDeposit(step.GetTotalEnergyDeposit());
333   fGhostStep->SetNonIonizingEnergyDeposit(step.GetNonIonizingEnergyDeposit());
334   fGhostStep->SetControlFlag(step.GetControlFlag());
335   fGhostStep->SetSecondary((const_cast<G4Step&>(step)).GetfSecondary());
336 
337   *fGhostPreStepPoint = *(step.GetPreStepPoint());
338   *fGhostPostStepPoint = *(step.GetPostStepPoint());
339 
340   fGhostPreStepPoint->SetStepStatus(prevStat);
341   if (fOnBoundary) {
342     fGhostPostStepPoint->SetStepStatus(fGeomBoundary);
343   }
344   else if (fGhostPostStepPoint->GetStepStatus() == fGeomBoundary) {
345     fGhostPostStepPoint->SetStepStatus(fPostStepDoItProc);
346   }
347 
348   if (iParallelWorld == 1) {
349     G4StepStatus prevStatHyp = fpHyperStep->GetPostStepPoint()->GetStepStatus();
350 
351     fpHyperStep->SetTrack(step.GetTrack());
352     fpHyperStep->SetStepLength(step.GetStepLength());
353     fpHyperStep->SetTotalEnergyDeposit(step.GetTotalEnergyDeposit());
354     fpHyperStep->SetNonIonizingEnergyDeposit(step.GetNonIonizingEnergyDeposit());
355     fpHyperStep->SetControlFlag(step.GetControlFlag());
356 
357     *(fpHyperStep->GetPreStepPoint()) = *(fpHyperStep->GetPostStepPoint());
358     *(fpHyperStep->GetPostStepPoint()) = *(step.GetPostStepPoint());
359 
360     fpHyperStep->GetPreStepPoint()->SetStepStatus(prevStatHyp);
361   }
362 
363   if (fOnBoundary) {
364     fpHyperStep->GetPostStepPoint()->SetStepStatus(fGeomBoundary);
365   }
366 }
367 
368 void G4ParallelWorldProcess::SwitchMaterial(G4StepPoint* realWorldStepPoint)
369 {
370   if (realWorldStepPoint->GetStepStatus() == fWorldBoundary) return;
371   G4VPhysicalVolume* thePhys = fNewGhostTouchable->GetVolume();
372   if (thePhys != nullptr) {
373     G4Material* ghostMaterial = thePhys->GetLogicalVolume()->GetMaterial();
374     if (ghostMaterial != nullptr) {
375       G4Region* ghostRegion = thePhys->GetLogicalVolume()->GetRegion();
376       G4ProductionCuts* prodCuts = realWorldStepPoint->GetMaterialCutsCouple()->GetProductionCuts();
377       if (ghostRegion != nullptr) {
378         G4ProductionCuts* ghostProdCuts = ghostRegion->GetProductionCuts();
379         if (ghostProdCuts != nullptr) prodCuts = ghostProdCuts;
380       }
381       const G4MaterialCutsCouple* ghostMCCouple =
382         G4ProductionCutsTable::GetProductionCutsTable()->GetMaterialCutsCouple(ghostMaterial,
383                                                                                prodCuts);
384       if (ghostMCCouple != nullptr) {
385         realWorldStepPoint->SetMaterial(ghostMaterial);
386         realWorldStepPoint->SetMaterialCutsCouple(ghostMCCouple);
387         *(fpHyperStep->GetPostStepPoint()) = *(fGhostPostStepPoint);
388         fpHyperStep->GetPostStepPoint()->SetMaterial(ghostMaterial);
389         fpHyperStep->GetPostStepPoint()->SetMaterialCutsCouple(ghostMCCouple);
390       }
391       else {
392         G4cout << "!!! MaterialCutsCouple is not found for " << ghostMaterial->GetName() << "."
393                << G4endl << "    Material in real world ("
394                << realWorldStepPoint->GetMaterial()->GetName() << ") is used." << G4endl;
395       }
396     }
397   }
398 }
399 
400 G4bool G4ParallelWorldProcess::IsAtRestRequired(G4ParticleDefinition* partDef)
401 {
402   G4int pdgCode = partDef->GetPDGEncoding();
403   if (pdgCode == 0) {
404     G4String partName = partDef->GetParticleName();
405     if (partName == "geantino") return false;
406     if (partName == "chargedgeantino") return false;
407   }
408   else {
409     if (pdgCode == 11 || pdgCode == 2212) return false;  // electrons and proton
410     pdgCode = std::abs(pdgCode);
411     if (pdgCode == 22) return false;  // gamma and optical photons
412     if (pdgCode == 12 || pdgCode == 14 || pdgCode == 16) return false;  // all neutronos
413   }
414   return true;
415 }
416