Commit 922c715f authored by Frank Gaede's avatar Frank Gaede
Browse files

change to layered sandwich calorimeter

parent c856d598
......@@ -52,13 +52,15 @@ public:
virtual void ConstructSDandField();
const G4VPhysicalVolume* GetCristal(int aNumCrystal)
{return fCrystalPhys[aNumCrystal];};
const G4VPhysicalVolume* GetLayer(int aNumLayer)
{return fLayerPhys[aNumLayer];};
private:
G4LogicalVolume* fCrystalLog;
G4VPhysicalVolume* fCrystalPhys[100];
G4LogicalVolume* fLayerLogCalo;
G4LogicalVolume* fLayerLogA;
G4LogicalVolume* fLayerLogS;
G4VPhysicalVolume* fLayerPhys[100];
G4Region* fRegion;
static G4ThreadLocal GFlashShowerModel* fFastShowerModel;
......
......@@ -62,7 +62,7 @@ class SimpleCaloHit : public G4VHit
private:
G4double fEdep;
G4ThreeVector fPos;
G4int fCrystalNumber;
G4int fLayerNumber;
G4ThreeVector fStart;
G4RotationMatrix fRot;
const G4LogicalVolume* fLogV;
......@@ -76,10 +76,10 @@ class SimpleCaloHit : public G4VHit
{ return fEdep; };
inline void SetPos(G4ThreeVector xyz)
{ fPos = xyz; };
inline G4int GetCrystalNum()
{ return fCrystalNumber; };
inline void SetCrystalNum(G4int num)
{ fCrystalNumber=num; };
inline G4int GetLayerNum()
{ return fLayerNumber; };
inline void SetLayerNum(G4int num)
{ fLayerNumber=num; };
inline G4ThreeVector GetPos()
{ return fPos; };
inline void SetStart(G4ThreeVector xyz)
......
......@@ -69,7 +69,7 @@ G4ThreadLocal GFlashHitMaker* SimpleCaloDetectorConstruction::fHitMaker = 0;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
SimpleCaloDetectorConstruction::SimpleCaloDetectorConstruction()
:G4VUserDetectorConstruction(), fCrystalLog(nullptr), fCrystalPhys{}, fRegion(nullptr)
:G4VUserDetectorConstruction(), fLayerLogCalo{}, fLayerLogA{}, fLayerLogS{}, fLayerPhys{}, fRegion(nullptr)
{
G4cout<<"SimpleCaloDetectorConstruction::Detector constructor"<<G4endl;
}
......@@ -129,30 +129,27 @@ G4VPhysicalVolume* SimpleCaloDetectorConstruction::Construct()
//------------------------------
// Calorimeter segments
//------------------------------
// Simplified `CMS-like` PbWO4 crystal calorimeter
// Simplified layered sandwich calorimeter
G4int nbOfCrystals = 10; // this are the crystals PER ROW in this example
// cube of 10 x 10 crystals
// don't change it @the moment, since
// the readout in event action assumes this
// dimensions and is not automatically adapted
// in this version of the example :-(
// Simplified `CMS-like` PbWO4 crystal calorimeter
G4double calo_xside = 31*cm;
G4double calo_yside = 31*cm;
G4double calo_zside = 24*cm;
G4int nbOfLayers = 100; // don't change it @the moment, since
// the readout in event action assumes this
// dimensions and is not automatically adapted
// in this version of the example :-(
G4double crystalWidth = 3*cm;
G4double crystalLength = 24*cm;
G4double calo_xside = 100*cm;
G4double calo_yside = 100*cm;
G4double calo_zAbs = 2*mm;
G4double calo_zSens = 0.5*mm;
G4double calo_zside = nbOfLayers * ( calo_zAbs + calo_zSens ) ;
calo_xside = (crystalWidth*nbOfCrystals)+1*cm;
calo_yside = (crystalWidth*nbOfCrystals)+1*cm;
calo_zside = crystalLength;
G4Box* calo_box= new G4Box("CMS calorimeter", // its name
calo_xside/2., // size
// calo envelope
G4Box* calo_box= new G4Box("Simple calorimeter", // its name
calo_xside/2., // size
calo_yside/2.,
calo_zside/2.);
G4LogicalVolume* caloLog
= new G4LogicalVolume(calo_box, // its solid
air, // its material
......@@ -172,52 +169,73 @@ G4VPhysicalVolume* SimpleCaloDetectorConstruction::Construct()
false,
1);
// Crystals
G4VSolid* crystal_box
= new G4Box("Crystal", // its name
crystalWidth/2,
crystalWidth/2,
crystalLength/2);
// size
fCrystalLog
= new G4LogicalVolume(crystal_box, // its solid
pbWO4, // its material
"CrystalLog"); // its name
// Layers
// absorber
G4VSolid* layer_boxA
= new G4Box("LayerA", // its name
calo_xside/2,
calo_yside/2,
calo_zAbs/2);
fLayerLogA
= new G4LogicalVolume(layer_boxA, // its solid
pbWO4, // its material
"LayerLogA"); // its name
// sensitive
G4VSolid* layer_boxS
= new G4Box("LayerS", // its name
calo_xside/2,
calo_yside/2,
calo_zSens/2);
fLayerLogS
= new G4LogicalVolume(layer_boxS, // its solid
pbWO4, // its material
"LayerLogS"); // its name
for (G4int i=0; i<nbOfCrystals; i++)
{
for (G4int j=0; j<nbOfCrystals; j++)
{
G4int n = i*10+j;
G4ThreeVector crystalPos((i*crystalWidth)-135,
(j*crystalWidth)-135,0 );
fCrystalPhys[n]
= new G4PVPlacement(0, // no rotation
crystalPos, // translation
fCrystalLog,
"crystal", // its name
caloLog,
false,
i);
}
}
G4cout << "There are " << nbOfCrystals <<
" crystals per row in the calorimeter, so in total "<<
nbOfCrystals*nbOfCrystals << " crystals" << G4endl;
G4cout << "They have width of " << crystalWidth /cm <<
" cm and a length of " << crystalLength /cm
<<" cm. The Material is "<< pbWO4 << G4endl;
for (G4int j=0; j<nbOfLayers; j++)
{
G4int n = j;
G4ThreeVector layerPosA( 0,0, ( j - nbOfLayers / 2 ) * (calo_zAbs + calo_zSens ) + calo_zAbs/2 ) ;
G4ThreeVector layerPosS( 0,0, ( j - nbOfLayers / 2 ) * (calo_zAbs + calo_zSens ) + calo_zAbs + calo_zSens/2 ) ;
new G4PVPlacement(0, // no rotation
layerPosA, // translation
fLayerLogA,
"layerA", // its name
caloLog,
false,
j);
fLayerPhys[n]
= new G4PVPlacement(0, // no rotation
layerPosS, // translation
fLayerLogS,
"layerS", // its name
caloLog,
false,
j);
}
G4cout << "There are " << nbOfLayers <<
" layers per row in the calorimeter " << G4endl;
G4cout << "They have absorber thickness of " << calo_zAbs /cm <<
" cm and a sensitive thickness " << calo_zSens /cm
<<" cm. The absorber material is "<< pbWO4
<<" The sensitive material is "<< pbWO4
<< G4endl;
experimentalHallLog->SetVisAttributes(G4VisAttributes::GetInvisible());
G4VisAttributes* caloVisAtt = new G4VisAttributes(G4Colour(1.0,1.0,1.0));
G4VisAttributes* crystalVisAtt = new G4VisAttributes(G4Colour(1.0,1.0,0.0));
G4VisAttributes* layerVisAtt = new G4VisAttributes(G4Colour(1.0,1.0,0.0));
caloLog->SetVisAttributes(caloVisAtt);
fCrystalLog->SetVisAttributes(crystalVisAtt);
fLayerLogS->SetVisAttributes(layerVisAtt);
// define the fParameterisation region
fRegion = new G4Region("crystals");
fRegion = new G4Region("layers");
caloLog->SetRegion(fRegion);
fRegion->AddRootLogicalVolume(caloLog);
......@@ -233,11 +251,12 @@ void SimpleCaloDetectorConstruction::ConstructSDandField()
SimpleCaloSensitiveDetector* CaloSD
= new SimpleCaloSensitiveDetector("Calorimeter",this);
SDman->AddNewDetector(CaloSD);
fCrystalLog->SetSensitiveDetector(CaloSD);
fLayerLogS->SetSensitiveDetector(CaloSD);
// fLayerLogCalo->SetSensitiveDetector(CaloSD);
// Get nist material manager
G4NistManager* nistManager = G4NistManager::Instance();
G4Material* pbWO4 = nistManager->FindOrBuildMaterial("G4_PbWO4");
G4Material* pbWO4 = nistManager->FindOrBuildMaterial("G4_PbWO4");
// -- fast simulation models:
// **********************************************
// * Initializing shower modell
......
......@@ -112,10 +112,10 @@ void SimpleCaloEventAction::EndOfEventAction(const G4Event *evt)
G4ThreeVector mom=pparticle->GetMomentum()/pparticle->GetMomentum().mag();
//@@@ SimpleCaloEventAction: Magicnumber
G4double energyincrystal[100];
G4int hitsincrystal[100];
for (int i=0;i<100;i++) energyincrystal[i]=0.;
for (int i=0;i<100;i++) hitsincrystal[i]=0.;
G4double energyinlayer[100];
G4int hitsinlayer[100];
for (int i=0;i<100;i++) energyinlayer[i]=0.;
for (int i=0;i<100;i++) hitsinlayer[i]=0.;
//@@@ SimpleCaloEventAction: Magicnumber
/// For all Hits in sensitive detector
......@@ -125,14 +125,14 @@ void SimpleCaloEventAction::EndOfEventAction(const G4Event *evt)
if (estep >0.0)
{
totE += (*THC)[i]->GetEdep()/GeV;
G4int num=(*THC)[i]->GetCrystalNum();
G4int num=(*THC)[i]->GetLayerNum();
energyincrystal[num]+=(*THC)[i]->GetEdep()/GeV;
hitsincrystal[num]++;
energyinlayer[num]+=(*THC)[i]->GetEdep()/GeV;
hitsinlayer[num]++;
//G4cout << num << G4endl;
// G4cout << " Crystal Nummer " << (*THC)[i]->GetCrystalNum() << G4endl;
// G4cout << (*THC)[i]->GetCrystalNum() /10 <<
// " "<<(*THC)[i]->GetCrystalNum()%10 << G4endl;
// G4cout << " Layer Nummer " << (*THC)[i]->GetLayerNum() << G4endl;
// G4cout << (*THC)[i]->GetLayerNum() /10 <<
// " "<<(*THC)[i]->GetLayerNum()%10 << G4endl;
G4ThreeVector hitpos=(*THC)[i]->GetPos();
G4ThreeVector l (hitpos.x(), hitpos.y(), hitpos.z());
......@@ -144,55 +144,22 @@ void SimpleCaloEventAction::EndOfEventAction(const G4Event *evt)
G4ThreeVector radial = vtx.cross(l);
}
}
G4cout << " energy deposited by later: " << G4endl;
G4double max = 0;
G4int index = 0;
//Find crystal with maximum energy
//Find layer with maximum energy
for (int i=0;i<100;i++)
{
//G4cout << i <<" " << energyincrystal[i] << G4endl;
if (max <energyincrystal[i])
G4cout << " " << i << " " << energyinlayer[i] << G4endl;
if (max <energyinlayer[i])
{
max=energyincrystal[i];
max=energyinlayer[i];
index = i;
}
}
//G4cout << index <<" NMAX " << index << G4endl;
// 3x3 matrix of crystals around the crystal with the maximum energy deposit
G4double e3x3 =
energyincrystal[index]+energyincrystal[index+1]+energyincrystal[index-1]+
energyincrystal[index-10]+energyincrystal[index-9]+energyincrystal[index-11]+
energyincrystal[index+10]+energyincrystal[index+11]+energyincrystal[index+9];
// 5x5 matrix of crystals around the crystal with the maximum energy deposit
G4double e5x5 =
energyincrystal[index]+energyincrystal[index+1]+energyincrystal[index-1]+
energyincrystal[index+2]+energyincrystal[index-2]+
energyincrystal[index-10]+energyincrystal[index-9]+energyincrystal[index-11]+
energyincrystal[index-8]+energyincrystal[index-12]+
energyincrystal[index+10]+energyincrystal[index+11]+energyincrystal[index+9]+
energyincrystal[index+12]+energyincrystal[index+8];
// 3x3 matrix of crystals around the crystal with the maximum energy deposit
G4int num3x3 =
hitsincrystal[index]+hitsincrystal[index+1]+hitsincrystal[index-1]+
hitsincrystal[index-10]+hitsincrystal[index-9]+hitsincrystal[index-11]+
hitsincrystal[index+10]+hitsincrystal[index+11]+hitsincrystal[index+9];
// 5x5 matrix of crystals around the crystal with the maximum energy deposit
G4int num5x5 =
hitsincrystal[index]+hitsincrystal[index+1]+hitsincrystal[index-1]+
hitsincrystal[index+2]+hitsincrystal[index-2]+
hitsincrystal[index-10]+hitsincrystal[index-9]+hitsincrystal[index-11]+
hitsincrystal[index-8]+hitsincrystal[index-12]+
hitsincrystal[index+10]+hitsincrystal[index+11]+hitsincrystal[index+9]+
hitsincrystal[index+12]+hitsincrystal[index+8];
G4cout << " e1 " << energyincrystal[index]
<< " e3x3 " << e3x3<< " GeV e5x5 " <<e5x5 <<G4endl;
G4cout << " max energy in layer " << index << " : " << max << G4endl;
G4cout << " num1 " << hitsincrystal[index]
<< " num3x3 " << num3x3<< " num5x5 " <<num5x5 <<G4endl;
}
G4cout << " Total energy deposited in the calorimeter: " << totE << " (GeV)" << G4endl;
......
......@@ -68,7 +68,7 @@ SimpleCaloHit::SimpleCaloHit(const SimpleCaloHit &right)
fStart =right.fStart;
fRot = right.fRot;
fLogV = right.fLogV;
fCrystalNumber = right.fCrystalNumber;
fLayerNumber = right.fLayerNumber;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......@@ -81,8 +81,8 @@ const SimpleCaloHit & SimpleCaloHit::operator=(const SimpleCaloHit &right)
fPos = right.fPos;
fRot = right.fRot;
fLogV = right.fLogV;
fCrystalNumber = right.fCrystalNumber;
fCrystalNumber = right.fCrystalNumber;
fLayerNumber = right.fLayerNumber;
fLayerNumber = right.fLayerNumber;
return *this;
}
......
......@@ -89,12 +89,12 @@ G4bool SimpleCaloSensitiveDetector::ProcessHits(G4Step* aStep,G4TouchableHistory
fCaloHitsCollection->insert(caloHit);
if (ROhist){;}
G4VPhysicalVolume* physVol = theTouchable->GetVolume();
G4int crystalnum=0;
G4int layernum=0;
for(int i=0;i<100;i++) //@@@@@@@ SimpleCaloSensitiveDetector:vorsichty
{
if(physVol == fDetector->GetCristal(i)) crystalnum= i;
if(physVol == fDetector->GetLayer(i)) layernum= i;
}
caloHit->SetCrystalNum(crystalnum);
caloHit->SetLayerNum(layernum);
std::cout << " ### processhit: - G4 " << caloHit->GetEdep() << " "
<< caloHit->GetPos()[0] << " "
......@@ -121,17 +121,18 @@ G4bool SimpleCaloSensitiveDetector::ProcessHits(G4GFlashSpot*aSpot ,G4TouchableH
fCaloHitsCollection->insert(caloHit);
if (ROhist){;}
//cout <<pCurrentVolume->GetName() << endl;
G4int crystalnum=0;
G4int layernum=0;
for(int i=0;i<100;i++) //@@@@@@@ SimpleCaloSensitiveDetector:vorsichty
{
if(pCurrentVolume == fDetector->GetCristal(i)) crystalnum= i;
if(pCurrentVolume == fDetector->GetLayer(i)) layernum= i;
}
caloHit->SetCrystalNum(crystalnum);
caloHit->SetLayerNum(layernum);
std::cout << " ### processhit: - GFlash " << caloHit->GetEdep() << " "
<< caloHit->GetPos()[0] << " "
<< caloHit->GetPos()[1] << " "
<< caloHit->GetPos()[2] << " "
<< caloHit->GetLayerNum()
<< std::endl ;
return true;
......
......@@ -5,7 +5,7 @@
# Sets some default verbose
/gps/position 0 0 0
/gps/particle e+
/gps/energy 2.0 GeV
/gps/energy 20. GeV
/gps/direction 0. 0. 1.0
#
# Use this open statement to create an OpenGL view:
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment