G4WeightCutOffProcess.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 // ----------------------------------------------------------------------
00030 // GEANT 4 class source file
00031 //
00032 // G4WeightCutOffProcess.cc
00033 //
00034 // ----------------------------------------------------------------------
00035 
00036 #include "G4WeightCutOffProcess.hh"
00037 //#include "G4VScorer.hh"
00038 #include "G4GeometryCellStep.hh"
00039 //#include "G4GCellFinder.hh"
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                       //                      const G4VGCellFinder &aGCellFinder,
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     //    fGCellFinder(aGCellFinder),
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 // SetParallelWorld 
00105 //
00106 //------------------------------------------------------
00107 void G4WeightCutOffProcess::
00108 SetParallelWorld(G4String parallelWorldName)
00109 {
00110 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00111 // Get pointers of the parallel world and its navigator
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 // Get pointer of navigator
00124 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00125   fGhostWorldName = parallelWorld->GetName();
00126   fGhostWorld = parallelWorld;
00127   fGhostNavigator = fTransportationManager->GetNavigator(fGhostWorld);
00128 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00129 }
00130 
00131 //------------------------------------------------------
00132 //
00133 // StartTracking
00134 //
00135 //------------------------------------------------------
00136 void G4WeightCutOffProcess::StartTracking(G4Track* trk)
00137 {
00138 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00139 // Activate navigator and get the navigator ID
00140 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00141 // G4cout << " G4ParallelWorldScoringProcess::StartTracking" << G4endl;
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 // G4cout << "G4ParallelWorldScoringProcess::StartTracking <<<<<<<<<<<<<<<<<< " << G4endl;
00155 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00156 // Let PathFinder initialize
00157 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00158     fPathFinder->PrepareNewTrack(trk->GetPosition(),trk->GetMomentumDirection());
00159 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00160 
00161 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00162 // Setup initial touchables for the first step
00163 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00164     fOldGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
00165     fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
00166     fNewGhostTouchable = fOldGhostTouchable;
00167     fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
00168     
00169     // Initialize 
00170     fGhostSafety = -1.;
00171     fOnBoundary = false;
00172   }
00173 }
00174 
00175 
00176 G4double G4WeightCutOffProcess::
00177 PostStepGetPhysicalInteractionLength(const G4Track &,
00178                                      G4double, G4ForceCondition* condition)
00179 {
00180 //   *condition = Forced;
00181 //   return kInfinity;
00182 
00183 //  *condition = StronglyForced;
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     //xbug?    fOnBoundary = false;
00197     CopyStep(aStep);
00198 
00199     if(fOnBoundary)
00200       {
00201 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00202 // Locate the point and get new touchable
00203 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00204   //??  fPathFinder->Locate(step.GetPostStepPoint()->GetPosition(),
00205   //??                      step.GetPostStepPoint()->GetMomentumDirection());
00206         fNewGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
00207 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00208       }
00209     else
00210       {
00211         // Do I need this ??????????????????????????????????????????????????????????
00212         // fGhostNavigator->LocateGlobalPointWithinVolume(track.GetPosition());
00213         // ?????????????????????????????????????????????????????????????????????????
00214         
00215         // fPathFinder->ReLocate(track.GetPosition());
00216         
00217         // reuse the touchable
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     //  G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(aStep);
00232     //  G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(fGhostStep);
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     //  G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(aStep);
00262     //  G4GeometryCell postCell = fGCellFinder.GetPostGeometryCell(fGhostStep);
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     //  else
00312     //  { fGhostSafety = -1.; }
00313     if (fGhostSafety < 0.) fGhostSafety = 0.0;
00314     
00315     // ------------------------------------------
00316     // Determination of the proposed STEP LENGTH:
00317     // ------------------------------------------
00318     if (currentMinimumStep <= fGhostSafety && currentMinimumStep > 0.)
00319       {
00320         // I have no chance to limit
00321         returnedStep = currentMinimumStep;
00322         fOnBoundary = false;
00323         proposedSafety = fGhostSafety - currentMinimumStep;
00324       }
00325     else // (currentMinimumStep > fGhostSafety: I may limit the Step)
00326       {
00327         G4FieldTrackUpdator::Update(&fFieldTrack,&track);
00328         //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00329         // ComputeStep
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             // Track is not on the boundary
00339             fOnBoundary = false;
00340             fGhostSafety = fGhostNavigator->ComputeSafety(endTrack.GetPosition());
00341           }
00342         else
00343           {
00344             // Track is on the boundary
00345             fOnBoundary = true;
00346             proposedSafety = fGhostSafety;
00347           }
00348         //xbug? proposedSafety = fGhostSafety;
00349         if(eLimited == kUnique || eLimited == kSharedOther) {
00350           *selection = CandidateForSelection;
00351         }else if (eLimited == kSharedTransport) { 
00352           returnedStep *= (1.0 + 1.0e-9);  
00353           // Expand to disable its selection in Step Manager comparison
00354         }
00355         
00356       }
00357 
00358   // ----------------------------------------------
00359   // Returns the fGhostSafety as the proposedSafety
00360   // The SteppingManager will take care of keeping
00361   // the smallest one.
00362   // ----------------------------------------------
00363     return returnedStep;
00364 
00365   } else {
00366     return DBL_MAX;
00367     //not sensible!    return -1.0;
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   // Dummy ParticleChange ie: does nothing
00390   // Expecting G4Transportation to move the track
00391   pParticleChange->Initialize(track);
00392   return pParticleChange; 
00393 
00394   //  return 0;
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 // Set StepStatus for ghost world
00409 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00410   if(fOnBoundary)
00411   { fGhostPostStepPoint->SetStepStatus(fGeomBoundary); }
00412   else if(fGhostPostStepPoint->GetStepStatus()==fGeomBoundary)
00413   { fGhostPostStepPoint->SetStepStatus(fPostStepDoItProc); }
00414 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00415 }

Generated on Mon May 27 17:50:24 2013 for Geant4 by  doxygen 1.4.7