lcdd/src

diff -N Cerenkov.cc
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ Cerenkov.cc	17 Jun 2013 22:07:33 -0000	1.1
@@ -0,0 +1,150 @@

+///////////////////////////////////////////////////////////////////////
+// Class to calculate Number of photons created by Cerenkov Radiation
+////////////////////////////////////////////////////////////////////////
+//  Hans Wenzel
+//
+//  The algorithm is stripped from the geant 4 class
+//  G4Cerenkov.cc 
+//  but since we are only interested in the number of photons created don't need to 
+//  actually create optical phtons and put them on the stack. 
+//--------------------------------------------------------------------------------------
+#include "Cerenkov.hh"
+#include "G4Step.hh"
+#include "G4ThreeVector.hh"
+#include<iomanip>
+
+Cerenkov::Cerenkov() {
+    Initialize();
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+Cerenkov::~Cerenkov() {
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void Cerenkov::Initialize() {
+    //
+    // now get the Cerenkov Angle Integrals for all materials that have the refraction index defined
+    //
+    const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
+    G4int numOfMaterials = G4Material::GetNumberOfMaterials();
+    for (G4int i = 0; i < numOfMaterials; i++) {
+        // Retrieve vector of refraction indices for the material
+        // from the material's optical properties table
+        G4Material* aMaterial = (*theMaterialTable)[i];
+        G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
+        G4PhysicsOrderedFreeVector* CerAngleIntegrals = new G4PhysicsOrderedFreeVector();
+        if (aMaterialPropertiesTable) {
+            G4MaterialPropertyVector* theRefractionIndexVector = aMaterialPropertiesTable->GetProperty("RINDEX");
+            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;
+                    CerAngleIntegrals-> 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;
+                        CerAngleIntegrals ->InsertValues(currentPM, currentCAI);
+                        prevPM = currentPM;
+                        prevCAI = currentCAI;
+                        prevRI = currentRI;
+                    }
+                }
+                G4cout << "Material:  " << aMaterial->GetName() << G4endl;
+                G4cout << "Refraction Index: " << G4endl;
+                G4cout << "=================" << G4endl;
+                theRefractionIndexVector->DumpValues();
+                G4cout << "Cerenkov angle Integrals: " << G4endl;
+                G4cout << "=========================" << G4endl;
+                CerAngleIntegrals ->DumpValues();
+                CAI.push_back(aMaterial->GetName());
+                CerenkovAngleIntegrals.push_back(CerAngleIntegrals);
+                RefractionIndeces.push_back(theRefractionIndexVector);
+            }
+        }
+    }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// This routine computes the number of Cerenkov photons produced per
+// GEANT-unit (millimeter) in the current medium.
+//             ^^^^^^^^^^
+
+G4double Cerenkov::GetAverageNumberOfPhotons(const G4double charge, const G4double beta, const G4String Material) const {
+    G4bool Ceren = false;
+    G4int MaterialIndex = -1;
+    //
+    // check if optical properties are defined for this material:
+    //
+    for (unsigned int ii = 0; ii < CAI.size(); ii++) {
+        if (CAI[ii] == Material) {
+            MaterialIndex = ii;
+            Ceren = true;
+            break;
+        }
+    }
+    if (!Ceren)return 0.0;
+    const G4double Rfact = 369.81 / (eV * cm);
+    if (beta <= 0.0)return 0.0;
+    if (abs(charge) == 0.0) return 0.0;
+
+    G4double BetaInverse = 1. / beta;
+    // Min and Max photon energies
+    G4double Pmin = RefractionIndeces[MaterialIndex]->GetMinLowEdgeEnergy();
+    G4double Pmax = RefractionIndeces[MaterialIndex]->GetMaxLowEdgeEnergy();
+
+    // Min and Max Refraction Indices
+    G4double nMin = RefractionIndeces[MaterialIndex]->GetMinValue();
+    G4double nMax = RefractionIndeces[MaterialIndex]->GetMaxValue();
+
+    // Max Cerenkov Angle Integral
+    G4double CAImax = CerenkovAngleIntegrals[MaterialIndex]->GetMaxValue();
+    G4double dp, ge;
+
+    // If n(Pmax) < 1/Beta -- no photons generated
+
+    if (nMax < BetaInverse) {
+        dp = 0;
+        ge = 0;
+    }// otherwise if n(Pmin) >= 1/Beta -- photons generated
+
+    else if (nMin > BetaInverse) {
+        dp = Pmax - Pmin;
+        ge = CAImax;
+    }// If n(Pmin) < 1/Beta, and n(Pmax) >= 1/Beta, then
+        // we need to find a P such that the value of n(P) == 1/Beta.
+        // Interpolation is performed by the GetEnergy() and
+        // Value() methods of the G4MaterialPropertiesTable and
+        // the GetValue() method of G4PhysicsVector.
+
+    else {
+        Pmin = RefractionIndeces[MaterialIndex]->GetEnergy(BetaInverse);
+        dp = Pmax - Pmin;
+
+        // need boolean for current implementation of G4PhysicsVector
+        // ==> being phased out
+        G4bool isOutRange;
+        G4double CAImin = CerenkovAngleIntegrals[MaterialIndex]->GetValue(Pmin, isOutRange);
+        ge = CAImax - CAImin;
+    }
+
+    // Calculate number of photons
+    G4double NumPhotons = Rfact * charge / eplus * charge / eplus *
+            (dp - ge * BetaInverse * BetaInverse);
+    return NumPhotons;
+}
+

lcdd/src

diff -u -r1.7 -r1.8
--- G4OpticalCalorimeter.cc	30 Oct 2009 23:16:59 -0000	1.7
+++ G4OpticalCalorimeter.cc	17 Jun 2013 22:07:33 -0000	1.8
@@ -1,4 +1,4 @@

-// $Header: /cvs/lcd/lcdd/src/G4OpticalCalorimeter.cc,v 1.7 2009/10/30 23:16:59 jeremy Exp $

+// $Header: /cvs/lcd/lcdd/src/G4OpticalCalorimeter.cc,v 1.8 2013/06/17 22:07:33 wenzel Exp $

 
 #include "G4OpticalCalorimeterSD.hh"
 

@@ -6,16 +6,16 @@

 #include "G4OpticalPhoton.hh"
 #include "G4TransportationManager.hh"
 #include "G4VProcess.hh"

-

+#include "G4Poisson.hh"

 // lcdd
 #include "G4Segmentation.hh"

-

+#include "Cerenkov.hh"

 using std::vector;
 
 /**
  * Constructor for the case only one Hit Collection is given.
  * In this case the only the energy deposited by aborbed Cerenkov 

- * hotons is stored. Probably this will be dropped in the 

+ * photons is stored. Probably this will be dropped in the 

  * near future.
  */ 
 G4OpticalCalorimeterSD::G4OpticalCalorimeterSD(G4String sdName,

@@ -27,7 +27,7 @@

                       sdSeg, 
                       compare)
 { 

-

+   CerenGenerator = 0;

 }
 
 /**

@@ -46,7 +46,8 @@

                       hcNames,
                       sdSeg, 
                       compare)

-{     

+{
+  CerenGenerator = 0;    

 }
   
 G4OpticalCalorimeterSD::~G4OpticalCalorimeterSD()

@@ -54,86 +55,77 @@

 
 G4bool G4OpticalCalorimeterSD::ProcessHits(G4Step* aStep, G4TouchableHistory* tahis)
 {

+  if (!CerenGenerator) CerenGenerator= new Cerenkov();
+    G4int NCerenPhotons = 0;
+    G4Track* theTrack = aStep->GetTrack();
+    const G4double charge = theTrack->GetDefinition()->GetPDGCharge();
+    G4StepPoint* pPreStepPoint = aStep->GetPreStepPoint();
+    G4StepPoint* pPostStepPoint = aStep->GetPostStepPoint();
+    G4double beta = 0.5 * (pPreStepPoint ->GetBeta() + pPostStepPoint->GetBeta());
+    const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
+    G4String thematerial = touch->GetVolume()->GetLogicalVolume()->GetMaterial()->GetName();
+    G4double MeanNumberOfPhotons = CerenGenerator->GetAverageNumberOfPhotons(charge, beta, thematerial);
+    if (MeanNumberOfPhotons > 0.0) {
+      G4double step_length = aStep->GetStepLength();
+      MeanNumberOfPhotons = MeanNumberOfPhotons * step_length;
+      NCerenPhotons = (G4int) G4Poisson(MeanNumberOfPhotons);
+    } else {
+      NCerenPhotons = 0;
+    }

-    G4Track * aTrack =  aStep->GetTrack();
-      
-    // check that particle is optical photon:
-    if(aTrack->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition())

+    if (NCerenPhotons <= 0)

       {
 	return G4CalorimeterSD::ProcessHits(aStep, tahis);
       }
     else 

-    {
-      G4SensitiveDetector::ProcessHits(aStep, 0);
-
-        G4String processname = aTrack->GetCreatorProcess()->G4VProcess::GetProcessName();
-
-        // deal only with Cerenkov photons
-        if(processname != "Cerenkov") 
-        {
-            aTrack->SetTrackStatus(fStopAndKill);
-            return false;
-        }

+      {
+	G4SensitiveDetector::ProcessHits(aStep, 0);

         G4ThreeVector myPoint        = aStep->GetPreStepPoint()->GetPosition();
         G4StepPoint*  apreStepPoint  = aStep->GetPreStepPoint();
         G4Navigator* theNavigator    = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();

-        G4VPhysicalVolume* myVolume  = theNavigator->LocateGlobalPointAndSetup(myPoint); 
-

+        G4VPhysicalVolume* myVolume  = theNavigator->LocateGlobalPointAndSetup(myPoint); 	

         if ( getVerbose() > 2 ) {
             G4cout << "Physical volume       = " << myVolume->GetName() << G4endl;
             G4cout << "Point of interaction  = " << myPoint<< G4endl;
             G4cout << "sdname " << GetName() << "  hcname " <<collectionName[0]<< G4endl;
         }

-

         // total photon energy

-        G4double theEdep = aTrack->GetTotalEnergy();
-

+	// G4double theEdep = aTrack->GetTotalEnergy();
+	G4double theEdep =double(NCerenPhotons);

         // get global cell pos from seg
         G4ThreeVector globalCellPos = m_segmentation->getGlobalHitPosPreStep(apreStepPoint);

-

         // reset the seg bins
         m_segmentation->resetBins();

-

         // set the seg bins
         m_segmentation->setBins(aStep);

-      

         // create id and pack into 64
         Id64bit id64 = makeId();     

-      

         // find hit by simple lkp of new hit with above info
         G4CalorimeterHit* thisHit = new G4CalorimeterHit(theEdep, globalCellPos);
         thisHit->setId64bit( id64.getId0(), id64.getId1() );
         G4CalorimeterHit* fndHit  = 0;

-

         // hit is not found?
         if (!(fndHit = findHit(thisHit, eCerenkov))) 

-        {   

+	  {   

             // add it to lkp map    
             hits_vector[eCerenkov].push_back(thisHit);

-

             // add to the HC
             m_hitsCollections[eCerenkov]->insert(thisHit);   

-
-        }

+	  }

         // found a hit
         else 

-        {   

+	  {   

             // don't need to insert thisHit, so delete it
             delete thisHit;

-            thisHit = 0;
-	

+            thisHit = 0;	    

             // incr total edep of the hit
             fndHit->incrEdep(theEdep);

-	

             // for setting contrib
             thisHit = fndHit;

-        }
-

+	  }

         // add McpHitContrib to this hit, setting info from step info
         thisHit->addMcpHitContrib( McpHitContrib( aStep ) );     

-        aTrack->SetTrackStatus(fStopAndKill); // don't step photon any further
-

+	//        aTrack->SetTrackStatus(fStopAndKill); // don't step photon any further

         return true;

-    }  // end optical photon treatment 
-

+    }  // end Cerenkov photon treatment 

 }