G4ScoreSplittingProcess.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: G4ScoreSplittingProcess.cc 69885 2013-05-17 07:49:30Z gcosmo $
00028 //
00029 
00030 #include "G4ScoreSplittingProcess.hh"
00031 
00032 #include "G4ios.hh"
00033 #include "G4SystemOfUnits.hh"
00034 #include "G4Step.hh"
00035 #include "G4VTouchable.hh"
00036 #include "G4VPhysicalVolume.hh"
00037 #include "G4ParticleChange.hh"
00038 #include "G4TransportationManager.hh"
00039 #include "G4RegularNavigationHelper.hh"
00040 #include "G4ParticleChange.hh"
00041 #include "G4StepPoint.hh"
00042 
00043 #include "G4SDManager.hh"
00044 #include "G4VSensitiveDetector.hh"
00045 
00046 #include "G4EnergySplitter.hh"
00047 #include "G4TouchableHistory.hh"
00048 
00049 //--------------------------------
00050 // Constructor with name and type:
00051 //--------------------------------
00052 G4ScoreSplittingProcess::
00053 G4ScoreSplittingProcess(const G4String& processName,G4ProcessType theType)
00054   :G4VProcess(processName,theType),
00055    fOldTouchableH(),  fNewTouchableH(), fInitialTouchableH(), fFinalTouchableH()
00056 {
00057   pParticleChange = &xParticleChange;
00058 
00059   fSplitStep = new G4Step();
00060   fSplitPreStepPoint =  fSplitStep->GetPreStepPoint();
00061   fSplitPostStepPoint = fSplitStep->GetPostStepPoint();
00062 
00063   if (verboseLevel>0)
00064   {
00065     G4cout << GetProcessName() << " is created " << G4endl;
00066   }
00067   fpEnergySplitter = new G4EnergySplitter(); 
00068 }
00069 
00070 // -----------
00071 // Destructor:
00072 // -----------
00073 G4ScoreSplittingProcess::~G4ScoreSplittingProcess()
00074 {
00075   delete fSplitStep;
00076   delete fpEnergySplitter; 
00077 }
00078 
00079 //------------------------------------------------------
00080 //
00081 // StartTracking
00082 //
00083 //------------------------------------------------------
00084 void G4ScoreSplittingProcess::StartTracking(G4Track* trk)
00085 {
00086 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00087 // Setup initial touchables for the first step
00088 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00089   const G4Step* pStep= trk->GetStep();
00090 
00091   fOldTouchableH = trk->GetTouchableHandle();
00092   *fSplitPreStepPoint = *(pStep->GetPreStepPoint()); // Best to copy, so as to initialise
00093   fSplitPreStepPoint->SetTouchableHandle(fOldTouchableH);
00094   fNewTouchableH = fOldTouchableH;
00095   *fSplitPostStepPoint= *(pStep->GetPostStepPoint()); // Best to copy, so as to initialise
00096   fSplitPostStepPoint->SetTouchableHandle(fNewTouchableH);
00097 
00099   fSplitPreStepPoint ->SetStepStatus(fUndefined);
00100   fSplitPostStepPoint->SetStepStatus(fUndefined);
00101 }
00102 
00103 
00104 //----------------------------------------------------------
00105 //
00106 //  PostStepGetPhysicalInteractionLength()
00107 //
00108 //----------------------------------------------------------
00109 G4double 
00110 G4ScoreSplittingProcess::PostStepGetPhysicalInteractionLength(
00111          const G4Track& /*track*/, 
00112          G4double   /*previousStepSize*/, 
00113          G4ForceCondition* condition)
00114 {
00115   // This process must be invoked anyway to score the hit
00116   //   - to do the scoring if the current volume is a regular structure, or
00117   //   - else to toggle the flag so that the SteppingManager does the scoring.
00118   *condition = StronglyForced;
00119 
00120   // Future optimisation: check whether in regular structure.  
00121   //  If it is in regular structure,  be StronglyForced
00122   //  If  not  in regular structure, 
00123   //         ask to be called only if SteppingControl is AvoidHitInvocation
00124   //         in order to reset it to NormalCondition
00125 
00126   return DBL_MAX;
00127 }
00128 
00129 //------------------------------------
00130 //
00131 //             PostStepDoIt()
00132 //
00133 //------------------------------------
00134 G4VParticleChange* G4ScoreSplittingProcess::PostStepDoIt(
00135      const G4Track& track,
00136      const G4Step& step)
00137 { 
00138   G4VPhysicalVolume*    pCurrentVolume= track.GetVolume(); 
00139   G4LogicalVolume*      pLogicalVolume= pCurrentVolume->GetLogicalVolume(); 
00140   G4VSensitiveDetector* ptrSD = pLogicalVolume->GetSensitiveDetector();
00141 
00142   pParticleChange->Initialize(track); 
00143   if(  ( ! pCurrentVolume->IsRegularStructure() ) || ( !ptrSD ) 
00144     || G4RegularNavigationHelper::Instance()->GetStepLengths().size() <= 1) {
00145      // Set the flag to make sure that Stepping Manager does the scoring
00146      pParticleChange->ProposeSteppingControl( NormalCondition );     
00147   } else { 
00148      G4ThreeVector preStepPosition, postStepPosition, direction, finalPostStepPosition;
00149      pParticleChange->ProposeSteppingControl( AvoidHitInvocation );
00150  
00151      G4double totalEnergyDeposit=  step.GetTotalEnergyDeposit();
00152      G4StepStatus  fullStepStatus= step.GetPostStepPoint()->GetStepStatus(); 
00153 
00154      CopyStepStart(step);
00155      fSplitPreStepPoint->SetSensitiveDetector(ptrSD);
00156      fOldTouchableH = fInitialTouchableH;
00157      fNewTouchableH=  fOldTouchableH; 
00158      *fSplitPostStepPoint= *(step.GetPreStepPoint()); 
00159      
00160      // Split the energy
00161      // ----------------
00162      G4int numberVoxelsInStep= fpEnergySplitter->SplitEnergyInVolumes( &step );
00163 
00164      preStepPosition= step.GetPreStepPoint()->GetPosition();
00165      finalPostStepPosition= step.GetPostStepPoint()->GetPosition(); 
00166      direction= (finalPostStepPosition - preStepPosition).unit(); 
00167 
00168      fFinalTouchableH= track.GetNextTouchableHandle(); 
00169 
00170      postStepPosition= preStepPosition;
00171      // Loop over the sub-parts of this step
00172      G4int iStep;
00173      for ( iStep=0; iStep < numberVoxelsInStep; iStep++ ){ 
00174         G4int     idVoxel=  -1;  // Voxel ID
00175         G4double  stepLength=0.0, energyLoss= 0.0;
00176 
00177         *fSplitPreStepPoint  = *fSplitPostStepPoint;
00178         fOldTouchableH = fNewTouchableH; 
00179 
00180         preStepPosition= postStepPosition;
00181         fSplitPreStepPoint->SetPosition( preStepPosition );
00182         fSplitPreStepPoint->SetTouchableHandle(fOldTouchableH);
00183         
00184         fpEnergySplitter->GetLengthAndEnergyDeposited( iStep, idVoxel, stepLength, energyLoss);
00185 
00186         // Correct the material, so that the track->GetMaterial gives correct answer
00187         pLogicalVolume->SetMaterial( fpEnergySplitter->GetVoxelMaterial( iStep) );  // idVoxel) ); 
00188 
00189         postStepPosition= preStepPosition + stepLength * direction; 
00190         fSplitPostStepPoint->SetPosition(postStepPosition); 
00191 
00192         // Load the Step with the new values
00193         fSplitStep->SetStepLength(stepLength);
00194         fSplitStep->SetTotalEnergyDeposit(energyLoss);
00195         if( iStep < numberVoxelsInStep -1 ){ 
00196           fSplitStep->GetPostStepPoint()->SetStepStatus( fGeomBoundary );
00197           G4int  nextVoxelId= -1;
00198           fpEnergySplitter->GetVoxelID( iStep+1, nextVoxelId );
00199 
00200           // Create new "next" touchable for each section ??
00201           G4VTouchable* fNewTouchablePtr= 
00202               CreateTouchableForSubStep( nextVoxelId, postStepPosition );
00203           fNewTouchableH= G4TouchableHandle(fNewTouchablePtr); 
00204           fSplitPostStepPoint->SetTouchableHandle( fNewTouchableH );
00205         } else { 
00206           fSplitStep->GetPostStepPoint()->SetStepStatus( fullStepStatus );
00207           fSplitPostStepPoint->SetTouchableHandle( fFinalTouchableH );
00208         }
00209 
00210 
00211         // As first approximation, split the NIEL in the same fractions as the energy deposit
00212         G4double eLossFraction; 
00213         eLossFraction= (totalEnergyDeposit>0.0) ? energyLoss / totalEnergyDeposit : 1.0 ; 
00214         fSplitStep->SetNonIonizingEnergyDeposit(step.GetNonIonizingEnergyDeposit()*eLossFraction);
00215         
00216         fSplitPostStepPoint->SetSensitiveDetector( ptrSD ); 
00217 
00218         // Call the Sensitive Detector
00219         ptrSD->Hit(fSplitStep);
00220 
00221         if (verboseLevel>1) Verbose(step);
00222      }
00223   }
00224 
00225   // This must change the Stepping Control
00226   return pParticleChange;
00227 }
00228 
00229 G4TouchableHistory*
00230 G4ScoreSplittingProcess::CreateTouchableForSubStep( G4int newVoxelNum, G4ThreeVector )
00231 {
00232   // G4cout << " Creating touchable handle for voxel-no " << newVoxelNum << G4endl;
00233 
00234   G4TouchableHistory*  oldTouchableHistory= dynamic_cast<G4TouchableHistory*>(fOldTouchableH());
00235   G4TouchableHistory*  ptrTouchableHistory= G4TransportationManager::GetTransportationManager()->
00236                                             GetNavigatorForTracking()->CreateTouchableHistory(oldTouchableHistory->GetHistory());
00237   
00238   // Change the history
00239   G4NavigationHistory* ptrNavHistory= const_cast<G4NavigationHistory*>(ptrTouchableHistory->GetHistory());
00240   G4VPhysicalVolume*   curPhysicalVol= ptrNavHistory->GetTopVolume();
00241 
00242   EVolume curVolumeType=  ptrNavHistory->GetTopVolumeType();
00243   if( curVolumeType == kParameterised )
00244   { 
00245     ptrNavHistory->BackLevel(); 
00246     // G4VPVParameterised parameterisedPV= pNewMother
00247     G4VPVParameterisation* curParamstn=  curPhysicalVol->GetParameterisation();
00248 
00249     // From G4ParameterisedNavigation::IdentifyAndPlaceSolid() inline method
00250     G4VSolid* sampleSolid = curParamstn->ComputeSolid(newVoxelNum, curPhysicalVol);
00251     sampleSolid->ComputeDimensions(curParamstn, newVoxelNum, curPhysicalVol);
00252     curParamstn->ComputeTransformation(newVoxelNum, curPhysicalVol);
00253 
00254     ptrNavHistory->NewLevel( curPhysicalVol, kParameterised, newVoxelNum );
00255   }
00256   else
00257   {
00258     G4cout << " Current volume type is not Parameterised. " << G4endl; 
00259     G4Exception("G4ScoreSplittingProcess::CreateTouchableForSubStep",
00260           "ErrorRegularParamaterisation", JustWarning,
00261          "Score Splitting Process is used for Regular Structure - but did not find one here.");
00262   }
00263   return ptrTouchableHistory; 
00264 }
00265 
00266 void G4ScoreSplittingProcess::CopyStepStart(const G4Step & step)
00267 {
00268   fSplitStep->SetTrack(step.GetTrack());
00269   fSplitStep->SetStepLength(step.GetStepLength());
00270   fSplitStep->SetTotalEnergyDeposit(step.GetTotalEnergyDeposit());
00271   fSplitStep->SetNonIonizingEnergyDeposit(step.GetNonIonizingEnergyDeposit());
00272   fSplitStep->SetControlFlag(step.GetControlFlag());
00273 
00274   *fSplitPreStepPoint  = *(step.GetPreStepPoint());
00275 
00276   fInitialTouchableH=  (step.GetPreStepPoint()) ->GetTouchableHandle();
00277   fFinalTouchableH =   (step.GetPostStepPoint())->GetTouchableHandle();
00278 }
00279 
00280 void G4ScoreSplittingProcess::Verbose(const G4Step& step) const
00281 {
00282   G4cout << "In mass geometry ------------------------------------------------" << G4endl;
00283   G4cout << " StepLength : " << step.GetStepLength()/mm << "      TotalEnergyDeposit : "
00284          << step.GetTotalEnergyDeposit()/MeV << G4endl;
00285   G4cout << " PreStepPoint : "
00286          << step.GetPreStepPoint()->GetPhysicalVolume()->GetName() << " - ";
00287   if(step.GetPreStepPoint()->GetProcessDefinedStep())
00288   { G4cout << step.GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName(); }
00289   else
00290   { G4cout << "NoProcessAssigned"; }
00291   G4cout << G4endl;
00292   G4cout << "                " << step.GetPreStepPoint()->GetPosition() << G4endl;
00293   G4cout << " PostStepPoint : ";
00294   if(step.GetPostStepPoint()->GetPhysicalVolume()) 
00295   { G4cout << step.GetPostStepPoint()->GetPhysicalVolume()->GetName(); }
00296   else
00297   { G4cout << "OutOfWorld"; }
00298   G4cout << " - ";
00299   if(step.GetPostStepPoint()->GetProcessDefinedStep())
00300   { G4cout << step.GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName(); }
00301   else
00302   { G4cout << "NoProcessAssigned"; }
00303   G4cout << G4endl;
00304   G4cout << "                 " << step.GetPostStepPoint()->GetPosition() << G4endl;
00305 
00306   G4cout << "In ghost geometry ------------------------------------------------" << G4endl;
00307   G4cout << " StepLength : " << fSplitStep->GetStepLength()/mm
00308          << "      TotalEnergyDeposit : "
00309          << fSplitStep->GetTotalEnergyDeposit()/MeV << G4endl;
00310   G4cout << " PreStepPoint : "
00311          << fSplitStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() << " ["
00312            << fSplitStep->GetPreStepPoint()->GetTouchable()->GetReplicaNumber()
00313            << " ]" << " - ";
00314   if(fSplitStep->GetPreStepPoint()->GetProcessDefinedStep())
00315   { G4cout << fSplitStep->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName(); }
00316   else
00317   { G4cout << "NoProcessAssigned"; }
00318   G4cout << G4endl;
00319   G4cout << "                " << fSplitStep->GetPreStepPoint()->GetPosition() << G4endl;
00320   G4cout << " PostStepPoint : ";
00321   if(fSplitStep->GetPostStepPoint()->GetPhysicalVolume()) 
00322   {
00323     G4cout << fSplitStep->GetPostStepPoint()->GetPhysicalVolume()->GetName() << " ["
00324            << fSplitStep->GetPostStepPoint()->GetTouchable()->GetReplicaNumber()
00325            << " ]";
00326   }
00327   else
00328   { G4cout << "OutOfWorld"; }
00329   G4cout << " - ";
00330   if(fSplitStep->GetPostStepPoint()->GetProcessDefinedStep())
00331   { G4cout << fSplitStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName(); }
00332   else
00333   { G4cout << "NoProcessAssigned"; }
00334   G4cout << G4endl;
00335   G4cout << "                 " << fSplitStep->GetPostStepPoint()->GetPosition() << " == "
00336          << fSplitStep->GetTrack()->GetMomentumDirection() 
00337          << G4endl;
00338 
00339 }
00340 
00341 
00342 //----------------------------------------------------------
00343 //
00344 //  AtRestGetPhysicalInteractionLength()
00345 //
00346 //----------------------------------------------------------
00347 G4double 
00348 G4ScoreSplittingProcess::AtRestGetPhysicalInteractionLength(
00349          const G4Track& /*track*/, 
00350          G4ForceCondition* condition)
00351 {
00352   *condition = NotForced;  // Was Forced
00353   return DBL_MAX;
00354 }
00355 
00356 
00357 //---------------------------------------
00358 //  AlongStepGetPhysicalInteractionLength
00359 //---------------------------------------
00360 G4double G4ScoreSplittingProcess::AlongStepGetPhysicalInteractionLength(
00361             const G4Track&   , // track, 
00362             G4double         , // previousStepSize, 
00363             G4double         , // currentMinimumStep,
00364             G4double&        , // proposedSafety, 
00365             G4GPILSelection* selection)
00366 {
00367   *selection = NotCandidateForSelection;
00368   return DBL_MAX;
00369 }
00370 
00371 //------------------------------------
00372 //             AlongStepDoIt()
00373 //------------------------------------
00374 
00375 G4VParticleChange* G4ScoreSplittingProcess::AlongStepDoIt(
00376     const G4Track& track, const G4Step& )
00377 {
00378   // Dummy ParticleChange ie: does nothing
00379   // Expecting G4Transportation to move the track
00380   dummyParticleChange.Initialize(track);
00381   return &dummyParticleChange;
00382 }
00383 
00384 //------------------------------------
00385 //             AtRestDoIt()
00386 //------------------------------------
00387 G4VParticleChange* G4ScoreSplittingProcess::AtRestDoIt(
00388      const G4Track& track,
00389      const G4Step&)
00390 { 
00391   pParticleChange->Initialize(track);
00392   return pParticleChange;
00393 }

Generated on Mon May 27 17:49:48 2013 for Geant4 by  doxygen 1.4.7