/*****************************************************************************
* 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 "G4Trap.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 "DSSDScorers.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 ResoTime = 0.2*ns ;
const double DE_ResoEnergy = 0.040*MeV ;
const double E_ResoEnergy = 0.018*MeV ;
// Trapezoid dimension PISTA1
//const double TrapezoidBaseLarge = 72.3*mm;
//const double TrapezoidBaseSmall = 41.0*mm;
//const double TrapezoidHeight = 57.7*mm;
//const double TrapezoidLength = 1*cm;
//const double FirstStageThickness = 300*um;
//const double SecondStageThickness = 1.5*mm;
//const double DistanceBetweenSi = 4*mm;
//const double FirstStageNbrOfStrips = 91;
//const double SecondStageNbrOfStrips = 57;
// Trapezoid dimension PISTA2
const double TrapezoidBaseLarge = 57*mm;
const double TrapezoidBaseSmall = 38.0*mm;
const double TrapezoidHeight = 31*mm;
const double TrapezoidLength = 1*cm;
const double FirstStageThickness = 300*um;
const double SecondStageThickness = 1.2*mm;
const double DistanceBetweenSi = 4*mm;
const double FirstStageNbrOfStrips = 49;
const double SecondStageNbrOfStrips = 45;
}
using namespace PISTA_NS;
//....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;
}
return m_TrapezoidDetector;
}
//....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"};
vector<string> sphe = {"R","Theta","Phi"};
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");
AddDetector(R,Theta,Phi);
}
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){
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();
///////////
// DE
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_BackDetector(DetNbr);
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();
///////////
// E
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_BackDetector(DetNbr);
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);
}
}*/
SecondStageScorer->clear();
}
//....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,
1,FirstStageNbrOfStrips);
G4VPrimitiveScorer* SecondStageScorer = new DSSDScorers::PS_Rectangle("SecondStageScorer",1,
TrapezoidBaseLarge,
TrapezoidHeight,
SecondStageNbrOfStrips,1);
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);
}
//....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;
}