Commit in lcdd/src on MAIN | |||
Cerenkov.cc | +150 | added 1.1 | |
G4OpticalCalorimeter.cc | +37 | -45 | 1.7 -> 1.8 |
+187 | -45 |
implement optimized optical calorimeter
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; +} +
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
}
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1