Skip to content
Snippets Groups Projects
MUST2Array.cc 80.02 KiB
/*****************************************************************************
 * Copyright (C) 2009-2013   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: Adrien MATTA  contact address: matta@ipno.in2p3.fr       *
 *                                                                           *
 * Creation Date  : January 2009                                             *
 * Last update    : May 2014                                                 *
 *---------------------------------------------------------------------------*
 * Decription:                                                               *
 *  This file describe the MUST2 charge particle Detector                    *
 *                                                                           *
 *---------------------------------------------------------------------------*
 * Comment:                                                                  *
 * MUST2 is a modular array made of Telescope (1 to 8 telescope). Each       *
 *  Telescope is made of Three Stage:                                        *
 *  - A 300um Silicium, double-sided strip                                   *
 *  - 16 Si(Li) pad                                                          *
 *  - 16 CsI scintillator Crystal                                            *
 *****************************************************************************/
#include"MUST2Array.hh"

//G4 Geometry object
#include "G4Trd.hh"
#include "G4Box.hh"
#include "G4Trap.hh"
#include "G4ExtrudedSolid.hh"     
#include "G4SubtractionSolid.hh" 

//G4 various object
#include "G4MaterialTable.hh"
#include "G4Element.hh"
#include "G4ElementTable.hh"
#include "G4RotationMatrix.hh"
#include "G4Transform3D.hh"
#include "G4PVPlacement.hh"
#include "G4VisAttributes.hh"
#include "G4Colour.hh"
#include "G4PVDivision.hh"

//G4 sensitive
#include "ObsoleteGeneralScorers.hh"
#include "MUST2Scorers.hh"

//Not G4
#include "RootOutput.h"
#include "CLHEP/Random/RandGauss.h"
#include "sstream"
#include "string"
#include <cmath>

using namespace std  ;
using namespace CLHEP   ;
using namespace MUST2   ;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// MUST2Array Specific Method
MUST2Array::MUST2Array(){
  m_Event = new TMust2Data() ;
  InitializeMaterial();

  // Common vis attribute:
  m_PCBVisAtt = new G4VisAttributes(G4Colour(0.2, 0.5, 0.2)) ; 
  m_SiliconVisAtt = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3)) ;
  m_SiLiVisAtt = new G4VisAttributes(G4Colour(0.3, 1, 0.3));
  m_CsIVisAtt = new G4VisAttributes(G4Colour(1, 0.5, 0));
  m_CoolingVisAtt = new G4VisAttributes(G4Colour(0.5, 0.25, 0.25));
}

MUST2Array::~MUST2Array(){
  delete m_MaterialAluminium;
  delete m_MaterialIron;
  delete m_MaterialCsI;
  delete m_MaterialVacuum ;
  delete m_MaterialMyl;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MUST2Array::AddTelescope(   G4ThreeVector X1_Y1,
    G4ThreeVector X128_Y1,
    G4ThreeVector X1_Y128,
    G4ThreeVector X128_Y128,
    bool wSi,
    bool wSiLi,
    bool wCsI)
{
  m_DefinitionType.push_back(true);

  m_X1_Y1     .push_back(X1_Y1);
  m_X128_Y1   .push_back(X128_Y1);
  m_X1_Y128   .push_back(X1_Y128);
  m_X128_Y128 .push_back(X128_Y128);
  m_wSi.push_back(wSi);
  m_wSiLi.push_back(wSiLi);
  m_wCsI.push_back(wCsI);

  m_R.push_back(0);
  m_Theta.push_back(0);
  m_Phi.push_back(0);
  m_beta_u.push_back(0);
  m_beta_v.push_back(0);
  m_beta_w.push_back(0);
}

void MUST2Array::AddTelescope(   G4double R,
    G4double Theta,
    G4double Phi,
    G4double beta_u,
    G4double beta_v,
    G4double beta_w,
    bool wSi,
    bool wSiLi,
    bool wCsI)
{
  G4ThreeVector empty = G4ThreeVector(0, 0, 0);

  m_DefinitionType.push_back(false);

  m_R.push_back(R);
  m_Theta.push_back(Theta);
  m_Phi.push_back(Phi);
  m_beta_u.push_back(beta_u);
  m_beta_v.push_back(beta_v);
  m_beta_w.push_back(beta_w);
  m_wSi.push_back(wSi);
  m_wSiLi.push_back(wSiLi);
  m_wCsI.push_back(wCsI);

  m_X1_Y1     .push_back(empty) ;
  m_X128_Y1   .push_back(empty) ;
  m_X1_Y128   .push_back(empty) ;
  m_X128_Y128 .push_back(empty) ;

}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4LogicalVolume* MUST2Array::GetLogicalVolumeMUST2Short(){
   ////////////////////////////////////////////////////////////////
  ////////////// Starting Volume Definition //////////////////////
  ////////////////////////////////////////////////////////////////
  G4Trd*           solidMM = new G4Trd("MUST2Telescope" + DetectorNumber, 0.5*FaceFront, 0.5*FaceBack, 0.5*FaceFront, 0.5*FaceBack, 0.5*Length);
  G4LogicalVolume* logicMM = new G4LogicalVolume(solidMM, m_MaterialIron, "MUST2Telescope" + DetectorNumber, 0, 0, 0);
  G4String Name = "MUST2Telescope" + DetectorNumber ;

  new G4PVPlacement(	G4Transform3D(*MMrot, MMpos),
      logicMM               ,
      Name                  ,
      world                 ,
      false                 ,
      0                     );

  if (m_non_sensitive_part_visiualisation){
    G4VisAttributes* FrameVisAtt = new G4VisAttributes(G4Colour(0.80, 0.80, 0.80));
    FrameVisAtt->SetForceWireframe(true);
    logicMM->SetVisAttributes(FrameVisAtt)  ;
  }

  else logicMM->SetVisAttributes(G4VisAttributes::Invisible)               ;

  G4ThreeVector positionVacBox = G4ThreeVector(0, 0, VacBox_PosZ);

  G4Trd* solidVacBox = 
    new G4Trd("solidVacBox", 0.5*ActiveWaferSize, 0.5*CsIFaceFront, 0.5*ActiveWaferSize, 0.5*CsIFaceFront, 0.5*VacBoxThickness);

  G4LogicalVolume* logicVacBox = new G4LogicalVolume(solidVacBox, m_MaterialVacuum, "logicVacBox", 0, 0, 0);


  new G4PVPlacement(0, positionVacBox, logicVacBox, Name + "_VacBox", logicMM, false, 0);

  logicVacBox->SetVisAttributes(G4VisAttributes::Invisible);


  ////////////////////////////////////////////////////////////////
  /////////////////Si Strip Construction//////////////////////////
  ////////////////////////////////////////////////////////////////
  // A no rotation matrix is always handy ;) 
  G4RotationMatrix* norotation = new  G4RotationMatrix();   

  // Create the front support frame for Si and SiLi
  vector<G4TwoVector> polygon;
  for(unsigned int i = 0 ; i < 4 ; i++){   
    G4TwoVector Point(PCBPointsX[i],PCBPointsY[i]);   
    polygon.push_back(Point);  
  } 

  // The PCB lenth is used to rescale different part of the detector
  G4double PCBlength = PCBPointsX[0] - PCBPointsX[1];

  G4ExtrudedSolid* solidFrontFrameBase = new G4ExtrudedSolid("Must2 Front Frame",
      polygon,   
      FrontFrameThickness*0.5,
      G4TwoVector(0,0),FrontFrameFrontSize/PCBlength,
      G4TwoVector(0,0),FrontFrameBackSize/PCBlength); 

  G4ExtrudedSolid* solidCavityShape1 = new G4ExtrudedSolid("Cavity1",
      polygon,   
      FrontFrameThickness*0.5+PCBThickness,
      G4TwoVector(0,0),CavityShape1FrontSize/PCBlength,
      G4TwoVector(0,0),CavityShape1BackSize/PCBlength); 

  G4SubtractionSolid* solidFrontFrame1 = new G4SubtractionSolid("Frame1",
      solidFrontFrameBase,
      solidCavityShape1,
      G4Transform3D(*norotation,G4ThreeVector(0,0,0))); 
  G4Box* solidCavityShape2 = new G4Box("Cavity2", 0.5*CavityShape2Width, (0.5+0.5)*CavityShape2Width, 0.5*CavityShape2Thickness);

  G4SubtractionSolid* solidFrontFrame2 = new G4SubtractionSolid("Frame2",
      solidFrontFrame1,
      solidCavityShape2,
      G4Transform3D(*norotation,G4ThreeVector(0,0.5*CavityShape2Width,CavityShape2Center))); 

  G4ExtrudedSolid* solidCutOut = new G4ExtrudedSolid("Cutout",
      polygon,   
      FrontFrameThickness*0.8,
      G4TwoVector(0,0),1.1*FrontFrameBackSize/PCBlength,
      G4TwoVector(0,0),1.1*FrontFrameBackSize/PCBlength); 


  G4SubtractionSolid* solidFrontFrame3 = new G4SubtractionSolid("Must2 Front Frame",
      solidFrontFrame2,
      solidCutOut,
      G4Transform3D(*norotation,G4ThreeVector(0,FrontFrameBackSize*0.5*1.1+66*mm,0))); 

  G4SubtractionSolid* solidFrontFrame = new G4SubtractionSolid("Must2 Front Frame",
      solidFrontFrame3,
      solidCutOut,
      G4Transform3D(*norotation,G4ThreeVector(-FrontFrameBackSize*0.5*1.1-66*mm,0,0))); 


  G4LogicalVolume* logicFrontFrame = new G4LogicalVolume(solidFrontFrame,m_MaterialAluminium,"Must2 Front Frame", 0, 0, 0); 

  new G4PVPlacement(G4Transform3D(*norotation,G4ThreeVector(0,0,Silicon_PosZ+FrontFrameThickness*0.5+0.5*PCBThickness)),
      logicFrontFrame, 
      "Must2 Front Fame",
      logicMM,
      false,
      0); 


  // Electronic bloc
  G4Box* solidMUFEE  = new G4Box("solidMUFEE", 0.5*MUFEESize, 0.5*MUFEESize, 0.5*MUFEEThickness);
  G4Box* solidCooling = new G4Box("solidCooling", 0.5*MUFEESize, 0.5*MUFEESize, 0.5*CoolingThickness);
  G4ExtrudedSolid* solidElecFrameShort = new G4ExtrudedSolid("SolidElecFrame",
      polygon,   
      CoolingThickness*0.5,
      G4TwoVector(0,0),FrontFrameBackSize/PCBlength,
      G4TwoVector(0,0),FrontFrameBackSize/PCBlength); 


  G4LogicalVolume* logicMUFEE   = new G4LogicalVolume(solidMUFEE,m_MaterialAluminium,"Must2 MUFEE", 0, 0, 0); 
  G4LogicalVolume* logicCooling = new G4LogicalVolume(solidCooling,m_MaterialAluminium,"Must2 Cooling", 0, 0, 0); 
  G4LogicalVolume* logicElecFrame = new G4LogicalVolume(solidElecFrameShort,m_MaterialAluminium,"Must2 Frame", 0, 0, 0); 

  logicMUFEE->SetVisAttributes(m_PCBVisAtt);
  logicCooling->SetVisAttributes(m_CoolingVisAtt);

  G4double CoolingPosition = FrontFrameThickness+Silicon_PosZ+0.5*CoolingThickness;

  new G4PVPlacement(G4Transform3D(*norotation,G4ThreeVector(0,0,CoolingPosition)),
      logicCooling, 
      "Must2 Cooling",
      logicMM,
      false,
      0); 

  new G4PVPlacement(G4Transform3D(*norotation,G4ThreeVector(0,0,CoolingPosition)),
      logicElecFrame, 
      "Must2 Fame",
      logicMM,
      false,
      0); 


  new G4PVPlacement(G4Transform3D(*norotation,G4ThreeVector(0,0,CoolingPosition-0.5*MUFEEThickness-0.5*CoolingThickness)),
      logicMUFEE, 
      "Must2 Front Fame",
      logicMM,
      false,
      0); 

  new G4PVPlacement(G4Transform3D(*norotation,G4ThreeVector(0,0,CoolingPosition+0.5*MUFEEThickness+0.5*CoolingThickness)),
      logicMUFEE, 
      "Must2 Front Fame",
      logicMM,
      false,
      0); 




  if (wSi) {

    // Volume containing the Silicon detector
    G4ExtrudedSolid* solidSilicon = new G4ExtrudedSolid(Name,
        polygon,   
        PCBThickness*0.5,
        G4TwoVector(0,0),1,
        G4TwoVector(0,0),1); 

    G4LogicalVolume* logicSilicon = new G4LogicalVolume(solidSilicon,m_MaterialVacuum, Name, 0, 0, 0); 

    new G4PVPlacement(G4Transform3D(*norotation,G4ThreeVector(0,0,Silicon_PosZ)),
        logicSilicon, 
        "Silicon",
        logicMM,
        false,
        0); 

    logicSilicon->SetVisAttributes(G4VisAttributes::Invisible);

    // PCB Base
    G4ExtrudedSolid* solidPCBBase = new G4ExtrudedSolid("PCBBase",
        polygon,   
        PCBThickness*0.5,
        G4TwoVector(0,0),1,
        G4TwoVector(0,0),1); 

    // Wafer Shape 
    G4Box* solidWaferShape = new G4Box("solidWaferShape", 0.5*WaferSize, 0.5*WaferSize, 0.6*PCBThickness);

    // PCB

    G4SubtractionSolid* solidPCB = new G4SubtractionSolid("PCB",
        solidPCBBase,
        solidWaferShape,
        G4Transform3D(*norotation,G4ThreeVector(0,0,0))); 

    G4LogicalVolume* logicPCB = new G4LogicalVolume(solidPCB, m_MaterialSilicon, "logicPCB", 0, 0, 0);

    new G4PVPlacement(0, G4ThreeVector(0,0,0), logicPCB, Name + "_Silicon", logicSilicon, false, 0);

    logicPCB->SetVisAttributes(m_PCBVisAtt); 


    // Wafer
    G4Box* solidWafer = new G4Box("solidWaferShape", 0.5*WaferSize, 0.5*WaferSize, 0.5*WaferThickness);
    G4LogicalVolume* logicWafer = new G4LogicalVolume(solidWafer, m_MaterialSilicon, "logicWafer", 0, 0, 0);
    new G4PVPlacement(0, G4ThreeVector(0,0,0), logicWafer, Name + "_Silicon", logicSilicon, false, 0);


    // Active Wafer
    G4Box* solidActiveWafer = new G4Box("solidActiveWaferShape", 0.5*ActiveWaferSize, 0.5*ActiveWaferSize, 0.5*WaferThickness);
    G4LogicalVolume* logicActiveWafer = new G4LogicalVolume(solidActiveWafer, m_MaterialSilicon, "logicActiveWafer", 0, 0, 0);
    new G4PVPlacement(0, G4ThreeVector(0,0,0), logicActiveWafer, Name + "_Silicon", logicSilicon, false, 0);

    G4ThreeVector positionAluStripFront = G4ThreeVector(0, 0, 0.5*WaferThickness+0.5*AluStripThickness);
    G4ThreeVector positionAluStripBack  = G4ThreeVector(0, 0, -0.5*WaferThickness-0.5*AluStripThickness);

    G4Box*           solidAluStrip = new G4Box("AluBox", 0.5*ActiveWaferSize, 0.5*ActiveWaferSize, 0.5*AluStripThickness);
    G4LogicalVolume* logicAluStrip = new G4LogicalVolume(solidAluStrip, m_MaterialAluminium, "logicAluStrip", 0, 0, 0);

    new G4PVPlacement(0, positionAluStripFront, logicAluStrip, "AluStripFront", logicSilicon, false, 0);
    new G4PVPlacement(0, positionAluStripBack, logicAluStrip, "AluStripBack", logicSilicon, false, 0);

    logicAluStrip->SetVisAttributes(G4VisAttributes::Invisible);

    G4ThreeVector  positionSilicon = G4ThreeVector(0, 0, Silicon_PosZ);

    new G4PVPlacement(0, G4ThreeVector(0, 0, Silicon_PosZ), logicSilicon, "Silicon", logicMM, false, 0);

    // Set Silicon strip sensible
    logicActiveWafer->SetSensitiveDetector(m_StripScorer);

    // Visualisation of Silicon Strip
    logicWafer->SetVisAttributes(m_SiliconVisAtt);
    logicActiveWafer->SetVisAttributes(m_SiliconVisAtt);
  }
  ////////////////////////////////////////////////////////////////
  //////////////////// SiLi  Construction ////////////////////////
  ////////////////////////////////////////////////////////////////

  if (wSiLi) {
    G4double SiLiSpace = 8 * mm;
    G4RotationMatrix* rotSiLi = new G4RotationMatrix(0,0,0); 
    G4Box* solidSiLi = new G4Box("SiLi", 0.5*ActiveWaferSize+0.5*SiLiSpace, 0.5*ActiveWaferSize, 0.5*SiLiThickness);
    G4LogicalVolume* logicSiLi = new G4LogicalVolume(solidSiLi, m_MaterialAluminium, Name + "_SiLi" , 0, 0, 0);

    logicSiLi->SetVisAttributes(G4VisAttributes::Invisible);

    new G4PVPlacement(  G4Transform3D(*rotSiLi, G4ThreeVector(0,0,0) ) 	,
        logicSiLi ,
        Name + "_SiLi",
        logicVacBox ,
        false ,
        0);

    // SiLi are placed inside of the VacBox...
    // Left/Right define when looking to detector from Si to CsI
    G4double SiLi_HighY_Upper 	= 19.86 * mm;
    G4double SiLi_HighY_Center = 25.39 * mm;
    G4double SiLi_WidthX_Left  = 22.85 * mm;
    G4double SiLi_WidthX_Right = 24.9  * mm;
    G4double SiLi_ShiftX    	= 0.775 * mm;

    //	SiLi are organized by two group of 8 Up(9 to 15) and Down(1 to 8).
    G4ThreeVector ShiftSiLiUp 	 = G4ThreeVector(-0.25 * ActiveWaferSize - 0.5 * SiLiSpace, 0, 0)	;
    G4ThreeVector ShiftSiLiDown  = G4ThreeVector(0.25  * ActiveWaferSize + 0.5 * SiLiSpace, 0, 0)	;

    // SiLi : left side of SiLi detector
    G4Box* solidSiLi_LT  = new G4Box("SiLi_LT"  , 0.5*SiLi_WidthX_Left  , 0.5*SiLi_HighY_Upper   , 0.5*SiLiThickness);
    G4Box* solidSiLi_RT  = new G4Box("SiLi_RT"  , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Upper   , 0.5*SiLiThickness);
    G4Box* solidSiLi_LC1 = new G4Box("SiLi_LC1" , 0.5*SiLi_WidthX_Left  , 0.5*SiLi_HighY_Center  , 0.5*SiLiThickness);
    G4Box* solidSiLi_RC1 = new G4Box("SiLi_RC1" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Center  , 0.5*SiLiThickness);
    G4Box* solidSiLi_LB  = new G4Box("SiLi_LB"  , 0.5*SiLi_WidthX_Left  , 0.5*SiLi_HighY_Upper   , 0.5*SiLiThickness);
    G4Box* solidSiLi_RB  = new G4Box("SiLi_RB"  , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Upper   , 0.5*SiLiThickness);
    G4Box* solidSiLi_LC2 = new G4Box("SiLi_LC2" , 0.5*SiLi_WidthX_Left  , 0.5*SiLi_HighY_Center  , 0.5*SiLiThickness);
    G4Box* solidSiLi_RC2 = new G4Box("SiLi_RC2" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Center  , 0.5*SiLiThickness);

    G4LogicalVolume* logicSiLi_LT    = new G4LogicalVolume(solidSiLi_LT   , m_MaterialSilicon , "SiLi_LT"  , 0 , 0 , 0);
    G4LogicalVolume* logicSiLi_RT    = new G4LogicalVolume(solidSiLi_RT   , m_MaterialSilicon , "SiLi_RT"  , 0 , 0 , 0);
    G4LogicalVolume* logicSiLi_LC1   = new G4LogicalVolume(solidSiLi_LC1  , m_MaterialSilicon , "SiLi_LC1" , 0 , 0 , 0);
    G4LogicalVolume* logicSiLi_RC1   = new G4LogicalVolume(solidSiLi_RC1  , m_MaterialSilicon , "SiLi_RC1" , 0 , 0 , 0);
    G4LogicalVolume* logicSiLi_LB    = new G4LogicalVolume(solidSiLi_LB   , m_MaterialSilicon , "SiLi_LB"  , 0 , 0 , 0);
    G4LogicalVolume* logicSiLi_RB    = new G4LogicalVolume(solidSiLi_RB   , m_MaterialSilicon , "SiLi_RB"  , 0 , 0 , 0);
    G4LogicalVolume* logicSiLi_LC2   = new G4LogicalVolume(solidSiLi_LC2  , m_MaterialSilicon , "SiLi_LC2" , 0 , 0 , 0);
    G4LogicalVolume* logicSiLi_RC2   = new G4LogicalVolume(solidSiLi_RC2  , m_MaterialSilicon , "SiLi_RC2" , 0 , 0 , 0);

    G4double interSiLi = 0.5 * mm;

    // Top
    G4ThreeVector positionSiLi_LT_up = G4ThreeVector(  -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
        0.5 * SiLi_HighY_Upper  + SiLi_HighY_Center + 1.5 * interSiLi,
        0);

    positionSiLi_LT_up += ShiftSiLiUp ;

    G4ThreeVector positionSiLi_RT_up = G4ThreeVector(  0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
        0.5 * SiLi_HighY_Upper  + SiLi_HighY_Center + 1.5 * interSiLi,
        0);

    positionSiLi_RT_up += ShiftSiLiUp ;

    G4ThreeVector positionSiLi_LC1_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
        0.5 * SiLi_HighY_Center + 0.5 * interSiLi,
        0);

    positionSiLi_LC1_up += ShiftSiLiUp ;

    G4ThreeVector positionSiLi_RC1_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
        0.5 * SiLi_HighY_Center + 0.5 * interSiLi,
        0);

    positionSiLi_RC1_up += ShiftSiLiUp ;

    G4ThreeVector positionSiLi_LB_up = G4ThreeVector(  -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
        -0.5 * SiLi_HighY_Upper  - SiLi_HighY_Center - 1.5 * interSiLi ,
        0);

    positionSiLi_LB_up += ShiftSiLiUp ;

    G4ThreeVector positionSiLi_RB_up = G4ThreeVector(  0.5 * SiLi_WidthX_Right - SiLi_ShiftX ,
        -0.5 * SiLi_HighY_Upper  - SiLi_HighY_Center - 1.5 * interSiLi ,
        0);

    positionSiLi_RB_up += ShiftSiLiUp ;

    G4ThreeVector positionSiLi_LC2_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
        -0.5 * SiLi_HighY_Center - 0.5 * interSiLi,
        0);

    positionSiLi_LC2_up += ShiftSiLiUp ;

    G4ThreeVector positionSiLi_RC2_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX ,
        -0.5 * SiLi_HighY_Center - 0.5 * interSiLi  ,
        0);

    positionSiLi_RC2_up += ShiftSiLiUp ;


    // Down
    G4ThreeVector positionSiLi_LT_down = G4ThreeVector(   -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
        0.5 * SiLi_HighY_Upper  + SiLi_HighY_Center + 1.5 * interSiLi,
        0);

    positionSiLi_LT_down += ShiftSiLiDown ;

    G4ThreeVector positionSiLi_RT_down = G4ThreeVector(   0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
        0.5 * SiLi_HighY_Upper  + SiLi_HighY_Center + 1.5 * interSiLi,
        0);

    positionSiLi_RT_down += ShiftSiLiDown ;

    G4ThreeVector positionSiLi_LC1_down = G4ThreeVector(  -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
        0.5 * SiLi_HighY_Center + 0.5 * interSiLi,
        0);

    positionSiLi_LC1_down += ShiftSiLiDown ;

    G4ThreeVector positionSiLi_RC1_down = G4ThreeVector(  0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
        0.5 * SiLi_HighY_Center + 0.5 * interSiLi ,
        0);

    positionSiLi_RC1_down += ShiftSiLiDown ;

    G4ThreeVector positionSiLi_LB_down = G4ThreeVector(   -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
        -0.5 * SiLi_HighY_Upper  - SiLi_HighY_Center - 1.5 * interSiLi,
        0);

    positionSiLi_LB_down += ShiftSiLiDown ;

    G4ThreeVector positionSiLi_RB_down = G4ThreeVector(   0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
        -0.5 * SiLi_HighY_Upper  - SiLi_HighY_Center - 1.5 * interSiLi,
        0);

    positionSiLi_RB_down += ShiftSiLiDown ;

    G4ThreeVector positionSiLi_LC2_down = G4ThreeVector(  -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
        -0.5 * SiLi_HighY_Center - 0.5 * interSiLi,
        0);

    positionSiLi_LC2_down += ShiftSiLiDown ;

    G4ThreeVector positionSiLi_RC2_down = G4ThreeVector(  0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
        -0.5 * SiLi_HighY_Center - 0.5 * interSiLi,
        0);

    positionSiLi_RC2_down += ShiftSiLiDown ;


    // up
    new G4PVPlacement(0 , positionSiLi_LT_up  , logicSiLi_LT  , Name + "_SiLi_Pad9"   , logicSiLi , false , 0)  ;
    new G4PVPlacement(0 , positionSiLi_RT_up  , logicSiLi_RT  , Name + "_SiLi_Pad10"  , logicSiLi , false , 0)  ;
    new G4PVPlacement(0 , positionSiLi_LC1_up , logicSiLi_LC1 , Name + "_SiLi_Pad11" 	, logicSiLi , false , 0);
    new G4PVPlacement(0 , positionSiLi_RC1_up , logicSiLi_RC1 , Name + "_SiLi_Pad12" 	, logicSiLi , false , 0);

    new G4PVPlacement(0 , positionSiLi_LB_up  , logicSiLi_LB  , Name + "_SiLi_Pad16"  , logicSiLi , false , 0)  ;
    new G4PVPlacement(0 , positionSiLi_RB_up  , logicSiLi_RB  , Name + "_SiLi_Pad15"  , logicSiLi , false , 0)  ;
    new G4PVPlacement(0 , positionSiLi_LC2_up , logicSiLi_LC2 , Name + "_SiLi_Pad14" 	, logicSiLi , false , 0);
    new G4PVPlacement(0 , positionSiLi_RC2_up , logicSiLi_RC2 , Name + "_SiLi_Pad13" 	, logicSiLi , false , 0);


    // down
    new G4PVPlacement(0 , positionSiLi_LT_down  , logicSiLi_LT  , Name + "_SiLi_Pad2"  , logicSiLi , false , 0) ;
    new G4PVPlacement(0 , positionSiLi_RT_down  , logicSiLi_RT  , Name + "_SiLi_Pad1"  , logicSiLi , false , 0) ;
    new G4PVPlacement(0 , positionSiLi_LC1_down , logicSiLi_LC1 , Name + "_SiLi_Pad4" 	, logicSiLi , false , 0) ;
    new G4PVPlacement(0 , positionSiLi_RC1_down , logicSiLi_RC1 , Name + "_SiLi_Pad3" 	, logicSiLi , false , 0) ;

    new G4PVPlacement(0 , positionSiLi_LB_down  , logicSiLi_LB  , Name + "_SiLi_Pad7"  , logicSiLi , false , 0) ;
    new G4PVPlacement(0 , positionSiLi_RB_down  , logicSiLi_RB  , Name + "_SiLi_Pad8"  , logicSiLi , false , 0) ;
    new G4PVPlacement(0 , positionSiLi_LC2_down , logicSiLi_LC2 , Name + "_SiLi_Pad5" 	, logicSiLi , false , 0) ;
    new G4PVPlacement(0 , positionSiLi_RC2_down , logicSiLi_RC2 , Name + "_SiLi_Pad6" 	, logicSiLi , false , 0) ;

    // Set SiLi sensible
    logicSiLi_LT->SetSensitiveDetector(m_SiLiScorer);
    logicSiLi_RT->SetSensitiveDetector(m_SiLiScorer);
    logicSiLi_LC1->SetSensitiveDetector(m_SiLiScorer);
    logicSiLi_RC1->SetSensitiveDetector(m_SiLiScorer);

    logicSiLi_LB->SetSensitiveDetector(m_SiLiScorer);
    logicSiLi_RB->SetSensitiveDetector(m_SiLiScorer);
    logicSiLi_LC2->SetSensitiveDetector(m_SiLiScorer);
    logicSiLi_RC2->SetSensitiveDetector(m_SiLiScorer);
    // Mark blue a SiLi to see telescope orientation

    logicSiLi_LT->SetVisAttributes(m_SiLiVisAtt);
    logicSiLi_RT->SetVisAttributes(m_SiLiVisAtt);
    logicSiLi_LC1->SetVisAttributes(m_SiLiVisAtt);
    logicSiLi_RC1->SetVisAttributes(m_SiLiVisAtt);

    logicSiLi_LB->SetVisAttributes(m_SiLiVisAtt);
    logicSiLi_RB->SetVisAttributes(m_SiLiVisAtt);
    logicSiLi_LC2->SetVisAttributes(m_SiLiVisAtt);
    logicSiLi_RC2->SetVisAttributes(m_SiLiVisAtt);


    delete rotSiLi;
  }




}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MUST2Array::VolumeMaker( G4int TelescopeNumber,
    G4ThreeVector MMpos,
    G4RotationMatrix* MMrot,
    bool wSi,
    bool wSiLi,
    bool wCsI,
    G4LogicalVolume* world)
{
  G4double NbrTelescopes = TelescopeNumber;
  G4String DetectorNumber;
  std::ostringstream Number;
  Number << NbrTelescopes ;
  DetectorNumber = Number.str();


   ////////////////////////////////////////////////////////////////
  ////////////// Starting Volume Definition //////////////////////
  ////////////////////////////////////////////////////////////////
  G4Trd*           solidMM = new G4Trd("MUST2Telescope" + DetectorNumber, 0.5*FaceFront, 0.5*FaceBack, 0.5*FaceFront, 0.5*FaceBack, 0.5*Length);
  G4LogicalVolume* logicMM = new G4LogicalVolume(solidMM, m_MaterialIron, "MUST2Telescope" + DetectorNumber, 0, 0, 0);
  G4String Name = "MUST2Telescope" + DetectorNumber ;

  new G4PVPlacement(	G4Transform3D(*MMrot, MMpos),
      logicMM               ,
      Name                  ,
      world                 ,
      false                 ,
      0                     );

  if (m_non_sensitive_part_visiualisation){
    G4VisAttributes* FrameVisAtt = new G4VisAttributes(G4Colour(0.80, 0.80, 0.80));
    FrameVisAtt->SetForceWireframe(true);
    logicMM->SetVisAttributes(FrameVisAtt)  ;
  }

  else logicMM->SetVisAttributes(G4VisAttributes::Invisible)               ;

  G4ThreeVector positionVacBox = G4ThreeVector(0, 0, VacBox_PosZ);

  G4Trd* solidVacBox = 
    new G4Trd("solidVacBox", 0.5*ActiveWaferSize, 0.5*CsIFaceFront, 0.5*ActiveWaferSize, 0.5*CsIFaceFront, 0.5*VacBoxThickness);

  G4LogicalVolume* logicVacBox = new G4LogicalVolume(solidVacBox, m_MaterialVacuum, "logicVacBox", 0, 0, 0);


  new G4PVPlacement(0, positionVacBox, logicVacBox, Name + "_VacBox", logicMM, false, 0);

  logicVacBox->SetVisAttributes(G4VisAttributes::Invisible);


  ////////////////////////////////////////////////////////////////
  /////////////////Si Strip Construction//////////////////////////
  ////////////////////////////////////////////////////////////////
  // A no rotation matrix is always handy ;) 
  G4RotationMatrix* norotation = new  G4RotationMatrix();   

  // Create the front support frame for Si and SiLi
  vector<G4TwoVector> polygon;
  for(unsigned int i = 0 ; i < 4 ; i++){   
    G4TwoVector Point(PCBPointsX[i],PCBPointsY[i]);   
    polygon.push_back(Point);  
  } 

  // The PCB lenth is used to rescale different part of the detector
  G4double PCBlength = PCBPointsX[0] - PCBPointsX[1];

  G4ExtrudedSolid* solidFrontFrameBase = new G4ExtrudedSolid("Must2 Front Frame",
      polygon,   
      FrontFrameThickness*0.5,
      G4TwoVector(0,0),FrontFrameFrontSize/PCBlength,
      G4TwoVector(0,0),FrontFrameBackSize/PCBlength); 

  G4ExtrudedSolid* solidCavityShape1 = new G4ExtrudedSolid("Cavity1",
      polygon,   
      FrontFrameThickness*0.5+PCBThickness,
      G4TwoVector(0,0),CavityShape1FrontSize/PCBlength,
      G4TwoVector(0,0),CavityShape1BackSize/PCBlength); 

  G4SubtractionSolid* solidFrontFrame1 = new G4SubtractionSolid("Frame1",
      solidFrontFrameBase,
      solidCavityShape1,
      G4Transform3D(*norotation,G4ThreeVector(0,0,0))); 

  G4Box* solidCavityShape2 = new G4Box("Cavity2", 0.5*CavityShape2Width, (0.5+0.5)*CavityShape2Width, 0.5*CavityShape2Thickness);

  G4SubtractionSolid* solidFrontFrame2 = new G4SubtractionSolid("Frame2",
      solidFrontFrame1,
      solidCavityShape2,
      G4Transform3D(*norotation,G4ThreeVector(0,0.5*CavityShape2Width,CavityShape2Center))); 

  G4ExtrudedSolid* solidCutOut = new G4ExtrudedSolid("Cutout",
      polygon,   
      FrontFrameThickness*0.8,
      G4TwoVector(0,0),1.1*FrontFrameBackSize/PCBlength,
      G4TwoVector(0,0),1.1*FrontFrameBackSize/PCBlength); 


  G4SubtractionSolid* solidFrontFrame3 = new G4SubtractionSolid("Must2 Front Frame",
      solidFrontFrame2,
      solidCutOut,
      G4Transform3D(*norotation,G4ThreeVector(0,FrontFrameBackSize*0.5*1.1+66*mm,0))); 

  G4SubtractionSolid* solidFrontFrame = new G4SubtractionSolid("Must2 Front Frame",
      solidFrontFrame3,
      solidCutOut,
      G4Transform3D(*norotation,G4ThreeVector(-FrontFrameBackSize*0.5*1.1-66*mm,0,0))); 


  G4LogicalVolume* logicFrontFrame = new G4LogicalVolume(solidFrontFrame,m_MaterialAluminium,"Must2 Front Frame", 0, 0, 0); 

  new G4PVPlacement(G4Transform3D(*norotation,G4ThreeVector(0,0,Silicon_PosZ+FrontFrameThickness*0.5+0.5*PCBThickness)),
      logicFrontFrame, 
      "Must2 Front Fame",
      logicMM,
      false,
      TelescopeNumber); 

  // Electronic bloc
  G4Box* solidMUFEE  = new G4Box("solidMUFEE", 0.5*MUFEESize, 0.5*MUFEESize, 0.5*MUFEEThickness);
  G4Box* solidCooling = new G4Box("solidCooling", 0.5*MUFEESize, 0.5*MUFEESize, 0.5*CoolingThickness);
  G4ExtrudedSolid* solidElecFrameShort = new G4ExtrudedSolid("SolidElecFrame",
      polygon,   
      CoolingThickness*0.5,
      G4TwoVector(0,0),FrontFrameBackSize/PCBlength,
      G4TwoVector(0,0),FrontFrameBackSize/PCBlength); 


  G4LogicalVolume* logicMUFEE   = new G4LogicalVolume(solidMUFEE,m_MaterialAluminium,"Must2 MUFEE", 0, 0, 0); 
  G4LogicalVolume* logicCooling = new G4LogicalVolume(solidCooling,m_MaterialAluminium,"Must2 Cooling", 0, 0, 0); 
  G4LogicalVolume* logicElecFrame = new G4LogicalVolume(solidElecFrameShort,m_MaterialAluminium,"Must2 Frame", 0, 0, 0); 

  logicMUFEE->SetVisAttributes(m_PCBVisAtt);
  logicCooling->SetVisAttributes(m_CoolingVisAtt);

  G4double CoolingPosition = FrontFrameThickness+Silicon_PosZ+0.5*CoolingThickness;

  new G4PVPlacement(G4Transform3D(*norotation,G4ThreeVector(0,0,CoolingPosition)),
      logicCooling, 
      "Must2 Cooling",
      logicMM,
      false,
      TelescopeNumber); 

  new G4PVPlacement(G4Transform3D(*norotation,G4ThreeVector(0,0,CoolingPosition)),
      logicElecFrame, 
      "Must2 Fame",
      logicMM,
      false,
      TelescopeNumber); 


  new G4PVPlacement(G4Transform3D(*norotation,G4ThreeVector(0,0,CoolingPosition-0.5*MUFEEThickness-0.5*CoolingThickness)),
      logicMUFEE, 
      "Must2 Front Fame",
      logicMM,
      false,
      TelescopeNumber); 

  new G4PVPlacement(G4Transform3D(*norotation,G4ThreeVector(0,0,CoolingPosition+0.5*MUFEEThickness+0.5*CoolingThickness)),
      logicMUFEE, 
      "Must2 Front Fame",
      logicMM,
      false,
      TelescopeNumber); 




  if (wSi) {

    // Volume containing the Silicon detector
    G4ExtrudedSolid* solidSilicon = new G4ExtrudedSolid(Name,
        polygon,   
        PCBThickness*0.5,
        G4TwoVector(0,0),1,
        G4TwoVector(0,0),1); 

    G4LogicalVolume* logicSilicon = new G4LogicalVolume(solidSilicon,m_MaterialVacuum, Name, 0, 0, 0); 

    new G4PVPlacement(G4Transform3D(*norotation,G4ThreeVector(0,0,Silicon_PosZ)),
        logicSilicon, 
        "Silicon",
        logicMM,
        false,
        TelescopeNumber); 

    logicSilicon->SetVisAttributes(G4VisAttributes::Invisible);

    // PCB Base
    G4ExtrudedSolid* solidPCBBase = new G4ExtrudedSolid("PCBBase",
        polygon,   
        PCBThickness*0.5,
        G4TwoVector(0,0),1,
        G4TwoVector(0,0),1); 

    // Wafer Shape 
    G4Box* solidWaferShape = new G4Box("solidWaferShape", 0.5*WaferSize, 0.5*WaferSize, 0.6*PCBThickness);

    // PCB

    G4SubtractionSolid* solidPCB = new G4SubtractionSolid("PCB",
        solidPCBBase,
        solidWaferShape,
        G4Transform3D(*norotation,G4ThreeVector(0,0,0))); 

    G4LogicalVolume* logicPCB = new G4LogicalVolume(solidPCB, m_MaterialSilicon, "logicPCB", 0, 0, 0);

    new G4PVPlacement(0, G4ThreeVector(0,0,0), logicPCB, Name + "_Silicon", logicSilicon, false, 0);

    logicPCB->SetVisAttributes(m_PCBVisAtt); 


    // Wafer
    G4Box* solidWafer = new G4Box("solidWaferShape", 0.5*WaferSize, 0.5*WaferSize, 0.5*WaferThickness);
    G4LogicalVolume* logicWafer = new G4LogicalVolume(solidWafer, m_MaterialSilicon, "logicWafer", 0, 0, 0);
    new G4PVPlacement(0, G4ThreeVector(0,0,0), logicWafer, Name + "_Silicon", logicSilicon, false, 0);


    // Active Wafer
    G4Box* solidActiveWafer = new G4Box("solidActiveWaferShape", 0.5*ActiveWaferSize, 0.5*ActiveWaferSize, 0.5*WaferThickness);
    G4LogicalVolume* logicActiveWafer = new G4LogicalVolume(solidActiveWafer, m_MaterialSilicon, "logicActiveWafer", 0, 0, 0);

    new G4PVPlacement(0, G4ThreeVector(0,0,0), logicActiveWafer, Name + "_Silicon", logicSilicon, false, 0);

    G4ThreeVector positionAluStripFront = G4ThreeVector(0, 0, 0.5*WaferThickness+0.5*AluStripThickness);
    G4ThreeVector positionAluStripBack  = G4ThreeVector(0, 0, -0.5*WaferThickness-0.5*AluStripThickness);

    G4Box*           solidAluStrip = new G4Box("AluBox", 0.5*ActiveWaferSize, 0.5*ActiveWaferSize, 0.5*AluStripThickness);
    G4LogicalVolume* logicAluStrip = new G4LogicalVolume(solidAluStrip, m_MaterialAluminium, "logicAluStrip", 0, 0, 0);

    new G4PVPlacement(0, positionAluStripFront, logicAluStrip, "AluStripFront", logicSilicon, false, 0);
    new G4PVPlacement(0, positionAluStripBack, logicAluStrip, "AluStripBack", logicSilicon, false, 0);

    logicAluStrip->SetVisAttributes(G4VisAttributes::Invisible);

    G4ThreeVector  positionSilicon = G4ThreeVector(0, 0, Silicon_PosZ);

    new G4PVPlacement(0, G4ThreeVector(0, 0, Silicon_PosZ), logicSilicon, "Silicon", logicMM, false, 0);

    // Set Silicon strip sensible
    logicActiveWafer->SetSensitiveDetector(m_StripScorer);

    // Visualisation of Silicon Strip
    logicWafer->SetVisAttributes(m_SiliconVisAtt);
    logicActiveWafer->SetVisAttributes(m_SiliconVisAtt);
  }

  ////////////////////////////////////////////////////////////////
  //////////////////// SiLi  Construction ////////////////////////
  ////////////////////////////////////////////////////////////////

  if (wSiLi) {
    G4double SiLiSpace = 8 * mm;
    G4RotationMatrix* rotSiLi = new G4RotationMatrix(0,0,0); 
    G4Box* solidSiLi = new G4Box("SiLi", 0.5*ActiveWaferSize+0.5*SiLiSpace, 0.5*ActiveWaferSize, 0.5*SiLiThickness);
    G4LogicalVolume* logicSiLi = new G4LogicalVolume(solidSiLi, m_MaterialAluminium, Name + "_SiLi" , 0, 0, 0);
    logicSiLi->SetVisAttributes(G4VisAttributes::Invisible);

    new G4PVPlacement(  G4Transform3D(*rotSiLi, G4ThreeVector(0,0,0) ) 	,
        logicSiLi ,
        Name + "_SiLi",
        logicVacBox ,
        false ,
        0);

    // SiLi are placed inside of the VacBox...
    // Left/Right define when looking to detector from Si to CsI
    G4double SiLi_HighY_Upper 	= 19.86 * mm;
    G4double SiLi_HighY_Center = 25.39 * mm;
    G4double SiLi_WidthX_Left  = 22.85 * mm;
    G4double SiLi_WidthX_Right = 24.9  * mm;
    G4double SiLi_ShiftX    	= 0.775 * mm;

    //	SiLi are organized by two group of 8 Up(9 to 15) and Down(1 to 8).
    G4ThreeVector ShiftSiLiUp 	 = G4ThreeVector(-0.25 * ActiveWaferSize - 0.5 * SiLiSpace, 0, 0)	;
    G4ThreeVector ShiftSiLiDown  = G4ThreeVector(0.25  * ActiveWaferSize + 0.5 * SiLiSpace, 0, 0)	;

    // SiLi : left side of SiLi detector
    G4Box* solidSiLi_LT  = new G4Box("SiLi_LT"  , 0.5*SiLi_WidthX_Left  , 0.5*SiLi_HighY_Upper   , 0.5*SiLiThickness);
    G4Box* solidSiLi_RT  = new G4Box("SiLi_RT"  , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Upper   , 0.5*SiLiThickness);
    G4Box* solidSiLi_LC1 = new G4Box("SiLi_LC1" , 0.5*SiLi_WidthX_Left  , 0.5*SiLi_HighY_Center  , 0.5*SiLiThickness);
    G4Box* solidSiLi_RC1 = new G4Box("SiLi_RC1" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Center  , 0.5*SiLiThickness);
    G4Box* solidSiLi_LB  = new G4Box("SiLi_LB"  , 0.5*SiLi_WidthX_Left  , 0.5*SiLi_HighY_Upper   , 0.5*SiLiThickness);
    G4Box* solidSiLi_RB  = new G4Box("SiLi_RB"  , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Upper   , 0.5*SiLiThickness);
    G4Box* solidSiLi_LC2 = new G4Box("SiLi_LC2" , 0.5*SiLi_WidthX_Left  , 0.5*SiLi_HighY_Center  , 0.5*SiLiThickness);
    G4Box* solidSiLi_RC2 = new G4Box("SiLi_RC2" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Center  , 0.5*SiLiThickness);

    G4LogicalVolume* logicSiLi_LT    = new G4LogicalVolume(solidSiLi_LT   , m_MaterialSilicon , "SiLi_LT"  , 0 , 0 , 0);
    G4LogicalVolume* logicSiLi_RT    = new G4LogicalVolume(solidSiLi_RT   , m_MaterialSilicon , "SiLi_RT"  , 0 , 0 , 0);
    G4LogicalVolume* logicSiLi_LC1   = new G4LogicalVolume(solidSiLi_LC1  , m_MaterialSilicon , "SiLi_LC1" , 0 , 0 , 0);
    G4LogicalVolume* logicSiLi_RC1   = new G4LogicalVolume(solidSiLi_RC1  , m_MaterialSilicon , "SiLi_RC1" , 0 , 0 , 0);
    G4LogicalVolume* logicSiLi_LB    = new G4LogicalVolume(solidSiLi_LB   , m_MaterialSilicon , "SiLi_LB"  , 0 , 0 , 0);
    G4LogicalVolume* logicSiLi_RB    = new G4LogicalVolume(solidSiLi_RB   , m_MaterialSilicon , "SiLi_RB"  , 0 , 0 , 0);
    G4LogicalVolume* logicSiLi_LC2   = new G4LogicalVolume(solidSiLi_LC2  , m_MaterialSilicon , "SiLi_LC2" , 0 , 0 , 0);
    G4LogicalVolume* logicSiLi_RC2   = new G4LogicalVolume(solidSiLi_RC2  , m_MaterialSilicon , "SiLi_RC2" , 0 , 0 , 0);

    G4double interSiLi = 0.5 * mm;

    // Top
    G4ThreeVector positionSiLi_LT_up = G4ThreeVector(  -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
        0.5 * SiLi_HighY_Upper  + SiLi_HighY_Center + 1.5 * interSiLi,
        0);

    positionSiLi_LT_up += ShiftSiLiUp ;

    G4ThreeVector positionSiLi_RT_up = G4ThreeVector(  0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
        0.5 * SiLi_HighY_Upper  + SiLi_HighY_Center + 1.5 * interSiLi,
        0);

    positionSiLi_RT_up += ShiftSiLiUp ;

    G4ThreeVector positionSiLi_LC1_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
        0.5 * SiLi_HighY_Center + 0.5 * interSiLi,
        0);

    positionSiLi_LC1_up += ShiftSiLiUp ;

    G4ThreeVector positionSiLi_RC1_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
        0.5 * SiLi_HighY_Center + 0.5 * interSiLi,
        0);

    positionSiLi_RC1_up += ShiftSiLiUp ;

    G4ThreeVector positionSiLi_LB_up = G4ThreeVector(  -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
        -0.5 * SiLi_HighY_Upper  - SiLi_HighY_Center - 1.5 * interSiLi ,
        0);

    positionSiLi_LB_up += ShiftSiLiUp ;

    G4ThreeVector positionSiLi_RB_up = G4ThreeVector(  0.5 * SiLi_WidthX_Right - SiLi_ShiftX ,
        -0.5 * SiLi_HighY_Upper  - SiLi_HighY_Center - 1.5 * interSiLi ,
        0);

    positionSiLi_RB_up += ShiftSiLiUp ;

    G4ThreeVector positionSiLi_LC2_up = G4ThreeVector( -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX ,
        -0.5 * SiLi_HighY_Center - 0.5 * interSiLi,
        0);

    positionSiLi_LC2_up += ShiftSiLiUp ;

    G4ThreeVector positionSiLi_RC2_up = G4ThreeVector( 0.5 * SiLi_WidthX_Right - SiLi_ShiftX ,
        -0.5 * SiLi_HighY_Center - 0.5 * interSiLi  ,
        0);

    positionSiLi_RC2_up += ShiftSiLiUp ;


    // Down
    G4ThreeVector positionSiLi_LT_down = G4ThreeVector(   -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
        0.5 * SiLi_HighY_Upper  + SiLi_HighY_Center + 1.5 * interSiLi,
        0);

    positionSiLi_LT_down += ShiftSiLiDown ;

    G4ThreeVector positionSiLi_RT_down = G4ThreeVector(   0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
        0.5 * SiLi_HighY_Upper  + SiLi_HighY_Center + 1.5 * interSiLi,
        0);

    positionSiLi_RT_down += ShiftSiLiDown ;

    G4ThreeVector positionSiLi_LC1_down = G4ThreeVector(  -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
        0.5 * SiLi_HighY_Center + 0.5 * interSiLi,
        0);

    positionSiLi_LC1_down += ShiftSiLiDown ;

    G4ThreeVector positionSiLi_RC1_down = G4ThreeVector(  0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
        0.5 * SiLi_HighY_Center + 0.5 * interSiLi ,
        0);

    positionSiLi_RC1_down += ShiftSiLiDown ;

    G4ThreeVector positionSiLi_LB_down = G4ThreeVector(   -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
        -0.5 * SiLi_HighY_Upper  - SiLi_HighY_Center - 1.5 * interSiLi,
        0);

    positionSiLi_LB_down += ShiftSiLiDown ;

    G4ThreeVector positionSiLi_RB_down = G4ThreeVector(   0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
        -0.5 * SiLi_HighY_Upper  - SiLi_HighY_Center - 1.5 * interSiLi,
        0);

    positionSiLi_RB_down += ShiftSiLiDown ;

    G4ThreeVector positionSiLi_LC2_down = G4ThreeVector(  -0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX,
        -0.5 * SiLi_HighY_Center - 0.5 * interSiLi,
        0);

    positionSiLi_LC2_down += ShiftSiLiDown ;

    G4ThreeVector positionSiLi_RC2_down = G4ThreeVector(  0.5 * SiLi_WidthX_Right - SiLi_ShiftX,
        -0.5 * SiLi_HighY_Center - 0.5 * interSiLi,
        0);

    positionSiLi_RC2_down += ShiftSiLiDown ;


    // up
    new G4PVPlacement(0 , positionSiLi_LT_up  , logicSiLi_LT  , Name + "_SiLi_Pad9"   , logicSiLi , false , 0)  ;
    new G4PVPlacement(0 , positionSiLi_RT_up  , logicSiLi_RT  , Name + "_SiLi_Pad10"  , logicSiLi , false , 0)  ;
    new G4PVPlacement(0 , positionSiLi_LC1_up , logicSiLi_LC1 , Name + "_SiLi_Pad11" 	, logicSiLi , false , 0);
    new G4PVPlacement(0 , positionSiLi_RC1_up , logicSiLi_RC1 , Name + "_SiLi_Pad12" 	, logicSiLi , false , 0);

    new G4PVPlacement(0 , positionSiLi_LB_up  , logicSiLi_LB  , Name + "_SiLi_Pad16"  , logicSiLi , false , 0)  ;
    new G4PVPlacement(0 , positionSiLi_RB_up  , logicSiLi_RB  , Name + "_SiLi_Pad15"  , logicSiLi , false , 0)  ;
    new G4PVPlacement(0 , positionSiLi_LC2_up , logicSiLi_LC2 , Name + "_SiLi_Pad14" 	, logicSiLi , false , 0);
    new G4PVPlacement(0 , positionSiLi_RC2_up , logicSiLi_RC2 , Name + "_SiLi_Pad13" 	, logicSiLi , false , 0);


    // down
    new G4PVPlacement(0 , positionSiLi_LT_down  , logicSiLi_LT  , Name + "_SiLi_Pad2"  , logicSiLi , false , 0) ;
    new G4PVPlacement(0 , positionSiLi_RT_down  , logicSiLi_RT  , Name + "_SiLi_Pad1"  , logicSiLi , false , 0) ;
    new G4PVPlacement(0 , positionSiLi_LC1_down , logicSiLi_LC1 , Name + "_SiLi_Pad4" 	, logicSiLi , false , 0) ;
    new G4PVPlacement(0 , positionSiLi_RC1_down , logicSiLi_RC1 , Name + "_SiLi_Pad3" 	, logicSiLi , false , 0) ;

    new G4PVPlacement(0 , positionSiLi_LB_down  , logicSiLi_LB  , Name + "_SiLi_Pad7"  , logicSiLi , false , 0) ;
    new G4PVPlacement(0 , positionSiLi_RB_down  , logicSiLi_RB  , Name + "_SiLi_Pad8"  , logicSiLi , false , 0) ;
    new G4PVPlacement(0 , positionSiLi_LC2_down , logicSiLi_LC2 , Name + "_SiLi_Pad5" 	, logicSiLi , false , 0) ;
    new G4PVPlacement(0 , positionSiLi_RC2_down , logicSiLi_RC2 , Name + "_SiLi_Pad6" 	, logicSiLi , false , 0) ;

    // Set SiLi sensible
    logicSiLi_LT->SetSensitiveDetector(m_SiLiScorer);
    logicSiLi_RT->SetSensitiveDetector(m_SiLiScorer);
    logicSiLi_LC1->SetSensitiveDetector(m_SiLiScorer);
    logicSiLi_RC1->SetSensitiveDetector(m_SiLiScorer);

    logicSiLi_LB->SetSensitiveDetector(m_SiLiScorer);
    logicSiLi_RB->SetSensitiveDetector(m_SiLiScorer);
    logicSiLi_LC2->SetSensitiveDetector(m_SiLiScorer);
    logicSiLi_RC2->SetSensitiveDetector(m_SiLiScorer);

    // Mark blue a SiLi to see telescope orientation

    logicSiLi_LT->SetVisAttributes(m_SiLiVisAtt);
    logicSiLi_RT->SetVisAttributes(m_SiLiVisAtt);
    logicSiLi_LC1->SetVisAttributes(m_SiLiVisAtt);
    logicSiLi_RC1->SetVisAttributes(m_SiLiVisAtt);

    logicSiLi_LB->SetVisAttributes(m_SiLiVisAtt);
    logicSiLi_RB->SetVisAttributes(m_SiLiVisAtt);
    logicSiLi_LC2->SetVisAttributes(m_SiLiVisAtt);
    logicSiLi_RC2->SetVisAttributes(m_SiLiVisAtt);


    delete rotSiLi;
  }

  ////////////////////////////////////////////////////////////////
  //////////////////// CsI  Construction//////////////////////////
  ////////////////////////////////////////////////////////////////

  if (wCsI) {

    G4ThreeVector positionCsI = G4ThreeVector(0, 0, CsI_PosZ);
    G4Trd* solidCsI = new G4Trd("csI", 0.5*CsIFaceFront, 0.5*CsIFaceBack, 0.5*CsIFaceFront, 0.5*CsIFaceBack, 0.5*CsIThickness);

    G4LogicalVolume* logicCsI = new G4LogicalVolume(solidCsI, m_MaterialAluminium, Name + "_CsI_Mylar", 0, 0, 0);
    new G4PVPlacement(0, positionCsI, logicCsI, Name + "_CsI_Mylar", logicMM, false, 0);

    G4ThreeVector   positionMylarCsI = G4ThreeVector(0, 0, MylarCsIThickness * 0.5 - CsIThickness * 0.5);

    G4Box*           solidMylarCsI = new G4Box("MylarCsIBox", 0.5*CsIFaceFront, 0.5*CsIFaceFront, 0.5*MylarCsIThickness);
    G4LogicalVolume* logicMylarCsI = new G4LogicalVolume(solidMylarCsI, m_MaterialMyl, Name + "_CsI_Mylar", 0, 0, 0);

    new G4PVPlacement(0, positionMylarCsI, logicMylarCsI, Name + "_CsI_Mylar", logicCsI, false, 0);


    logicCsI->SetVisAttributes(G4VisAttributes::Invisible);
    logicMylarCsI->SetVisAttributes(G4VisAttributes::Invisible);

    // Cristal1
    G4Trap* solidCristal1 = new G4Trap("Cristal1", 40.*mm / 2., 6.693896*deg, 41.97814*deg, 33.1*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 26.9*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg);
    G4LogicalVolume* logicCristal1 = new G4LogicalVolume(solidCristal1, m_MaterialCsI, Name + "_CsI_Cristal1", 0, 0, 0);

    // Cristal2
    G4Trap* solidCristal2 = new G4Trap("Cristal2", 40.*mm / 2., 17.8836*deg, (74.3122 + 180)*deg, 43.49*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 31.0377*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg);
    G4LogicalVolume* logicCristal2 = new G4LogicalVolume(solidCristal2, m_MaterialCsI, Name + "_CsI_Cristal2", 0, 0, 0);

    // Cristal3
    G4Trap* solidCristal3 = new G4Trap("Cristal3", 40.*mm / 2., 18.243*deg, 13.5988*deg, 33.11*mm / 2., 39.25*mm / 2., 39.25*mm / 2., 0.*deg, 26.91*mm / 2., 27.58*mm / 2., 27.58*mm / 2., 0.*deg);
    G4LogicalVolume* logicCristal3 = new G4LogicalVolume(solidCristal3, m_MaterialCsI, Name + "_CsI_Cristal3", 0, 0, 0);

    // Cristal4

    G4Trap* solidCristal4 = new G4Trap("Cristal4", 40.*mm / 2., 24.0482*deg, 44.1148*deg, 43.49*mm / 2., 39.19*mm / 2., 39.19*mm / 2., 0.*deg, 31.04*mm / 2., 27.52*mm / 2., 27.52*mm / 2., 0.*deg);
    G4LogicalVolume* logicCristal4 = new G4LogicalVolume(solidCristal4, m_MaterialCsI, Name + "_CsI_Cristal4", 0, 0, 0);


    // Cristal1s

    G4Trap* solidCristal1s = new G4Trap("Cristal1s", 40.*mm / 2., 6.693896*deg, -41.97814*deg, 33.1*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 26.9*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg);
    G4LogicalVolume* logicCristal1s = new G4LogicalVolume(solidCristal1s, m_MaterialCsI, Name + "_CsI_Cristal1s", 0, 0, 0);

    // Cristal2s

    G4Trap* solidCristal2s = new G4Trap("Cristal2s", 40.*mm / 2., 17.8836*deg, -(74.3122 + 180)*deg, 43.49*mm / 2., 37.39*mm / 2., 37.39*mm / 2., 0.*deg, 31.0377*mm / 2., 30.41*mm / 2., 30.41*mm / 2., 0.*deg);
    G4LogicalVolume* logicCristal2s = new G4LogicalVolume(solidCristal2s, m_MaterialCsI, Name + "_CsI_Cristal2s", 0, 0, 0);

    // Cristal3s

    G4Trap* solidCristal3s = new G4Trap("Cristal3s", 40.*mm / 2., 18.243*deg, -13.5988*deg, 33.11*mm / 2., 39.25*mm / 2., 39.25*mm / 2., 0.*deg, 26.91*mm / 2., 27.58*mm / 2., 27.58*mm / 2., 0.*deg);
    G4LogicalVolume* logicCristal3s = new G4LogicalVolume(solidCristal3s, m_MaterialCsI, Name + "_CsI_Cristal3s", 0, 0, 0);

    // Cristal4s

    G4Trap* solidCristal4s = new G4Trap("Cristal4s", 40.*mm / 2., 24.0482*deg, -44.1148*deg, 43.49*mm / 2., 39.19*mm / 2., 39.19*mm / 2., 0.*deg, 31.04*mm / 2., 27.52*mm / 2., 27.52*mm / 2., 0.*deg);
    G4LogicalVolume* logicCristal4s = new G4LogicalVolume(solidCristal4s, m_MaterialCsI, Name + "_CsI_Cristal4s", 0, 0, 0);

    G4double XEdge1 = 16.96 * mm + DistInterCsI * 0.5;
    G4double YEdge1 = 15.01 * mm + DistInterCsI * 0.5;
    G4double XEdge2 = 50.63 * mm + DistInterCsI * 1.5;
    G4double YEdge2 = 48.67 * mm + DistInterCsI * 1.5;

    G4ThreeVector positionCristal1 = G4ThreeVector(-XEdge1, YEdge1, 0 * mm);
    G4ThreeVector positionCristal2 = G4ThreeVector(-XEdge1, YEdge2, 0);
    G4ThreeVector positionCristal3 = G4ThreeVector(-XEdge2, YEdge1, 0);
    G4ThreeVector positionCristal4 = G4ThreeVector(-XEdge2, YEdge2, 0);

    new G4PVPlacement(Rotation(180., 0., 0.), positionCristal1, logicCristal1, Name + "_CsI_Cristal1", logicCsI, false, 1);
    new G4PVPlacement(Rotation(180., 0., 180.), positionCristal2, logicCristal2, Name + "_CsI_Cristal2", logicCsI, false, 2);
    new G4PVPlacement(Rotation(180., 0., 0.), positionCristal3, logicCristal3, Name + "_CsI_Cristal3", logicCsI, false, 3);
    new G4PVPlacement(Rotation(180., 0., 0.), positionCristal4, logicCristal4, Name + "_CsI_Cristal4", logicCsI, false, 4);


    G4ThreeVector positionCristal1b = G4ThreeVector(XEdge1, -YEdge1, 0 * mm);
    G4ThreeVector positionCristal2b = G4ThreeVector(XEdge1, -YEdge2, 0);
    G4ThreeVector positionCristal3b = G4ThreeVector(XEdge2, -YEdge1, 0);
    G4ThreeVector positionCristal4b = G4ThreeVector(XEdge2, -YEdge2, 0);

    new G4PVPlacement(Rotation(180., 0., 180.), positionCristal1b, logicCristal1, Name + "_CsI_Cristal5", logicCsI, false, 5);
    new G4PVPlacement(Rotation(180., 0., 0.), positionCristal2b, logicCristal2, Name + "_CsI_Cristal6", logicCsI, false, 6);
    new G4PVPlacement(Rotation(180., 0., 180.), positionCristal3b, logicCristal3, Name + "_CsI_Cristal7", logicCsI, false, 7);
    new G4PVPlacement(Rotation(180., 0., 180.), positionCristal4b, logicCristal4, Name + "_CsI_Cristal8", logicCsI, false, 8);

    G4ThreeVector positionCristal1s = G4ThreeVector(-XEdge1, -YEdge1, 0 * mm);
    G4ThreeVector positionCristal2s = G4ThreeVector(-XEdge1, -YEdge2, 0);
    G4ThreeVector positionCristal3s = G4ThreeVector(-XEdge2, -YEdge1, 0);
    G4ThreeVector positionCristal4s = G4ThreeVector(-XEdge2, -YEdge2, 0);

    new G4PVPlacement(Rotation(180., 0., 0.), positionCristal1s, logicCristal1s, Name + "_CsI_Cristal9", logicCsI, false, 9);
    new G4PVPlacement(Rotation(180., 0., 180.), positionCristal2s, logicCristal2s, Name + "_CsI_Cristal10", logicCsI, false, 10);
    new G4PVPlacement(Rotation(180., 0., 0.), positionCristal3s, logicCristal3s, Name + "_CsI_Cristal11", logicCsI, false, 11);
    new G4PVPlacement(Rotation(180., 0., 0.), positionCristal4s, logicCristal4s, Name + "_CsI_Cristal12", logicCsI, false, 12);

    G4ThreeVector positionCristal1sb = G4ThreeVector(XEdge1, YEdge1, 0 * mm);
    G4ThreeVector positionCristal2sb = G4ThreeVector(XEdge1, YEdge2, 0);
    G4ThreeVector positionCristal3sb = G4ThreeVector(XEdge2, YEdge1, 0);
    G4ThreeVector positionCristal4sb = G4ThreeVector(XEdge2, YEdge2, 0);

    new G4PVPlacement(Rotation(180., 0., 180.), positionCristal1sb, logicCristal1s, Name + "_CsI_Cristal13", logicCsI, false, 13);
    new G4PVPlacement(Rotation(180, 0, 0), positionCristal2sb, logicCristal2s, Name + "_CsI_Cristal14", logicCsI, false, 14);
    new G4PVPlacement(Rotation(180., 0., 180.), positionCristal3sb, logicCristal3s, Name + "_CsI_Cristal15", logicCsI, false, 15);
    new G4PVPlacement(Rotation(180., 0., 180.), positionCristal4sb, logicCristal4s, Name + "_CsI_Cristal16", logicCsI, false, 16);

    ///Set CsI sensible
    logicCristal1->SetSensitiveDetector(m_CsIScorer);
    logicCristal2->SetSensitiveDetector(m_CsIScorer);
    logicCristal3->SetSensitiveDetector(m_CsIScorer);
    logicCristal4->SetSensitiveDetector(m_CsIScorer);

    logicCristal1s->SetSensitiveDetector(m_CsIScorer);
    logicCristal2s->SetSensitiveDetector(m_CsIScorer);
    logicCristal3s->SetSensitiveDetector(m_CsIScorer);
    logicCristal4s->SetSensitiveDetector(m_CsIScorer);

    //Mark blue a CsI corner crystal to see telescope orientation
    logicCristal1  ->SetVisAttributes(m_CsIVisAtt);
    logicCristal2  ->SetVisAttributes(m_CsIVisAtt);
    logicCristal3  ->SetVisAttributes(m_CsIVisAtt);
    logicCristal4  ->SetVisAttributes(m_CsIVisAtt);
    logicCristal1s ->SetVisAttributes(m_CsIVisAtt);
    logicCristal2s ->SetVisAttributes(m_CsIVisAtt);
    logicCristal3s ->SetVisAttributes(m_CsIVisAtt);
    logicCristal4s ->SetVisAttributes(m_CsIVisAtt);

  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Virtual Method of VDetector class

// Read stream at Configfile to pick-up parameters of detector (Position,...)
// Called in DetecorConstruction::ReadDetextorConfiguration Method
void MUST2Array::ReadConfiguration(string Path){
  ifstream ConfigFile           ;
  ConfigFile.open(Path.c_str()) ;
  string LineBuffer          ;
  string DataBuffer          ;

  // A:X1_Y1     --> X:1    Y:1
  // B:X128_Y1   --> X:128  Y:1
  // C:X1_Y128   --> X:1    Y:128
  // D:X128_Y128 --> X:128  Y:128

  G4double Ax , Bx , Cx , Dx , Ay , By , Cy , Dy , Az , Bz , Cz , Dz           	;
  G4ThreeVector A , B , C , D                                          			;
  G4double Theta = 0 , Phi = 0 , R = 0 , beta_u = 0 , beta_v = 0 , beta_w = 0      ;
  int SI = 0 , SILI = 0 , CSI = 0                                         			;	

  bool check_A = false ;
  bool check_C = false  ;
  bool check_B = false ;
  bool check_D = false  ;

  bool check_Theta = false   ;
  bool check_Phi  = false  ;
  bool check_R     = false   ;
  //   bool check_beta = false  ;

  bool check_SI = false   ;
  bool check_SILI  = false  ;
  bool check_CSI     = false   ;
  bool check_VIS = false  ;

  bool ReadingStatus = false ;


  while (!ConfigFile.eof()) {

    getline(ConfigFile, LineBuffer);

    //	If line is a Start Up MUST2 bloc, Reading toggle to true      
    if (LineBuffer.compare(0, 11, "M2Telescope") == 0) 
    {
      G4cout << "///" << G4endl           ;
      G4cout << "Telescope found: " << G4endl   ;        
      ReadingStatus = true ;

    }

    //	Else don't toggle to Reading Block Status
    else ReadingStatus = false ;

    //	Reading Block
    while(ReadingStatus)
    {

      ConfigFile >> DataBuffer ;
      //	Comment Line 
      if(DataBuffer.compare(0, 1, "%") == 0) {
        ConfigFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );

      }

      //	Finding another telescope (safety), toggle out
      else if (DataBuffer.compare(0, 11, "M2Telescope") == 0) {
        cout << "WARNING: Another Telescope is find before standard sequence of Token, Error may occured in Telecope definition" << endl ;
        ReadingStatus = false ;
      }

      //	Position method

      else if (DataBuffer.compare(0, 6, "X1_Y1=") == 0) {
        check_A = true;
        ConfigFile >> DataBuffer ;
        Ax = atof(DataBuffer.c_str()) ;
        Ax = Ax * mm ;
        ConfigFile >> DataBuffer ;
        Ay = atof(DataBuffer.c_str()) ;
        Ay = Ay * mm ;
        ConfigFile >> DataBuffer ;
        Az = atof(DataBuffer.c_str()) ;
        Az = Az * mm ;

        A = G4ThreeVector(Ax, Ay, Az);
        cout << "X1 Y1 corner position : " << A << endl;

      }


      else if (DataBuffer.compare(0, 8, "X128_Y1=") == 0) {
        check_B = true;
        ConfigFile >> DataBuffer ;
        Bx = atof(DataBuffer.c_str()) ;
        Bx = Bx * mm ;
        ConfigFile >> DataBuffer ;
        By = atof(DataBuffer.c_str()) ;
        By = By * mm ;
        ConfigFile >> DataBuffer ;
        Bz = atof(DataBuffer.c_str()) ;
        Bz = Bz * mm ;

        B = G4ThreeVector(Bx, By, Bz);
        cout << "X128 Y1 corner position : " << B << endl;

      }


      else if (DataBuffer.compare(0, 8, "X1_Y128=") == 0) {
        check_C = true;
        ConfigFile >> DataBuffer ;
        Cx = atof(DataBuffer.c_str()) ;
        Cx = Cx * mm ;
        ConfigFile >> DataBuffer ;
        Cy = atof(DataBuffer.c_str()) ;
        Cy = Cy * mm ;
        ConfigFile >> DataBuffer ;
        Cz = atof(DataBuffer.c_str()) ;
        Cz = Cz * mm ;

        C = G4ThreeVector(Cx, Cy, Cz);
        cout << "X1 Y128 corner position : " << C << endl;

      }

      else if (DataBuffer.compare(0, 10, "X128_Y128=") == 0) {
        check_D = true;
        ConfigFile >> DataBuffer ;
        Dx = atof(DataBuffer.c_str()) ;
        Dx = Dx * mm ;
        ConfigFile >> DataBuffer ;
        Dy = atof(DataBuffer.c_str()) ;
        Dy = Dy * mm ;
        ConfigFile >> DataBuffer ;
        Dz = atof(DataBuffer.c_str()) ;
        Dz = Dz * mm ;

        D = G4ThreeVector(Dx, Dy, Dz);
        cout << "X128 Y128 corner position : " << D << endl;

      }

      //	End Position Method

      //	Angle method
      else if (DataBuffer.compare(0, 6, "THETA=") == 0) {
        check_Theta = true;
        ConfigFile >> DataBuffer ;
        Theta = atof(DataBuffer.c_str()) ;
        Theta = Theta * deg;
        cout << "Theta:  " << Theta / deg << endl;

      }

      //Angle method
      else if (DataBuffer.compare(0, 4, "PHI=") == 0) {
        check_Phi = true;
        ConfigFile >> DataBuffer ;
        Phi = atof(DataBuffer.c_str()) ;
        Phi = Phi * deg;
        cout << "Phi:  " << Phi / deg << endl;

      }

      //Angle method
      else if (DataBuffer.compare(0, 2, "R=") == 0) {
        check_R = true;
        ConfigFile >> DataBuffer ;
        R = atof(DataBuffer.c_str()) ;
        R = R * mm;
        cout << "R:  " << R / mm << endl;

      }

      //Angle method
      else if (DataBuffer.compare(0, 5, "BETA=") == 0) {
        //		            check_beta = true;
        ConfigFile >> DataBuffer ;
        beta_u = atof(DataBuffer.c_str()) ;
        beta_u = beta_u * deg   ;
        ConfigFile >> DataBuffer ;
        beta_v = atof(DataBuffer.c_str()) ;
        beta_v = beta_v * deg   ;
        ConfigFile >> DataBuffer ;
        beta_w = atof(DataBuffer.c_str()) ;
        beta_w = beta_w * deg   ;
        G4cout << "Beta:  " << beta_u / deg << " " << beta_v / deg << " " << beta_w / deg << G4endl  ;

      }

      //	End Angle Method


      //	Common part
      else if (DataBuffer.compare(0, 3, "SI=") == 0) {
        check_SI = true ;
        ConfigFile >> DataBuffer;
        SI = atof(DataBuffer.c_str()) ;
        G4cout << " SI= " << SI << G4endl ;

      }


      else if (DataBuffer.compare(0, 5, "SILI=") == 0) {
        check_SILI = true ;
        ConfigFile >> DataBuffer;
        SILI = atof(DataBuffer.c_str()) ;
        G4cout << " SILI= " << SILI << G4endl ;

      }


      else if (DataBuffer.compare(0, 4, "CSI=") == 0) {
        check_CSI = true ;
        ConfigFile >> DataBuffer;
        CSI = atof(DataBuffer.c_str()) ;
        G4cout << " CSI= " << CSI << G4endl ;

      }


      else if (DataBuffer.compare(0, 4, "VIS=") == 0) {
        check_VIS = true ;
        ConfigFile >> DataBuffer;
        if (DataBuffer.compare(0, 3, "all") == 0) 
        {
          m_non_sensitive_part_visiualisation = true;
          G4cout << " VIS= all" << G4endl ;
        }

        else 
          G4cout << " VIS= Sensible Only" << G4endl ;	             

      }

      // If no MUST2 Token and no comment, toggle out
      else 
      {ReadingStatus = false; G4cout << "other token " << DataBuffer << G4endl ;}



      /////////////////////////////////////////////////
      //	If All necessary information there, toggle out
      if ( ( check_SI && check_SILI && check_CSI && check_VIS ) && ( (check_A && check_B && check_C && check_D) || (check_Theta && check_Phi && check_R) ) ) 
      { 
        ReadingStatus = false; 

        ///Add The previously define telescope
        //With position method
        if ((check_A && check_B && check_C && check_D) || !(check_Theta && check_Phi && check_R)) 
        {
          AddTelescope(  A,
              B,
              C,
              D,
              SI == 1,
              SILI == 1,
              CSI == 1);
        }

        //with angle method
        else if ((check_Theta && check_Phi && check_R) || !(check_A && check_B && check_C && check_D)) 
        {
          AddTelescope(  R,
              Theta,
              Phi,
              beta_u,
              beta_v,
              beta_w,
              SI == 1,
              SILI == 1,
              CSI == 1);
        }


        check_A = false 	;
        check_C = false  	;
        check_B = false 	;
        check_D = false  	;

        check_Theta = false   	;
        check_Phi  = false  	;
        check_R    = false   	;
        //	      				check_beta = false  	;

        check_SI = false   		;
        check_SILI  = false  	;
        check_CSI   = false   	;
        check_VIS = false  		;	

      }
    }
  }
}


// Construct detector and inialise sensitive part.
// Called After DetecorConstruction::AddDetector Method
void MUST2Array::ConstructDetector(G4LogicalVolume* world){
  G4RotationMatrix* MMrot  		= NULL ;
  G4ThreeVector     MMpos  		= G4ThreeVector(0, 0, 0) ;
  G4ThreeVector     MMu       	= G4ThreeVector(0, 0, 0) ;
  G4ThreeVector     MMv       	= G4ThreeVector(0, 0, 0) ;
  G4ThreeVector     MMw       	= G4ThreeVector(0, 0, 0) ;
  G4ThreeVector     MMCenter 	= G4ThreeVector(0, 0, 0) ;
  bool           	Si     		= true ;
  bool          	 	SiLi   		= true ;
  bool           	CsI    		= true ;

  G4int NumberOfTelescope = m_DefinitionType.size() ;

  for (G4int i = 0 ; i < NumberOfTelescope ; i++) {
    // By Point
    if (m_DefinitionType[i]) {
      // (u,v,w) unitary vector associated to telescope referencial
      // (u,v) // to silicon plan
      // w perpendicular to (u,v) plan and pointing CsI
      MMu = m_X128_Y1[i] - m_X1_Y1[i]     ;
      MMu = MMu.unit()              ;

      MMv = m_X1_Y128[i] - m_X1_Y1[i]     ;
      MMv = MMv.unit()              ;

      MMw = MMv.cross(MMu)          ;
      // if (MMw.z() > 0)MMw = MMv.cross(MMu)  ;
      MMw = MMw.unit()                ;

      MMCenter = (m_X1_Y1[i] + m_X1_Y128[i] + m_X128_Y1[i] + m_X128_Y128[i]) / 4 ;

      // Passage Matrix from Lab Referential to Telescope Referential
      MMrot = new G4RotationMatrix(MMv, MMu, MMw);
      MMpos = MMw * Length * 0.5 + MMCenter ;
    }

    // By Angle
    else {
      G4double Theta = m_Theta[i]   ;
      G4double Phi   = m_Phi[i]  ;

      // (u,v,w) unitary vector associated to telescope referencial
      // (u,v) // to silicon plan
      // w perpendicular to (u,v) plan and pointing ThirdStage
      // Phi is angle between X axis and projection in (X,Y) plan
      // Theta is angle between  position vector and z axis
      G4double wX = m_R[i] * sin(Theta / rad) * cos(Phi / rad) ;
      G4double wY = m_R[i] * sin(Theta / rad) * sin(Phi / rad) ;
      G4double wZ = m_R[i] * cos(Theta / rad);
      MMw = G4ThreeVector(wX, wY, wZ) ;

      // vector corresponding to the center of the module
      G4ThreeVector CT = MMw;

      // vector parallel to one axis of silicon plane
      G4double ii = cos(Theta / rad) * cos(Phi / rad);
      G4double jj = cos(Theta / rad) * sin(Phi / rad);
      G4double kk = -sin(Theta / rad);
      G4ThreeVector Y = G4ThreeVector(ii, jj, kk);

      MMw = MMw.unit();
      MMu = MMw.cross(Y);
      MMv = MMw.cross(MMu);
      MMv = MMv.unit();
      MMu = MMu.unit();

      // Passage Matrix from Lab Referential to Telescope Referential
      // MUST2
      MMrot = new G4RotationMatrix(MMu, MMv, MMw);
      // Telescope is rotate of Beta angle around MMv axis.
      MMrot->rotate(m_beta_u[i], MMu);
      MMrot->rotate(m_beta_v[i], MMv);
      MMrot->rotate(m_beta_w[i], MMw);
      // translation to place Telescope
      MMpos = MMw * Length * 0.5 + CT ;
    }

    Si  = m_wSi[i]    ;
    SiLi = m_wSiLi[i] ;
    CsI  = m_wCsI[i]  ;

    VolumeMaker(i + 1 , MMpos , MMrot , Si , SiLi , CsI , world);
  }

  delete MMrot ;
}

// Add Detector branch to the EventTree.
// Called After DetecorConstruction::AddDetector Method

void MUST2Array::InitializeRootOutput(){
  RootOutput *pAnalysis = RootOutput::getInstance();
  TTree *pTree = pAnalysis->GetTree();
  pTree->Branch("MUST2", "TMust2Data", &m_Event) ;
}

// Read sensitive part and fill the Root tree.
// Called at in the EventAction::EndOfEventAvtion
void MUST2Array::ReadSensitive(const G4Event* event){
  G4String DetectorNumber 	;
  m_Event->Clear()			;

  //////////////////////////////////////////////////////////////////////////////////////
  //////////////////////// Used to Read Event Map of detector //////////////////////////
  //////////////////////////////////////////////////////////////////////////////////////

  // Si
  std::map<G4int, G4int*>::iterator DetectorNumber_itr;
  std::map<G4int, G4double*>::iterator Energy_itr;
  std::map<G4int, G4double*>::iterator Time_itr;
  std::map<G4int, G4int*>::iterator X_itr;
  std::map<G4int, G4int*>::iterator Y_itr;
  std::map<G4int, G4double*>::iterator Pos_X_itr;
  std::map<G4int, G4double*>::iterator Pos_Y_itr;
  std::map<G4int, G4double*>::iterator Pos_Z_itr;
  std::map<G4int, G4double*>::iterator Ang_Theta_itr;
  std::map<G4int, G4double*>::iterator Ang_Phi_itr;

  G4THitsMap<G4int>*	  DetectorNumberHitMap;
  G4THitsMap<G4double>* EnergyHitMap;
  G4THitsMap<G4double>* TimeHitMap;
  G4THitsMap<G4int>* XHitMap;
  G4THitsMap<G4int>* YHitMap;
  G4THitsMap<G4double>* PosXHitMap;
  G4THitsMap<G4double>* PosYHitMap;
  G4THitsMap<G4double>* PosZHitMap;
  G4THitsMap<G4double>* AngThetaHitMap;
  G4THitsMap<G4double>* AngPhiHitMap;

  // Si(Li)
  std::map<G4int, G4double*>::iterator SiLiEnergy_itr;
  std::map<G4int, G4int*>::iterator SiLiPadNbr_itr;
  G4THitsMap<G4double>* 				 SiLiEnergyHitMap;
  G4THitsMap<G4int>* 				    SiLiPadNbrHitMap;

  // CsI
  std::map<G4int, G4double*>::iterator CsIEnergy_itr;
  std::map<G4int, G4int*>::iterator    CsICristalNbr_itr;
  G4THitsMap<G4double>* 				    CsIEnergyHitMap;
  G4THitsMap<G4int>* 				       CsICristalNbrHitMap;
  //////////////////////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////////////////////

  // Read the Scorer associate to the Silicon Strip

  //Detector Number
  G4int StripDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/DetectorNumber");
  DetectorNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripDetCollectionID));
  DetectorNumber_itr =  DetectorNumberHitMap->GetMap()->begin();

  //Energy
  G4int StripEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripEnergy") ;
  EnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripEnergyCollectionID)) ;
  Energy_itr = EnergyHitMap->GetMap()->begin();

  //Time of Flight
  G4int StripTimeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripTime");
  TimeHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripTimeCollectionID));
  Time_itr = TimeHitMap->GetMap()->begin();

  //Strip Number X
  G4int StripXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripNumberX") ;
  XHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripXCollectionID));
  X_itr = XHitMap->GetMap()->begin();

  //Strip Number Y
  G4int StripYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripNumberY");
  YHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripYCollectionID));
  Y_itr = YHitMap->GetMap()->begin();

  //Interaction Coordinate X
  G4int InterCoordXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordX");
  PosXHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordXCollectionID));
  Pos_X_itr = PosXHitMap->GetMap()->begin();

  //Interaction Coordinate Y
  G4int InterCoordYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordY");
  PosYHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordYCollectionID));
  Pos_Y_itr = PosYHitMap->GetMap()->begin();

  //Interaction Coordinate Z
  G4int InterCoordZCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordZ");
  PosZHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordZCollectionID));
  Pos_Z_itr = PosXHitMap->GetMap()->begin();

  //Interaction Coordinate Angle Theta
  G4int InterCoordAngThetaCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordAngTheta");
  AngThetaHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngThetaCollectionID));
  Ang_Theta_itr = AngThetaHitMap->GetMap()->begin();

  //Interaction Coordinate Angle Phi
  G4int InterCoordAngPhiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordAngPhi");	
  AngPhiHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngPhiCollectionID));
  Ang_Phi_itr = AngPhiHitMap->GetMap()->begin();			

  // Read the Scorer associate to the SiLi
  //Energy
  G4int SiLiEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_SiLiScorer/SiLiEnergy");
  SiLiEnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(SiLiEnergyCollectionID));
  SiLiEnergy_itr = SiLiEnergyHitMap->GetMap()->begin();

  G4int SiLiPadNbrCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_SiLiScorer/SiLiPadNbr");
  SiLiPadNbrHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(SiLiPadNbrCollectionID));
  SiLiPadNbr_itr = SiLiPadNbrHitMap->GetMap()->begin();

  // Read the Scorer associate to the CsI crystal
  //Energy
  G4int CsIEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_CsIScorer/CsIEnergy");
  CsIEnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(CsIEnergyCollectionID));
  CsIEnergy_itr = CsIEnergyHitMap->GetMap()->begin();

  G4int CsICristalNbrCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_CsIScorer/CsICristalNbr");
  CsICristalNbrHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(CsICristalNbrCollectionID));
  CsICristalNbr_itr = CsICristalNbrHitMap->GetMap()->begin();


  /////////////////////

  G4int sizeN = DetectorNumberHitMap->entries();
  G4int sizeE = EnergyHitMap->entries();
  G4int sizeT = TimeHitMap->entries();
  G4int sizeX = XHitMap->entries();
  G4int sizeY = YHitMap->entries();

  // Loop on Telescope Number entry
  for (G4int l = 0 ; l < sizeN ; l++) {
    G4int N         =  *(DetectorNumber_itr->second);
    G4int NTrackID  =   DetectorNumber_itr->first - N;


    if (N > 0) {

      m_Event->SetMMStripXEDetectorNbr(N) ;
      m_Event->SetMMStripYEDetectorNbr(N) ;
      m_Event->SetMMStripXTDetectorNbr(N) ;
      m_Event->SetMMStripYTDetectorNbr(N) ;

      //  Energy
      Energy_itr = EnergyHitMap->GetMap()->begin();
      for (G4int h = 0 ; h < sizeE ; h++) {
        G4int ETrackID  =   Energy_itr->first - N;
        G4double E      = *(Energy_itr->second);

        if (ETrackID == NTrackID) {
          m_Event->SetMMStripXEEnergy(RandGauss::shoot(E, ResoStrip));
          m_Event->SetMMStripYEEnergy(RandGauss::shoot(E, ResoStrip));
        }

        Energy_itr++;
      }


      //  Time
      Time_itr = TimeHitMap->GetMap()->begin();
      for (G4int h = 0 ; h < sizeT ; h++) {
        G4int TTrackID  =   Time_itr->first - N;
        G4double T      = *(Time_itr->second);

        if (TTrackID == NTrackID) {
          m_Event->SetMMStripXTTime(RandGauss::shoot(T, ResoTimeMust)) ;
          m_Event->SetMMStripYTTime(RandGauss::shoot(T, ResoTimeMust)) ;
        }

        Time_itr++;
      }


      // X
      X_itr = XHitMap->GetMap()->begin();
      for (G4int h = 0 ; h < sizeX ; h++) {
        G4int XTrackID  =   X_itr->first  - N;
        G4int X         = *(X_itr->second);
        if (XTrackID == NTrackID) {
          m_Event->SetMMStripXEStripNbr(X);
          m_Event->SetMMStripXTStripNbr(X);
        }

        X_itr++;
      }

      // Y
      Y_itr = YHitMap->GetMap()->begin()  ;
      for (G4int h = 0 ; h < sizeY ; h++) {
        G4int YTrackID  =   Y_itr->first  - N ;
        G4int Y         = *(Y_itr->second);
        if (YTrackID == NTrackID) {
          m_Event->SetMMStripYEStripNbr(Y);
          m_Event->SetMMStripYTStripNbr(Y);
        }

        Y_itr++;
      }

      // Pos X
      Pos_X_itr = PosXHitMap->GetMap()->begin();
      for (G4int h = 0 ; h < PosXHitMap->entries() ; h++) {
        G4int PosXTrackID =   Pos_X_itr->first  - N ;
        G4double PosX     = *(Pos_X_itr->second) ;
        if (PosXTrackID == NTrackID) {
          ms_InterCoord->SetDetectedPositionX(PosX) ;
        }
        Pos_X_itr++;
      }

      // Pos Y
      Pos_Y_itr = PosYHitMap->GetMap()->begin();
      for (G4int h = 0 ; h < PosYHitMap->entries() ; h++) {
        G4int PosYTrackID =   Pos_Y_itr->first  - N ;
        G4double PosY     = *(Pos_Y_itr->second) ;
        if (PosYTrackID == NTrackID) {
          ms_InterCoord->SetDetectedPositionY(PosY) ;
        }
        Pos_Y_itr++;
      }

      // Pos Z
      Pos_Z_itr = PosZHitMap->GetMap()->begin();
      for (G4int h = 0 ; h < PosZHitMap->entries() ; h++) {
        G4int PosZTrackID =   Pos_Z_itr->first   - N  ;
        G4double PosZ     = *(Pos_Z_itr->second)      ;
        if (PosZTrackID == NTrackID) {
          ms_InterCoord->SetDetectedPositionZ(PosZ) ;
        }
        Pos_Z_itr++;
      }

      // Angle Theta
      Ang_Theta_itr = AngThetaHitMap->GetMap()->begin();
      for (G4int h = 0 ; h < AngThetaHitMap->entries(); h++) {
        G4int AngThetaTrackID =   Ang_Theta_itr->first  - N   ;
        G4double AngTheta     = *(Ang_Theta_itr->second)      ;
        if (AngThetaTrackID == NTrackID) {
          ms_InterCoord->SetDetectedAngleTheta(AngTheta) ;
        }
        Ang_Theta_itr++;
      }

      // Angle Phi
      Ang_Phi_itr = AngPhiHitMap->GetMap()->begin();
      for (G4int h = 0 ; h < AngPhiHitMap->entries() ; h++) {
        G4int AngPhiTrackID =   Ang_Phi_itr->first  - N   ;
        G4double AngPhi     = *(Ang_Phi_itr->second)      ;
        if (AngPhiTrackID == NTrackID) {
          ms_InterCoord->SetDetectedAnglePhi(AngPhi) ;
        }
        Ang_Phi_itr++;
      }

      // Si(Li)
      SiLiEnergy_itr = SiLiEnergyHitMap->GetMap()->begin() ;
      SiLiPadNbr_itr = SiLiPadNbrHitMap->GetMap()->begin() ;

      for (G4int h = 0 ; h < SiLiEnergyHitMap->entries() ; h++) {
        G4int SiLiEnergyTrackID =   SiLiEnergy_itr->first  -N ;
        G4double SiLiEnergy     = *(SiLiEnergy_itr->second)   ;

        G4int PadNbr = *(SiLiPadNbr_itr->second)   ;

        if (SiLiEnergyTrackID == NTrackID) {
          m_Event->SetMMSiLiEEnergy(RandGauss::shoot(SiLiEnergy, ResoSiLi)) ;
          m_Event->SetMMSiLiEPadNbr(PadNbr);
          m_Event->SetMMSiLiTPadNbr(PadNbr);
          m_Event->SetMMSiLiTTime(1);
          m_Event->SetMMSiLiTDetectorNbr(N);
          m_Event->SetMMSiLiEDetectorNbr(N);
        }

        SiLiEnergy_itr++	;
        SiLiPadNbr_itr++	;
      }

      // CsI    
      CsIEnergy_itr = CsIEnergyHitMap->GetMap()->begin()  				;
      CsICristalNbr_itr = CsICristalNbrHitMap->GetMap()->begin()  ;

      for (G4int h = 0 ; h < CsIEnergyHitMap->entries() ; h++) {
        G4int CsIEnergyTrackID  =   CsIEnergy_itr->first-N ;
        G4double CsIEnergy      = *(CsIEnergy_itr->second) ;

        G4int CristalNumber  = *(CsICristalNbr_itr->second) ;
        if (CsIEnergyTrackID == NTrackID) {
          m_Event->SetMMCsIEEnergy(RandGauss::shoot(CsIEnergy, ResoCsI*sqrt(CsIEnergy)));
          m_Event->SetMMCsIECristalNbr(CristalNumber);
          m_Event->SetMMCsITCristalNbr(CristalNumber);
          m_Event->SetMMCsITTime(1);
          m_Event->SetMMCsITDetectorNbr(N);
          m_Event->SetMMCsIEDetectorNbr(N);
        }

        CsIEnergy_itr++			;
        CsICristalNbr_itr++	;
      }

    }

    DetectorNumber_itr++;
  }

  // clear map for next event
  DetectorNumberHitMap	->clear();
  EnergyHitMap   			->clear() ;
  TimeHitMap     			->clear() ;
  XHitMap        			->clear() ;
  YHitMap        			->clear() ;
  SiLiEnergyHitMap			->clear() ;
  SiLiPadNbrHitMap			->clear() ;
  CsIEnergyHitMap			->clear() ;
  CsICristalNbrHitMap		->clear() ;
  PosXHitMap					->clear() ; 
  PosYHitMap					->clear() ;
  PosZHitMap					->clear() ;
  AngThetaHitMap			->clear() ;
  AngPhiHitMap				->clear() ;

}


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MUST2Array::InitializeScorers() { 
  //	Silicon Associate Scorer
  m_StripScorer = new G4MultiFunctionalDetector("MUST2_StripScorer");

  G4VPrimitiveScorer* DetNbr 									= new OBSOLETEGENERALSCORERS::PSDetectorNumber("DetectorNumber","MUST2Telescope", 0);
  G4VPrimitiveScorer* Energy 									= new OBSOLETEGENERALSCORERS::PSEnergy("StripEnergy","MUST2Telescope", 0);			
  G4VPrimitiveScorer* TOF 										= new OBSOLETEGENERALSCORERS::PSTOF("StripTime","MUST2Telescope", 0);          					 		 

  G4VPrimitiveScorer* StripPositionX							= new PSStripNumberX("StripNumberX", 0, ActiveWaferSize, 128);
  G4VPrimitiveScorer* StripPositionY							= new PSStripNumberY("StripNumberY", 0, ActiveWaferSize, 128);  		

  G4VPrimitiveScorer* InteractionCoordinatesX 				= new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesX("InterCoordX","MUST2Telescope", 0);
  G4VPrimitiveScorer* InteractionCoordinatesY				= new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesY("InterCoordY","MUST2Telescope", 0);
  G4VPrimitiveScorer* InteractionCoordinatesZ  			= new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesZ("InterCoordZ","MUST2Telescope", 0);	 		 
  G4VPrimitiveScorer* InteractionCoordinatesAngleTheta	= new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesAngleTheta("InterCoordAngTheta","MUST2Telescope", 0);
  G4VPrimitiveScorer* InteractionCoordinatesAnglePhi    = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesAnglePhi("InterCoordAngPhi","MUST2Telescope", 0) ;	    


  //and register it to the multifunctionnal detector
  m_StripScorer->RegisterPrimitive(DetNbr);
  m_StripScorer->RegisterPrimitive(Energy);
  m_StripScorer->RegisterPrimitive(TOF);
  m_StripScorer->RegisterPrimitive(StripPositionX);
  m_StripScorer->RegisterPrimitive(StripPositionY);
  m_StripScorer->RegisterPrimitive(InteractionCoordinatesX);
  m_StripScorer->RegisterPrimitive(InteractionCoordinatesY);
  m_StripScorer->RegisterPrimitive(InteractionCoordinatesZ);
  m_StripScorer->RegisterPrimitive(InteractionCoordinatesAngleTheta);
  m_StripScorer->RegisterPrimitive(InteractionCoordinatesAnglePhi);

  //	SiLi Associate Scorer
  m_SiLiScorer	= new G4MultiFunctionalDetector("MUST2_SiLiScorer");
  G4VPrimitiveScorer* SiLiEnergy 			= new OBSOLETEGENERALSCORERS::PSEnergy("SiLiEnergy","MUST2Telescope", 0) ;
  G4VPrimitiveScorer* SiLiPadNbr 			= new PSPadOrCristalNumber("SiLiPadNbr",0) ;
  m_SiLiScorer->RegisterPrimitive(SiLiEnergy);
  m_SiLiScorer->RegisterPrimitive(SiLiPadNbr);

  //	CsI Associate Scorer 
  m_CsIScorer	= new G4MultiFunctionalDetector("MUST2_CsIScorer");
  G4VPrimitiveScorer* CsIEnergy 		= new OBSOLETEGENERALSCORERS::PSEnergy("CsIEnergy","MUST2Telescope", 0) 	;
  G4VPrimitiveScorer* CsICristalNbr 	= new PSPadOrCristalNumber("CsICristalNbr",0) ;
  m_CsIScorer->RegisterPrimitive(CsIEnergy) ;
  m_CsIScorer->RegisterPrimitive(CsICristalNbr) ;

  //	Add All Scorer to the Global Scorer Manager
  G4SDManager::GetSDMpointer()->AddNewDetector(m_StripScorer) ;
  G4SDManager::GetSDMpointer()->AddNewDetector(m_SiLiScorer)	 ;
  G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIScorer)	 ;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MUST2Array::InitializeMaterial(){

  ////////////////////////////////////////////////////////////////
  /////////////////Element  Definition ///////////////////////////
  ////////////////////////////////////////////////////////////////
  G4String symbol               ;
  G4double density = 0. , a = 0, z = 0   ;
  G4int ncomponents = 0, natoms = 0  ;

  G4Element* H   = new G4Element("Hydrogen" , symbol = "H"  , z = 1  , a = 1.01   * g / mole);
  G4Element* C   = new G4Element("Carbon"   , symbol = "C"  , z = 6  , a = 12.011 * g / mole);
  G4Element* N   = new G4Element("Nitrogen" , symbol = "N"  , z = 7  , a = 14.01  * g / mole);
  G4Element* O   = new G4Element("Oxigen"   , symbol = "O"  , z = 8  , a = 16.00  * g / mole);
  G4Element* I   = new G4Element("Iode"     , symbol = "I"  , z = 53 , a = 126.9  * g / mole);
  G4Element* Cs  = new G4Element("Cesium"   , symbol = "Cs" , z = 55 , a = 132.9  * g / mole);
  G4Element* Fe  = new G4Element("Iron"     , symbol = "Fe" , z = 26 , a = 55.847 * g / mole);

  ////////////////////////////////////////////////////////////////
  /////////////////Material Definition ///////////////////////////
  ////////////////////////////////////////////////////////////////

  // Si
  a = 28.0855 * g / mole;
  density = 2.321 * g / cm3;
  m_MaterialSilicon = new G4Material("Si", z = 14., a, density);

  // Al
  density = 2.702 * g / cm3;
  a = 26.98 * g / mole;
  m_MaterialAluminium = new G4Material("Aluminium", z = 13., a, density);

  // Iron
  density = 7.874 * g / cm3;
  a = 55.847 * g / mole;
  m_MaterialIron = new G4Material("Iron", z = 26., a, density);

  // CsI
  density = 4.51 * g / cm3;
  m_MaterialCsI = new G4Material("CsI", density, ncomponents = 2);
  m_MaterialCsI->AddElement(Cs , natoms = 1);
  m_MaterialCsI->AddElement(I  , natoms = 1);

  //  Vacuum
  density = 0.000000001 * mg / cm3;
  m_MaterialVacuum = new G4Material("Vacuum", density, ncomponents = 2);
  m_MaterialVacuum->AddElement(N, .7);
  m_MaterialVacuum->AddElement(O, .3);

  //  Mylar
  density = 1.397 * g / cm3;
  m_MaterialMyl = new G4Material("Mylar", density, ncomponents = 3);
  m_MaterialMyl->AddElement(C, natoms = 10);
  m_MaterialMyl->AddElement(H, natoms = 8);
  m_MaterialMyl->AddElement(O, natoms = 4);

}


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4RotationMatrix* Rotation(double tetaX, double tetaY, double tetaZ){
  double PI = 3.141592653589793238;
  double radX = tetaX * PI / 180.;
  double radY = tetaY * PI / 180.;
  double radZ = tetaZ * PI / 180.;

  G4ThreeVector col1 = G4ThreeVector(cos(radZ) * cos(radY),
      -sin(radZ) * cos(radY),
      -sin(radY));
  G4ThreeVector col2 = G4ThreeVector(sin(radZ) * cos(radX) - cos(radZ) * sin(radY) * sin(radX),
      cos(radZ) * cos(radX) + sin(radZ) * sin(radY) * sin(radX),
      -cos(radY) * sin(radX));
  G4ThreeVector col3 = G4ThreeVector(sin(radZ) * sin(radX) + cos(radZ) * sin(radY) * sin(radX),
      cos(radZ) * sin(radX) - sin(radZ) * sin(radY) * cos(radX),
      cos(radY) * cos(radX));

  return (new G4RotationMatrix(col1, col2, col3));
}