ScintillatorConstruction.cc

Go to the documentation of this file.

#include "ScintillatorConstruction.h"

#include <fwk/CentralConfig.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/MathConstants.h>

//G4 Classes
#include <G4RunManager.hh>
#include <G4RegionStore.hh>
#include <G4Box.hh>
#include <G4Tubs.hh>
#include <G4LogicalVolume.hh>
#include <G4ThreeVector.hh>
#include <G4PVPlacement.hh>
#include <G4Transform3D.hh>
#include <G4RotationMatrix.hh>
#include <G4VisAttributes.hh>
#include <G4Colour.hh>
#include <G4Material.hh>
#include <G4Element.hh>
#include <G4MaterialPropertiesTable.hh>

#include <utl/TabularStream.h>
#include <utl/ConfigUtils.h>


using namespace fwk;
using namespace utl;
using namespace EdepSimulatorAG;
 
ScintillatorConstruction::ScintillatorConstruction ( )  
{
  fLogicalWorld  = 0;
  fPhysicalWorld = 0;
  fLogicalScintillator   = 0;
  fPhysicalScintillator  = 0;
  L_LogScintillator = 0;
  L_PhyScintillator = 0;
  R_LogScintillator = 0;
  R_PhyScintillator = 0;

  INFO("\n++++\nGetting parameters from XML file...");
  Branch config(fwk::CentralConfig::GetInstance()->GetTopBranch("MdEdepSimulator"));

  fLog.Configure(config);
  //Set energy units to MeV (otherwise default is eV)
  fUnits.SetEnergyDefault( utl::MeV, "MeV");
  fUnits.SetLengthDefault( utl::cm , "cm");
  fUnits.Configure(config);
   
  LoadConfig( config.GetChild("scintillator"), "density", fScintillatorDensity, 1.06*utl::g/utl::cm3 );

  std::vector<G4double> vScintiSize, vScintiSizeDefault;

  vScintiSizeDefault.push_back(   4.01*utl::cm );
  vScintiSizeDefault.push_back( 400.00*utl::cm );
  vScintiSizeDefault.push_back(   1.00*utl::cm );

  LoadConfig( config.GetChild("scintillator"), "size", vScintiSize,vScintiSizeDefault );

  fScint_Side_X = vScintiSize.at(0)/2;
  fScint_Side_Y = vScintiSize.at(1)/2;
  fScint_Side_Z = vScintiSize.at(2)/2;

  fWorld_Side_X = 3*fScint_Side_X;//To host 3 strips (central + 2 neighbours)
  fWorld_Side_Y = fScint_Side_Y;
  fWorld_Side_Z = fScint_Side_Z;

  //Convert from AugerUnits to CLHEP units
  fScintillatorDensity  *= ((CLHEP::g/CLHEP::cm3)/(utl::g/utl::cm3));

  fWorld_Side_X *= (CLHEP::m / utl::m); 
  fWorld_Side_Y *= (CLHEP::m / utl::m);
  fWorld_Side_Z *= (CLHEP::m / utl::m);

  fScint_Side_X *= (CLHEP::m / utl::m);
  fScint_Side_Y *= (CLHEP::m / utl::m);
  fScint_Side_Z *= (CLHEP::m / utl::m);

  G4cerr << "X World size     " << 2*fWorld_Side_X/CLHEP::m << " (m)" << G4endl;
  G4cerr << "Y World size     " << 2*fWorld_Side_Y/CLHEP::m << " (m)" << G4endl;
  G4cerr << "Z World size     " << 2*fWorld_Side_Z/CLHEP::m << " (m)" << G4endl;

  G4cerr << "Scintillator Density " << fScintillatorDensity/(CLHEP::g/CLHEP::cm3) << " g/cm3"<< G4endl;

  G4cerr << "X Scinti size     " << 2*fScint_Side_X/CLHEP::m << " (m)" << G4endl;
  G4cerr << "Y Scinti size     " << 2*fScint_Side_Y/CLHEP::m << " (m)" << G4endl;
  G4cerr << "Z Scinti size     " << 2*fScint_Side_Z/CLHEP::m << " (m)" << G4endl;

  //Creates materials after unit conversion
  CreateMaterials();
}

ScintillatorConstruction::~ScintillatorConstruction ( ){}

G4VPhysicalVolume *
ScintillatorConstruction::Construct ( )
{
  if(!fPhysicalWorld){
    G4Colour  red   ( 1, 0, 0 );
    G4Colour  green ( 0, 1, 0 );
    G4Colour  blue  ( 0, 0, 1 );

    G4ThreeVector placement;
    
    INFO("\n++++\nCreating World for scintillator simulation");
    //World Box
    G4Box *WorldBox = new G4Box("WorldBox", fWorld_Side_X, fWorld_Side_Y, fWorld_Side_Z);
    fLogicalWorld  = new G4LogicalVolume ( WorldBox, fWorldMaterial, "LogWorld", 0, 0, 0 );
    fPhysicalWorld = new G4PVPlacement ( 0, G4ThreeVector (  ), fLogicalWorld, "WorldPhys", 0, false, 0 );
    G4VisAttributes *fLogicalWorldVisAtt = new G4VisAttributes ( red );
    fLogicalWorld->SetVisAttributes ( fLogicalWorldVisAtt );
    //fLogicalWorld->SetVisAttributes ( G4VisAttributes::Invisible );
    
    INFO("\n++++\nCreating scintillator + first neighbours");
    //Scintillator Box
    G4Box *ScintillatorBox = new G4Box ( "ScintillatorBox", fScint_Side_X, fScint_Side_Y, fScint_Side_Z );
    fLogicalScintillator = new G4LogicalVolume ( ScintillatorBox, fScintillatorMaterial, "ScintillatorLog" );
    
    G4double scintillator_pos_x;
    G4double scintillator_pos_y;
    G4double scintillator_pos_z;

    //Central scintillator (where particle is injected)
    scintillator_pos_x = 0.0 * utl::m;
    scintillator_pos_y = 0.0 * utl::m;
    scintillator_pos_z = 0.0 * utl::m;//-fScint_Side_Z;
    placement = G4ThreeVector ( scintillator_pos_x, scintillator_pos_y, scintillator_pos_z );
    fPhysicalScintillator   = new G4PVPlacement ( 0, placement, fLogicalScintillator, "CentralScint", fLogicalWorld, false, 0 );//CENTRAL
    fLogicalScintillator->SetVisAttributes ( new G4VisAttributes ( green ) );
    
    //First neighbouring scintillators
    L_LogScintillator = new G4LogicalVolume ( ScintillatorBox, fScintillatorMaterial, "L_ScintillatorLog" );//LEFT
    placement = G4ThreeVector ( scintillator_pos_x - 2*fScint_Side_X, scintillator_pos_y, scintillator_pos_z );
    L_PhyScintillator = new G4PVPlacement ( 0, placement, L_LogScintillator, "LeftScint", fLogicalWorld, false, 0 );
    L_LogScintillator  ->SetVisAttributes ( new G4VisAttributes ( blue ) );
    
    R_LogScintillator = new G4LogicalVolume ( ScintillatorBox, fScintillatorMaterial, "R_ScintillatorLog" );//RIGHT
    placement = G4ThreeVector ( scintillator_pos_x + 2*fScint_Side_X, scintillator_pos_y, scintillator_pos_z );
    R_PhyScintillator = new G4PVPlacement ( 0, placement, R_LogScintillator, "RightScint", fLogicalWorld, false, 0 );
    R_LogScintillator  ->SetVisAttributes ( new G4VisAttributes ( blue ) );
  }
  return fPhysicalWorld;
}

void 
ScintillatorConstruction::CreateMaterials (  )
{
  INFO("\n++++\nCreating materials...");
  G4double  a(0), z(0);
  G4int     nel(0); 
 
  H = new G4Element ( "Hydrogen" , "H",  z = 1., a = 1.010 * CLHEP::g / CLHEP::mole );
  C = new G4Element ( "Carbon"   , "C",  z = 6., a = 12.01 * CLHEP::g / CLHEP::mole );
  N = new G4Element ( "Nitrogen" , "N",  z = 7., a = 14.01 * CLHEP::g / CLHEP::mole );
  O = new G4Element ( "Oxygen"   , "O",  z = 8., a = 15.99 * CLHEP::g / CLHEP::mole );
  
  G4Material *Polystyrene = new G4Material ( "Polystyrene",  fScintillatorDensity, nel = 2 );
  Polystyrene->AddElement ( C, 8 );
  Polystyrene->AddElement ( H, 8 );
    
  //Malargue Scintillator 
  fScintillatorMaterial = new G4Material ( "Scintillator", fScintillatorDensity , nel = 1 );
  fScintillatorMaterial->AddMaterial ( Polystyrene,  100.0 * CLHEP::perCent );

  //World material (Air)
  fWorldMaterial = new G4Material ( "Air_MdEdep", 1.29 * CLHEP::mg / CLHEP::cm3, nel = 2 );
  fWorldMaterial->AddElement ( N, 70. * CLHEP::perCent );
  fWorldMaterial->AddElement ( O, 30. * CLHEP::perCent );

  return;
}