Newer
Older
/*****************************************************************************
* Copyright (C) 2009-2020 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
*****************************************************************************/
/*****************************************************************************
* Original Author: Pierre Morfouace contact address: pierre.morfouace2@cea.fr *
* *
* Creation Date : May 2020 *
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
* This class describe PISTA simulation *
* *
*---------------------------------------------------------------------------*
* Comment: *
* *
*****************************************************************************/
// C++ headers
#include <sstream>
#include <cmath>
#include <limits>
//G4 Geometry object
#include "G4Tubs.hh"
#include "G4Box.hh"
//G4 sensitive
#include "G4SDManager.hh"
#include "G4MultiFunctionalDetector.hh"
//G4 various object
#include "G4Material.hh"
#include "G4Transform3D.hh"
#include "G4PVPlacement.hh"
#include "G4VisAttributes.hh"
#include "G4Colour.hh"
// NPTool header
#include "PISTA.hh"
#include "InteractionScorers.hh"
#include "RootOutput.h"
#include "MaterialManager.hh"
#include "NPSDetectorFactory.hh"
#include "NPOptionManager.h"
#include "NPSHitsMap.hh"
// CLHEP header
#include "CLHEP/Random/RandGauss.h"
using namespace std;
using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
namespace PISTA_NS{
// Energy and time Resolution
const double EnergyThreshold = 0.1*MeV;
const double DE_ResoEnergy = 0.040*MeV ;
const double E_ResoEnergy = 0.018*MeV ;
const double TrapezoidBaseLarge = 72.3*mm;
const double TrapezoidBaseSmall = 41.0*mm;
const double TrapezoidHeight = 57.7*mm;
const double FirstStageNbrOfStrips = 91;
const double SecondStageNbrOfStrips = 57;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// PISTA Specific Method
PISTA::PISTA(){
m_Event = new TPISTAData() ;
m_FirstStageScorer = 0;
m_SecondStageScorer = 0;
m_TrapezoidDetector = 0;
}
PISTA::~PISTA(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PISTA::AddDetector(G4ThreeVector A, G4ThreeVector B, G4ThreeVector C, G4ThreeVector D){
m_DefinitionType.push_back(true);
m_A.push_back(A);
m_B.push_back(B);
m_C.push_back(C);
m_D.push_back(D);
m_R.push_back(0);
m_Theta.push_back(0);
m_Phi.push_back(0);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PISTA::AddDetector(double R, double Theta, double Phi){
m_DefinitionType.push_back(false);
m_R.push_back(R);
m_Theta.push_back(Theta);
m_Phi.push_back(Phi);
G4ThreeVector empty = G4ThreeVector(0,0,0);
m_A.push_back(empty);
m_B.push_back(empty);
m_C.push_back(empty);
m_D.push_back(empty);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4LogicalVolume* PISTA::BuildTrapezoidDetector(){
if(!m_TrapezoidDetector){
// Definittion of the volume containing the sensitive detectors
G4Trap* solidTrapezoid = new G4Trap("PISTA",
TrapezoidLength*0.5, 0*deg, 0*deg,
TrapezoidHeight*0.5, TrapezoidBaseLarge*0.5, TrapezoidBaseSmall*0.5, 0*deg,
TrapezoidHeight*0.5, TrapezoidBaseLarge*0.5, TrapezoidBaseSmall*0.5, 0*deg);
G4LogicalVolume* logicTrapezoid = new G4LogicalVolume(solidTrapezoid,
MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"),
"PISTA",
0,0,0);
G4VisAttributes* TrapezoidVisAtt = new G4VisAttributes(G4Colour(0,0,0,0.5));
TrapezoidVisAtt->SetForceWireframe(true);
logicTrapezoid->SetVisAttributes(TrapezoidVisAtt);
// First stage silicon detector
G4ThreeVector positionFirstStage = G4ThreeVector(0,0,-0.5*TrapezoidLength + 0.5*FirstStageThickness);
G4Trap* solidFirstStage = new G4Trap("solidFirstSatge",
FirstStageThickness*0.5, 0*deg, 0*deg,
TrapezoidHeight*0.5, TrapezoidBaseLarge*0.5, TrapezoidBaseSmall*0.5, 0*deg,
TrapezoidHeight*0.5, TrapezoidBaseLarge*0.5, TrapezoidBaseSmall*0.5, 0*deg);
G4LogicalVolume* logicFirstStage = new G4LogicalVolume(solidFirstStage,
MaterialManager::getInstance()->GetMaterialFromLibrary("Si"),
"logicFirstStage",
0,0,0);
new G4PVPlacement(0,
positionFirstStage,
logicFirstStage,
"PISTA_FirstStage",
logicTrapezoid,
false,
0);
// Set First Stage sensitive
logicFirstStage->SetSensitiveDetector(m_FirstStageScorer);
// Visualisation of First Stage strips
//G4VisAttributes* FirstStageVisAtt = new G4VisAttributes(G4Colour(0.2,0.8,0.5));
G4VisAttributes* FirstStageVisAtt = new G4VisAttributes(G4Colour(0.5,0.5,0.55,0.8));
logicFirstStage->SetVisAttributes(FirstStageVisAtt);
//////
// Second stage silicon detector
G4ThreeVector positionSecondStage = G4ThreeVector(0,0,-0.5*TrapezoidLength+DistanceBetweenSi+0.5*SecondStageThickness);
G4Trap* solidSecondStage = new G4Trap("solidSecondSatge",
SecondStageThickness*0.5, 0*deg, 0*deg,
TrapezoidHeight*0.5, TrapezoidBaseLarge*0.5, TrapezoidBaseSmall*0.5, 0*deg,
TrapezoidHeight*0.5, TrapezoidBaseLarge*0.5, TrapezoidBaseSmall*0.5, 0*deg);
G4LogicalVolume* logicSecondStage = new G4LogicalVolume(solidSecondStage,
MaterialManager::getInstance()->GetMaterialFromLibrary("Si"),
"logicSecondStage",
0,0,0);
new G4PVPlacement(0,
positionSecondStage,
logicSecondStage,
"PISTA_SecondStage",
logicTrapezoid,
false,
0);
// Set Second Stage sensitive
logicSecondStage->SetSensitiveDetector(m_SecondStageScorer);
// Visualisation of Second Stage strips
G4VisAttributes* SecondStageVisAtt = new G4VisAttributes(G4Colour(0.5,0.5,0.5));
logicSecondStage->SetVisAttributes(SecondStageVisAtt);
m_TrapezoidDetector = logicTrapezoid;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Virtual Method of NPS::VDetector class
// Read stream at Configfile to pick-up parameters of detector (Position,...)
// Called in DetecorConstruction::ReadDetextorConfiguration Method
void PISTA::ReadConfiguration(NPL::InputParser parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("PISTA");
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detectors found " << endl;
vector<string> cart = {"POS_A", "POS_B", "POS_C", "POS_D"};
for(unsigned int i = 0 ; i < blocks.size() ; i++){
if(blocks[i]->HasTokenList(cart)){
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << endl << "//// PISTA " << i+1 << endl;
G4ThreeVector A = NPS::ConvertVector(blocks[i]->GetTVector3("POS_A","mm"));
G4ThreeVector B = NPS::ConvertVector(blocks[i]->GetTVector3("POS_B","mm"));
G4ThreeVector C = NPS::ConvertVector(blocks[i]->GetTVector3("POS_C","mm"));
G4ThreeVector D = NPS::ConvertVector(blocks[i]->GetTVector3("POS_D","mm"));
AddDetector(A,B,C,D);
}
else if(blocks[i]->HasTokenList(sphe)){
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << endl << "//// PISTA " << i+1 << endl;
double R = blocks[i]->GetDouble("R","mm");
double Theta = blocks[i]->GetDouble("Theta","deg");
double Phi = blocks[i]->GetDouble("Phi","deg");
}
else{
cout << "NPS ERROR: check your input file formatting " << endl;
exit(1);
}
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Construct detector and inialise sensitive part.
// Called After DetecorConstruction::AddDetector Method
void PISTA::ConstructDetector(G4LogicalVolume* world){
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
int NumberOfTelescopes = m_DefinitionType.size();
G4RotationMatrix* Rot = NULL;
G4ThreeVector Det_pos = G4ThreeVector(0,0,0);
G4ThreeVector u = G4ThreeVector(0,0,0);
G4ThreeVector v = G4ThreeVector(0,0,0);
G4ThreeVector w = G4ThreeVector(0,0,0);
G4ThreeVector center = G4ThreeVector(0,0,0);
for (int i = 0 ; i < NumberOfTelescopes ; i++) {
if(m_DefinitionType[i]){
// (u,v,w) unitary vector associated to telescope referencial
// (u,v) // to silicon plan
// w perpendicular to (u,v) plan
u = m_B[i] - m_A[i];
u = u.unit();
//v = m_D[i] - m_A[i];
v = (m_C[i] + m_D[i])/2 - (m_A[i] + m_B[i])/2;
v = v.unit();
w = u.cross(v);
w = w.unit();
center = (m_A[i] + m_B[i] + m_C[i] + m_D[i])/4;
Rot = new G4RotationMatrix(u,v,w);
Det_pos = w * TrapezoidLength * 0.5 + center;
}
else{
G4double wX = m_R[i] * sin(m_Theta[i] ) * cos(m_Phi[i] ) ;
G4double wY = m_R[i] * sin(m_Theta[i] ) * sin(m_Phi[i] ) ;
//G4double wZ = TrapezoidHeight*0.5 + m_R[i] * cos(m_Theta[i] ) ;
G4double wZ = m_R[i] * cos(m_Theta[i] ) ;
Det_pos = G4ThreeVector(wX, wY, wZ) ;
// So the face of the detector is at R instead of the middle
Det_pos+=Det_pos.unit()*PISTA_NS::TrapezoidLength*0.5;
// Building Detector reference frame
G4double ii = cos(m_Theta[i]) * cos(m_Phi[i]);
G4double jj = cos(m_Theta[i]) * sin(m_Phi[i]);
G4double kk = -sin(m_Theta[i]);
G4ThreeVector Y(ii,jj,kk);
w = Det_pos.unit();
u = w.cross(Y);
v = w.cross(u);
v = v.unit();
u = u.unit();
Rot = new G4RotationMatrix(u,v,w);
new G4PVPlacement(G4Transform3D(*Rot,Det_pos),
BuildTrapezoidDetector(),
"PISTA",world,false,i+1);
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Add Detector branch to the EventTree.
// Called After DetecorConstruction::AddDetector Method
void PISTA::InitializeRootOutput(){
RootOutput *pAnalysis = RootOutput::getInstance();
TTree *pTree = pAnalysis->GetTree();
if(!pTree->FindBranch("PISTA")){
pTree->Branch("PISTA", "TPISTAData", &m_Event) ;
}
pTree->SetBranchAddress("PISTA", &m_Event) ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Read sensitive part and fill the Root tree.
// Called at in the EventAction::EndOfEventAvtion
void PISTA::ReadSensitive(const G4Event* ){
m_Event->Clear();
///////////
DSSDScorers::PS_Rectangle* FirstStageScorer= (DSSDScorers::PS_Rectangle*) m_FirstStageScorer->GetPrimitive(0);
unsigned int sizeDEFront = FirstStageScorer->GetWidthMult();
for(unsigned int i = 0 ; i < sizeDEFront ; i++){
double EnergyFront = RandGauss::shoot(FirstStageScorer->GetEnergyWidth(i), DE_ResoEnergy);
if(EnergyFront>EnergyThreshold){
double Time = RandGauss::shoot(FirstStageScorer->GetTimeWidth(i), ResoTime);
int DetNbr = FirstStageScorer->GetDetectorWidth(i);
int StripFront = 92-FirstStageScorer->GetStripWidth(i);
m_Event->SetPISTA_DE_DetectorNbr(DetNbr);
m_Event->SetPISTA_DE_StripNbr(StripFront);
m_Event->SetPISTA_DE_StripEnergy(EnergyFront);
m_Event->SetPISTA_DE_StripTime(Time);
}
}
unsigned int sizeDEBack = FirstStageScorer->GetLengthMult();
for(unsigned int i = 0 ; i < sizeDEBack ; i++){
double EnergyBack = RandGauss::shoot(FirstStageScorer->GetEnergyLength(i), DE_ResoEnergy);
if(EnergyBack>EnergyThreshold){
double Time = RandGauss::shoot(FirstStageScorer->GetTimeLength(i), ResoTime);
int DetNbr = FirstStageScorer->GetDetectorLength(i);
m_Event->SetPISTA_DE_BackEnergy(EnergyBack);
m_Event->SetPISTA_DE_BackTime(Time);
/*for(unsigned int i = 0 ; i < sizeDEFront ; i++){
double EnergyFront = RandGauss::shoot(FirstStageScorer->GetEnergyWidth(i), DE_ResoEnergy);
double EnergyBack = RandGauss::shoot(FirstStageScorer->GetEnergyLength(i), DE_ResoEnergy);
if(EnergyFront>EnergyThreshold){
double Time = RandGauss::shoot(FirstStageScorer->GetTimeLength(i), ResoTime);
int DetNbr = FirstStageScorer->GetDetectorWidth(i);
int StripFront = 92-FirstStageScorer->GetStripWidth(i);
m_Event->SetPISTA_DE(DetNbr, StripFront, EnergyFront, EnergyBack, Time, Time);
m_Event->SetPISTA_DE_BackDetector(DetNbr);
}
}*/
FirstStageScorer->clear();
///////////
DSSDScorers::PS_Rectangle* SecondStageScorer= (DSSDScorers::PS_Rectangle*) m_SecondStageScorer->GetPrimitive(0);
unsigned int sizeEFront = SecondStageScorer->GetLengthMult();
for(unsigned int i = 0 ; i < sizeEFront ; i++){
double EnergyFront = RandGauss::shoot(SecondStageScorer->GetEnergyLength(i), E_ResoEnergy);
if(EnergyFront>EnergyThreshold){
double Time = RandGauss::shoot(SecondStageScorer->GetTimeLength(i), ResoTime);
int DetNbr = SecondStageScorer->GetDetectorLength(i);
int StripFront = SecondStageScorer->GetStripLength(i);
m_Event->SetPISTA_E_DetectorNbr(DetNbr);
m_Event->SetPISTA_E_StripNbr(StripFront);
m_Event->SetPISTA_E_StripEnergy(EnergyFront);
m_Event->SetPISTA_E_StripTime(Time);
}
}
unsigned int sizeEBack = SecondStageScorer->GetWidthMult();
for(unsigned int i = 0 ; i < sizeEBack ; i++){
double EnergyBack = RandGauss::shoot(SecondStageScorer->GetEnergyWidth(i), E_ResoEnergy);
if(EnergyBack>EnergyThreshold){
double Time = RandGauss::shoot(SecondStageScorer->GetTimeWidth(i), ResoTime);
int DetNbr = SecondStageScorer->GetDetectorWidth(i);
m_Event->SetPISTA_E_BackEnergy(EnergyBack);
m_Event->SetPISTA_E_BackTime(Time);
/*for(unsigned int i = 0 ; i < sizeEFront ; i++){
double EnergyFront = RandGauss::shoot(SecondStageScorer->GetEnergyLength(i), E_ResoEnergy);
double EnergyBack = RandGauss::shoot(SecondStageScorer->GetEnergyWidth(i), E_ResoEnergy);
if(EnergyFront>EnergyThreshold){
double Time = RandGauss::shoot(SecondStageScorer->GetTimeLength(i), ResoTime);
int DetNbr = SecondStageScorer->GetDetectorLength(i);
int StripFront = SecondStageScorer->GetStripLength(i);
m_Event->SetPISTA_E(DetNbr, StripFront, EnergyFront, EnergyBack, Time, Time);
m_Event->SetPISTA_E_BackDetector(DetNbr);
}
}*/
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
////////////////////////////////////////////////////////////////
void PISTA::InitializeScorers() {
// This check is necessary in case the geometry is reloaded
bool already_exist = false;
m_FirstStageScorer = CheckScorer("FirstStageScorer",already_exist) ;
m_SecondStageScorer = CheckScorer("SecondStageScorer",already_exist) ;
if(already_exist)
return ;
// Otherwise the scorer is initialised
G4VPrimitiveScorer* FirstStageScorer = new DSSDScorers::PS_Rectangle("FirstStageScorer",1,
TrapezoidBaseLarge,
TrapezoidHeight,
G4VPrimitiveScorer* SecondStageScorer = new DSSDScorers::PS_Rectangle("SecondStageScorer",1,
TrapezoidBaseLarge,
TrapezoidHeight,
G4VPrimitiveScorer* InteractionFirstStage = new InteractionScorers::PS_Interactions("InteractionFirstStage",ms_InterCoord,0);
//G4VPrimitiveScorer* InteractionSecondStage = new InteractionScorers::PS_Interactions("InteractionSecondStage",ms_InterCoord,0);
// Register it to the multifunctionnal detector
m_FirstStageScorer->RegisterPrimitive(FirstStageScorer);
m_FirstStageScorer->RegisterPrimitive(InteractionFirstStage);
m_SecondStageScorer->RegisterPrimitive(SecondStageScorer);
//m_SecondStageScorer->RegisterPrimitive(InteractionSecondStage);
G4SDManager::GetSDMpointer()->AddNewDetector(m_FirstStageScorer);
G4SDManager::GetSDMpointer()->AddNewDetector(m_SecondStageScorer);
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
NPS::VDetector* PISTA::Construct(){
return (NPS::VDetector*) new PISTA();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern"C" {
class proxy_nps_PISTA{
public:
proxy_nps_PISTA(){
NPS::DetectorFactory::getInstance()->AddToken("PISTA","PISTA");
NPS::DetectorFactory::getInstance()->AddDetector("PISTA",PISTA::Construct);
}
};
proxy_nps_PISTA p_nps_PISTA;
}