PISTA.cc

/*****************************************************************************
 * 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
  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;
}
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;
}