Allow to limit materials in which Cerenkov process will happen to exclude some...

Allow to limit materials in which Cerenkov process will happen to exclude some which do have refractive index

source/processes/electromagnetic/xrays/include/G4Cerenkov.hh

@@ -167,6 +167,12 @@ class G4Cerenkov : public G4VProcess
   G4double GetAverageNumberOfPhotons(const G4double charge, const G4double beta,
                                      const G4Material* aMaterial,
                                      G4MaterialPropertyVector* Rindex) const;
+  inline void AddActiveMaterial( std::string material );
+
+ private:
+
+  void BuildThePhysicsTable();
+
 
  protected:
   G4PhysicsTable* thePhysicsTable;
@@ -179,6 +185,10 @@ class G4Cerenkov : public G4VProcess
   G4bool fStackingFlag;
 
   G4int fNumPhotons;
+
+  std::vector<std::string> fActiveMaterials; /// List of materials in which process should be active. This is for
+                                             /// optimisation to avoid photons we know will not reach sensitive
+                                             /// detector.
 };
 
 inline G4bool G4Cerenkov::GetTrackSecondariesFirst() const
@@ -207,4 +217,8 @@ inline G4PhysicsTable* G4Cerenkov::GetPhysicsTable() const
   return thePhysicsTable;
 }
 
+inline void G4Cerenkov::AddActiveMaterial( std::string material )
+{
+  fActiveMaterials.push_back( material );
+}
 #endif /* G4Cerenkov_h */

source/processes/electromagnetic/xrays/src/G4Cerenkov.cc

@@ -377,7 +377,105 @@ G4VParticleChange* G4Cerenkov::PostStepDoIt(const G4Track& aTrack,
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 void G4Cerenkov::PreparePhysicsTable(const G4ParticleDefinition&)
 {
-  Initialise();
+  if (thePhysicsTable) return;
+
+  const G4MaterialTable* theMaterialTable=
+  G4Material::GetMaterialTable();
+  G4int numOfMaterials = G4Material::GetNumberOfMaterials();
+
+  // create new physics table
+	
+  thePhysicsTable = new G4PhysicsTable(numOfMaterials);
+
+  // loop for materials
+
+  for (G4int i=0 ; i < numOfMaterials; i++) {
+
+      G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector = 0;
+
+      // Retrieve vector of refraction indices for the material
+      // from the material's optical properties table 
+
+      G4Material* aMaterial = (*theMaterialTable)[i];
+
+      // If we got list of active materials then we activate process only
+      // in materials in the list. If list is empty, use all materials with
+      // refractive index defined. This is used for optimisation to avoid
+      // generation of photons in places from where they will not make it
+      // to sensitive detector.
+      // To avoid long nested if statements, if we ignore material, just
+      // insert it here and skip rest of the loop.
+      if ( fActiveMaterials.size() and std::find( fActiveMaterials.begin(), fActiveMaterials.end(),
+                                                  aMaterial->GetName() ) == fActiveMaterials.end() ) {
+        thePhysicsTable->insertAt(i,aPhysicsOrderedFreeVector);
+        continue;
+      }
+
+      G4MaterialPropertiesTable* aMaterialPropertiesTable =
+                                      aMaterial->GetMaterialPropertiesTable();
+
+      if (aMaterialPropertiesTable) {
+         aPhysicsOrderedFreeVector = new G4PhysicsOrderedFreeVector();
+         G4MaterialPropertyVector* theRefractionIndexVector = 
+                              aMaterialPropertiesTable->GetProperty(kRINDEX);
+
+         if (theRefractionIndexVector) {
+
+            // Retrieve the first refraction index in vector
+            // of (photon energy, refraction index) pairs 
+
+            G4double currentRI = (*theRefractionIndexVector)[0];
+
+            if (currentRI > 1.0) {
+
+               // Create first (photon energy, Cerenkov Integral)
+               // pair  
+
+               G4double currentPM = theRefractionIndexVector->Energy(0);
+               G4double currentCAI = 0.0;
+
+               aPhysicsOrderedFreeVector->InsertValues(currentPM , currentCAI);
+
+               // Set previous values to current ones prior to loop
+
+               G4double prevPM  = currentPM;
+               G4double prevCAI = currentCAI;
+               G4double prevRI  = currentRI;
+
+               // loop over all (photon energy, refraction index)
+               // pairs stored for this material  
+
+               for (size_t ii = 1;
+                           ii < theRefractionIndexVector->GetVectorLength();
+                           ++ii) {
+                   currentRI = (*theRefractionIndexVector)[ii];
+                   currentPM = theRefractionIndexVector->Energy(ii);
+
+                   currentCAI = 0.5*(1.0/(prevRI*prevRI) +
+                                     1.0/(currentRI*currentRI));
+
+                   currentCAI = prevCAI + (currentPM - prevPM) * currentCAI;
+
+                   aPhysicsOrderedFreeVector->
+                                         InsertValues(currentPM, currentCAI);
+
+                   prevPM  = currentPM;
+                   prevCAI = currentCAI;
+                   prevRI  = currentRI;
+               }
+
+            }
+         }
+      }
+
+      // The Cerenkov integral for a given material
+      // will be inserted in thePhysicsTable
+      // according to the position of the material in
+      // the material table. 
+
+      thePhysicsTable->insertAt(i,aPhysicsOrderedFreeVector); 
+
+  }
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......