00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include "G4WeightCutOffProcess.hh"
00037
00038 #include "G4GeometryCellStep.hh"
00039
00040 #include "G4TouchableHandle.hh"
00041 #include "G4VIStore.hh"
00042
00043 #include "G4Step.hh"
00044 #include "G4Navigator.hh"
00045 #include "G4VTouchable.hh"
00046 #include "G4VPhysicalVolume.hh"
00047 #include "G4ParticleChange.hh"
00048 #include "G4PathFinder.hh"
00049 #include "G4TransportationManager.hh"
00050 #include "G4StepPoint.hh"
00051 #include "G4FieldTrackUpdator.hh"
00052
00053
00054 G4WeightCutOffProcess::
00055 G4WeightCutOffProcess(G4double wsurvival,
00056 G4double wlimit,
00057 G4double isource,
00058 G4VIStore *istore,
00059
00060 const G4String &aName, G4bool para)
00061 : G4VProcess(aName),
00062 fParticleChange(new G4ParticleChange),
00063 fWeightSurvival(wsurvival),
00064 fWeightLimit(wlimit),
00065 fSourceImportance(isource),
00066 fIStore(istore),
00067
00068 fGhostNavigator(0), fNavigatorID(-1), fFieldTrack('0')
00069 {
00070 if (!fParticleChange)
00071 {
00072 G4Exception("G4WeightCutOffProcess::G4WeightCutOffProcess()",
00073 "FatalError", FatalException,
00074 "Failed to allocate G4ParticleChange !");
00075 }
00076
00077 G4VProcess::pParticleChange = fParticleChange;
00078
00079 fGhostStep = new G4Step();
00080 fGhostPreStepPoint = fGhostStep->GetPreStepPoint();
00081 fGhostPostStepPoint = fGhostStep->GetPostStepPoint();
00082
00083 fTransportationManager = G4TransportationManager::GetTransportationManager();
00084 fPathFinder = G4PathFinder::GetInstance();
00085
00086 if (verboseLevel>0)
00087 {
00088 G4cout << GetProcessName() << " is created " << G4endl;
00089 }
00090
00091 paraflag = para;
00092
00093 }
00094
00095 G4WeightCutOffProcess::~G4WeightCutOffProcess()
00096 {
00097 delete fParticleChange;
00098 delete fGhostStep;
00099 }
00100
00101
00102
00103
00104
00105
00106
00107 void G4WeightCutOffProcess::
00108 SetParallelWorld(G4String parallelWorldName)
00109 {
00110
00111
00112
00113 fGhostWorldName = parallelWorldName;
00114 fGhostWorld = fTransportationManager->GetParallelWorld(fGhostWorldName);
00115 fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
00116
00117 }
00118
00119 void G4WeightCutOffProcess::
00120 SetParallelWorld(G4VPhysicalVolume* parallelWorld)
00121 {
00122
00123
00124
00125 fGhostWorldName = parallelWorld->GetName();
00126 fGhostWorld = parallelWorld;
00127 fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
00128
00129 }
00130
00131
00132
00133
00134
00135
00136 void G4WeightCutOffProcess::StartTracking(G4Track* trk)
00137 {
00138
00139
00140
00141
00142
00143 if(paraflag) {
00144 if(fGhostNavigator)
00145 { fNavigatorID = fTransportationManager->ActivateNavigator(fGhostNavigator); }
00146 else
00147 {
00148 G4Exception("G4WeightCutOffProcess::StartTracking",
00149 "ProcParaWorld000",FatalException,
00150 "G4WeightCutOffProcess is used for tracking without having a parallel world assigned");
00151 }
00152
00153
00154
00155
00156
00157
00158 fPathFinder->PrepareNewTrack(trk->GetPosition(),trk->GetMomentumDirection());
00159
00160
00161
00162
00163
00164 fOldGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
00165 fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
00166 fNewGhostTouchable = fOldGhostTouchable;
00167 fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
00168
00169
00170 fGhostSafety = -1.;
00171 fOnBoundary = false;
00172 }
00173 }
00174
00175
00176 G4double G4WeightCutOffProcess::
00177 PostStepGetPhysicalInteractionLength(const G4Track &,
00178 G4double, G4ForceCondition* condition)
00179 {
00180
00181
00182
00183
00184 *condition = Forced;
00185 return DBL_MAX;
00186 }
00187
00188 G4VParticleChange *
00189 G4WeightCutOffProcess::PostStepDoIt(const G4Track& aTrack,
00190 const G4Step &aStep)
00191 {
00192 fParticleChange->Initialize(aTrack);
00193
00194 if(paraflag) {
00195 fOldGhostTouchable = fGhostPostStepPoint->GetTouchableHandle();
00196
00197 CopyStep(aStep);
00198
00199 if(fOnBoundary)
00200 {
00201
00202
00203
00204
00205
00206 fNewGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
00207
00208 }
00209 else
00210 {
00211
00212
00213
00214
00215
00216
00217
00218 fNewGhostTouchable = fOldGhostTouchable;
00219 }
00220
00221 fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
00222 fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
00223
00224 }
00225
00226 if(paraflag) {
00227 G4GeometryCell postCell(*(fGhostPostStepPoint->GetPhysicalVolume()),
00228 fGhostPostStepPoint->GetTouchable()->GetReplicaNumber());
00229
00230
00231
00232
00233 G4double R = fSourceImportance;
00234 if (fIStore)
00235 {
00236 G4double i = fIStore->GetImportance(postCell);
00237 if (i>0)
00238 {
00239 R/=i;
00240 }
00241 }
00242 G4double w = aTrack.GetWeight();
00243 if (w<R*fWeightLimit)
00244 {
00245 G4double ws = fWeightSurvival*R;
00246 G4double p = w/(ws);
00247 if (G4UniformRand()<p)
00248 {
00249 fParticleChange->ProposeTrackStatus(fStopAndKill);
00250 }
00251 else
00252 {
00253 fParticleChange->ProposeWeight(ws);
00254 }
00255 }
00256 } else {
00257
00258 G4GeometryCell postCell(*(aStep.GetPostStepPoint()->GetPhysicalVolume()),
00259 aStep.GetPostStepPoint()->GetTouchable()->GetReplicaNumber());
00260
00261
00262
00263 G4double R = fSourceImportance;
00264 if (fIStore)
00265 {
00266 G4double i = fIStore->GetImportance(postCell);
00267 if (i>0)
00268 {
00269 R/=i;
00270 }
00271 }
00272 G4double w = aTrack.GetWeight();
00273 if (w<R*fWeightLimit)
00274 {
00275 G4double ws = fWeightSurvival*R;
00276 G4double p = w/(ws);
00277 if (G4UniformRand()<p)
00278 {
00279 fParticleChange->ProposeTrackStatus(fStopAndKill);
00280 }
00281 else
00282 {
00283 fParticleChange->ProposeWeight(ws);
00284 }
00285 }
00286 }
00287
00288 return fParticleChange;
00289
00290 }
00291
00292 const G4String &G4WeightCutOffProcess::GetName() const
00293 {
00294 return theProcessName;
00295 }
00296
00297 G4double G4WeightCutOffProcess::
00298 AlongStepGetPhysicalInteractionLength(
00299 const G4Track& track, G4double previousStepSize, G4double currentMinimumStep,
00300 G4double& proposedSafety, G4GPILSelection* selection)
00301 {
00302 if(paraflag) {
00303 static G4FieldTrack endTrack('0');
00304 static ELimited eLimited;
00305
00306 *selection = NotCandidateForSelection;
00307 G4double returnedStep = DBL_MAX;
00308
00309 if (previousStepSize > 0.)
00310 { fGhostSafety -= previousStepSize; }
00311
00312
00313 if (fGhostSafety < 0.) fGhostSafety = 0.0;
00314
00315
00316
00317
00318 if (currentMinimumStep <= fGhostSafety && currentMinimumStep > 0.)
00319 {
00320
00321 returnedStep = currentMinimumStep;
00322 fOnBoundary = false;
00323 proposedSafety = fGhostSafety - currentMinimumStep;
00324 }
00325 else
00326 {
00327 G4FieldTrackUpdator::Update(&fFieldTrack,&track);
00328
00329
00330
00331 returnedStep
00332 = fPathFinder->ComputeStep(fFieldTrack,currentMinimumStep,fNavigatorID,
00333 track.GetCurrentStepNumber(),fGhostSafety,eLimited,
00334 endTrack,track.GetVolume());
00335
00336 if(eLimited == kDoNot)
00337 {
00338
00339 fOnBoundary = false;
00340 fGhostSafety = fGhostNavigator->ComputeSafety(endTrack.GetPosition());
00341 }
00342 else
00343 {
00344
00345 fOnBoundary = true;
00346 proposedSafety = fGhostSafety;
00347 }
00348
00349 if(eLimited == kUnique || eLimited == kSharedOther) {
00350 *selection = CandidateForSelection;
00351 }else if (eLimited == kSharedTransport) {
00352 returnedStep *= (1.0 + 1.0e-9);
00353
00354 }
00355
00356 }
00357
00358
00359
00360
00361
00362
00363 return returnedStep;
00364
00365 } else {
00366 return DBL_MAX;
00367
00368 }
00369
00370 }
00371
00372
00373 G4double G4WeightCutOffProcess::
00374 AtRestGetPhysicalInteractionLength(const G4Track& ,
00375 G4ForceCondition*)
00376 {
00377 return -1.0;
00378 }
00379
00380 G4VParticleChange*
00381 G4WeightCutOffProcess::AtRestDoIt(const G4Track&, const G4Step&)
00382 {
00383 return 0;
00384 }
00385
00386 G4VParticleChange*
00387 G4WeightCutOffProcess::AlongStepDoIt(const G4Track& track, const G4Step&)
00388 {
00389
00390
00391 pParticleChange->Initialize(track);
00392 return pParticleChange;
00393
00394
00395 }
00396
00397 void G4WeightCutOffProcess::CopyStep(const G4Step & step)
00398 {
00399 fGhostStep->SetTrack(step.GetTrack());
00400 fGhostStep->SetStepLength(step.GetStepLength());
00401 fGhostStep->SetTotalEnergyDeposit(step.GetTotalEnergyDeposit());
00402 fGhostStep->SetControlFlag(step.GetControlFlag());
00403
00404 *fGhostPreStepPoint = *(step.GetPreStepPoint());
00405 *fGhostPostStepPoint = *(step.GetPostStepPoint());
00406
00407
00408
00409
00410 if(fOnBoundary)
00411 { fGhostPostStepPoint->SetStepStatus(fGeomBoundary); }
00412 else if(fGhostPostStepPoint->GetStepStatus()==fGeomBoundary)
00413 { fGhostPostStepPoint->SetStepStatus(fPostStepDoItProc); }
00414
00415 }