Finished Clover Position

CMakeLists.txt

@@ -0,0 +1,65 @@
+#----------------------------------------------------------------------------
+# Setup the project
+#
+cmake_minimum_required(VERSION 3.8...3.18)
+if(${CMAKE_VERSION} VERSION_LESS 3.12)
+  cmake_policy(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION})
+endif()
+project(B4c)
+
+#----------------------------------------------------------------------------
+# Find Geant4 package, activating all available UI and Vis drivers by default
+# You can set WITH_GEANT4_UIVIS to OFF via the command line or ccmake/cmake-gui
+# to build a batch mode only executable
+#
+option(WITH_GEANT4_UIVIS "Build example with Geant4 UI and Vis drivers" ON)
+if(WITH_GEANT4_UIVIS)
+  find_package(Geant4 REQUIRED ui_all vis_all)
+else()
+  find_package(Geant4 REQUIRED)
+endif()
+
+#----------------------------------------------------------------------------
+# Setup Geant4 include directories and compile definitions
+# Setup include directory for this project
+#
+include(${Geant4_USE_FILE})
+include_directories(${PROJECT_SOURCE_DIR}/include)
+
+#----------------------------------------------------------------------------
+# Locate sources and headers for this project
+# NB: headers are included so they will show up in IDEs
+#
+file(GLOB sources ${PROJECT_SOURCE_DIR}/*.cc)
+file(GLOB headers ${PROJECT_SOURCE_DIR}/*.hh)
+
+#----------------------------------------------------------------------------
+# Add the executable, and link it to the Geant4 libraries
+#
+add_executable(clover clover.cc ${sources} ${headers})
+target_link_libraries(clover ${Geant4_LIBRARIES})
+
+#----------------------------------------------------------------------------
+# Copy all scripts to the build directory, i.e. the directory in which we
+# build Clover. This is so that we can run the executable directly because it
+# relies on these scripts being in the current working directory.
+
+set(CLOVER_SCRIPTS
+  init_vis.mac
+  vis.mac
+  my.in
+  analysis.C
+)
+
+foreach(_script ${CLOVER_SCRIPTS})
+  configure_file(
+    ${PROJECT_SOURCE_DIR}/${_script}
+    ${PROJECT_BINARY_DIR}/${_script}
+    COPYONLY
+    )
+endforeach()
+
+#----------------------------------------------------------------------------
+# Install the executable to 'bin' directory under CMAKE_INSTALL_PREFIX
+#
+install(TARGETS clover DESTINATION bin)

CloverActionInitialization.cc

@@ -0,0 +1,40 @@
+
+
+#include "CloverActionInitialization.hh"
+#include "CloverPrimaryGeneratorAction.hh"
+#include "CloverRunAction.hh"
+#include "CloverEventAction.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverActionInitialization::CloverActionInitialization()
+ : G4VUserActionInitialization()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverActionInitialization::~CloverActionInitialization()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverActionInitialization::BuildForMaster() const
+{
+  auto eventAction = new CloverEventAction;
+  SetUserAction(new CloverRunAction(eventAction));
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverActionInitialization::Build() const
+{
+
+
+  auto eventAction = new CloverEventAction;
+  SetUserAction(eventAction);
+
+  SetUserAction(new CloverPrimaryGeneratorAction());
+  SetUserAction(new CloverRunAction(eventAction));
+}  
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

CloverActionInitialization.hh

@@ -0,0 +1,22 @@
+
+#ifndef CloverActionInitialization_h
+#define CloverActionInitialization_h 1
+
+#include "G4VUserActionInitialization.hh"
+
+/// Action initialization class.
+///
+
+class CloverActionInitialization : public G4VUserActionInitialization
+{
+  public:
+    CloverActionInitialization();
+    virtual ~CloverActionInitialization();
+
+    virtual void BuildForMaster() const;
+    virtual void Build() const;
+};
+
+#endif
+
+    

CloverAnalysis.hh

@@ -0,0 +1,12 @@
+
+#ifndef CloverAnalysis_h
+#define CloverAnalysis_h 1
+
+#include "g4root.hh"
+//#include "g4csv.hh"
+//#include "g4xml.hh"
+
+//#include "G4GenericAnalysisManager.hh"
+//using G4AnalysisManager = G4GenericAnalysisManager;
+
+#endif

CloverCrystalHit.cc

@@ -0,0 +1,85 @@
+
+#include "CloverCrystalHit.hh"
+
+#include "G4UnitsTable.hh"
+#include "G4VVisManager.hh"
+#include "G4Circle.hh"
+#include "G4Colour.hh"
+#include "G4VisAttributes.hh"
+
+#include <iomanip>
+
+G4ThreadLocal G4Allocator<CloverCrystalHit>* CloverCrystalHitAllocator = 0;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverCrystalHit::CloverCrystalHit()
+ : G4VHit(),
+   fTrackID(-1),
+   fCrysalID(-1),
+   fEdep(0.),
+   fPos(G4ThreeVector()),
+   fStepLength(0.),
+   fBeamEnergy(0.)
+{
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverCrystalHit::~CloverCrystalHit() {}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverCrystalHit::CloverCrystalHit(const CloverCrystalHit& right)
+  : G4VHit()
+{
+  fTrackID   = right.fTrackID;
+  fCrysalID  = right.fCrysalID;
+  fEdep      = right.fEdep;
+  fPos       = right.fPos;
+
+  fStepLength = right.fStepLength;
+  fBeamEnergy = right.fBeamEnergy;
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+const CloverCrystalHit& CloverCrystalHit::operator=(const CloverCrystalHit& right)
+{
+  fTrackID   = right.fTrackID;
+  fCrysalID  = right.fCrysalID;
+  fEdep      = right.fEdep;
+  fPos       = right.fPos;
+
+  fStepLength = right.fStepLength;
+  fBeamEnergy = right.fBeamEnergy;
+
+  return *this;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4bool CloverCrystalHit::operator==(const CloverCrystalHit& right) const
+{
+  return ( this == &right ) ? true : false;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverCrystalHit::Print()
+{
+  G4cout
+     << "Beam Energy : " << fBeamEnergy
+     << "  trackID: " << fTrackID << " Crystal: " << fCrysalID 
+     << "  Edep: "
+     << std::setw(7) << G4BestUnit(fEdep,"Energy")
+     << "  StepLen: "
+     << std::setw(7) << G4BestUnit(fStepLength,"Length")
+     << "  Position: "
+     << std::setw(7) << G4BestUnit( fPos,"Length")
+     << G4endl;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

CloverCrystalHit.hh

@@ -0,0 +1,91 @@
+
+#ifndef CloverCrystalHit_h
+#define CloverCrystalHit_h 1
+
+#include "G4VHit.hh"
+#include "G4THitsCollection.hh"
+#include "G4Allocator.hh"
+#include "G4ThreeVector.hh"
+#include "G4Threading.hh"
+
+#include <vector>
+
+/// Calorimeter hit class
+///
+/// It defines data members to store the the energy deposit and track lengths
+/// of charged particles in a selected volume:
+/// - fEdep, fTrackLength
+
+class CloverCrystalHit : public G4VHit
+{
+  public:
+    CloverCrystalHit();
+    CloverCrystalHit(const CloverCrystalHit&);
+    virtual ~CloverCrystalHit();
+
+    // operators
+    const CloverCrystalHit& operator=(const CloverCrystalHit&);
+    G4bool operator==(const CloverCrystalHit&) const;
+
+    inline void* operator new(size_t);
+    inline void  operator delete(void*);
+
+    // methods from base class
+    virtual void Draw() {}
+    virtual void Print();
+
+    // Set methods
+    void SetTrackID    (G4int track)      { fTrackID = track; };
+    void SetCrystalID  (G4int id)         { fCrysalID = id; };
+    void SetEdep       (G4double de)      { fEdep += de; };
+    void SetPos        (G4ThreeVector xyz){ fPos = xyz; };
+    
+    void SetStepLength (G4double sl)      { fStepLength += sl; };
+    void SetBeamEnergy (G4double be)      { fBeamEnergy = be;};
+
+    // Get methods
+    G4int GetTrackID() const     { return fTrackID; };
+    G4int GetChamberNb() const   { return fCrysalID; };
+    G4double GetEdep() const     { return fEdep; };
+    G4ThreeVector GetPos() const { return fPos; };
+
+    G4double GetStepLength () const  { return fStepLength ; };
+    G4double GetBeamEnergy () const  { return fBeamEnergy ; };
+
+  private:
+
+    G4int         fTrackID;
+    G4int         fCrysalID;
+    G4double      fEdep;
+    G4ThreeVector fPos;
+    G4double      fStepLength;
+    G4double      fBeamEnergy;
+
+
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+using CloverCrystalHitsCollection = G4THitsCollection<CloverCrystalHit>;
+
+extern G4ThreadLocal G4Allocator<CloverCrystalHit>* CloverCrystalHitAllocator;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+inline void* CloverCrystalHit::operator new(size_t)
+{
+  if (!CloverCrystalHitAllocator) CloverCrystalHitAllocator = new G4Allocator<CloverCrystalHit>;
+  
+  return (void *) CloverCrystalHitAllocator->MallocSingle();
+}
+
+inline void CloverCrystalHit::operator delete(void *hit)
+{
+  if (!CloverCrystalHitAllocator) CloverCrystalHitAllocator = new G4Allocator<CloverCrystalHit>;
+  
+  CloverCrystalHitAllocator->FreeSingle((CloverCrystalHit*) hit);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif

CloverCrystalSD.cc

@@ -0,0 +1,106 @@
+
+#include "CloverCrystalSD.hh"
+
+#include "G4HCofThisEvent.hh"
+#include "G4Step.hh"
+#include "G4ThreeVector.hh"
+#include "G4SDManager.hh"
+#include "G4ios.hh"
+#include "G4SystemOfUnits.hh"
+
+///....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverCrystalSD::CloverCrystalSD(const G4String& name, const G4String& hitsCollectionName, const G4int nCrystal)
+ : G4VSensitiveDetector(name),
+   fHitsCollection(nullptr),
+   fNDet(nCrystal)
+{
+  collectionName.insert(hitsCollectionName);
+}
+
+///....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverCrystalSD::~CloverCrystalSD() 
+{ 
+}
+
+///....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverCrystalSD::Initialize(G4HCofThisEvent* hce)
+{
+  // Create hits collection
+  fHitsCollection = new CloverCrystalHitsCollection(SensitiveDetectorName, collectionName[0]); 
+
+  // Add this collection in hce
+  G4int hcID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);  
+  hce->AddHitsCollection( hcID, fHitsCollection );
+
+  for (G4int i=0; i<fNDet; i++ ) {
+    fHitsCollection->insert(new CloverCrystalHit());
+  }
+
+  //G4cout << "######### size of fHitCollection : " <<  fHitsCollection->GetSize() << G4endl;
+
+}
+
+///....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4bool CloverCrystalSD::ProcessHits(G4Step* step, G4TouchableHistory* /*history*/)
+{
+
+  // energy deposit
+  G4double edep = step->GetTotalEnergyDeposit();
+  
+  // step length
+  G4double stepLength = 0.;
+  if ( step->GetTrack()->GetDefinition()->GetPDGCharge() != 0. ) {
+    stepLength = step->GetStepLength();    
+  }
+
+  if ( edep==0. && stepLength == 0. ) return false;
+
+  G4int crystalID = step->GetPreStepPoint()->GetTouchableHandle() ->GetCopyNumber();
+  //----------- save hit in each crystal
+  CloverCrystalHit * hit = (*fHitsCollection)[crystalID];
+
+  hit->SetTrackID  (step->GetTrack()->GetTrackID());
+  hit->SetEdep(edep);
+  hit->SetPos (step->GetPostStepPoint()->GetPosition());
+  hit->SetStepLength( stepLength);
+  hit->SetCrystalID(crystalID);
+
+
+  //---------- Save indivual hit
+  CloverCrystalHit* newHit = new CloverCrystalHit();
+  
+  newHit->SetTrackID  (step->GetTrack()->GetTrackID());
+  newHit->SetCrystalID(crystalID);
+  newHit->SetEdep(edep);
+  newHit->SetPos (step->GetPostStepPoint()->GetPosition());
+  newHit->SetStepLength( stepLength);
+  
+  fHitsCollection->insert( newHit );
+
+
+  return true;
+}
+
+///....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverCrystalSD::EndOfEvent(G4HCofThisEvent*)
+{
+  if ( verboseLevel > 1 ) { 
+     G4int nofHits = fHitsCollection->GetSize();
+     G4cout
+       << G4endl 
+       << ">>>>>>>>>>Hits Collection: in this event they are " << nofHits 
+       << " hits in the tracker chambers: " << G4endl;
+     for ( G4int i=0; i<nofHits; ++i ) {
+       if( i == fNDet ) G4cout << "---------------------------" << G4endl;
+       G4cout << i << " |";
+       (*fHitsCollection)[i]->Print();
+     }
+  }
+}
+
+///....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

CloverCrystalSD.hh

@@ -0,0 +1,32 @@
+
+#ifndef CloverCrystalSD_h
+#define CloverCrystalSD_h 1
+
+#include "G4VSensitiveDetector.hh"
+
+#include "CloverCrystalHit.hh"
+
+#include <vector>
+
+class G4Step;
+class G4HCofThisEvent;
+
+class CloverCrystalSD : public G4VSensitiveDetector
+{
+  public:
+    CloverCrystalSD(const G4String& name, const G4String& hitsCollectionName, const G4int nCrystal);
+    virtual ~CloverCrystalSD();
+  
+    // methods from base class
+    virtual void   Initialize(G4HCofThisEvent* hitCollection);
+    virtual G4bool ProcessHits(G4Step* step, G4TouchableHistory* history);
+    virtual void   EndOfEvent(G4HCofThisEvent* hitCollection);
+
+  private:
+    CloverCrystalHitsCollection* fHitsCollection;
+
+    G4int fNDet;
+};
+
+#endif
+

CloverDetectorConstruction.cc

@@ -0,0 +1,336 @@
+
+#include "CloverDetectorConstruction.hh"
+#include "CloverCrystalSD.hh"
+#include "Detector_Parameters.hh"
+
+#include "G4Material.hh"
+#include "G4NistManager.hh"
+
+#include "G4Box.hh"
+#include "G4Tubs.hh"
+#include "G4Sphere.hh"
+#include "G4IntersectionSolid.hh"
+#include "G4SubtractionSolid.hh"
+
+#include "G4LogicalVolume.hh"
+#include "G4PVPlacement.hh"
+#include "G4PVReplica.hh"
+#include "G4GlobalMagFieldMessenger.hh"
+#include "G4AutoDelete.hh"
+
+#include "G4RotationMatrix.hh"
+
+#include "G4SDManager.hh"
+
+#include "G4VisAttributes.hh"
+#include "G4Colour.hh"
+
+#include "G4PhysicalConstants.hh"
+#include "G4SystemOfUnits.hh"
+
+#include <vector>
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4ThreadLocal 
+G4GlobalMagFieldMessenger* CloverDetectorConstruction::fMagFieldMessenger = 0; 
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverDetectorConstruction::CloverDetectorConstruction()
+ : G4VUserDetectorConstruction(),
+   fNumOfCrystal(0),
+   fLogicCrystal(NULL), fCrystalMaterial(NULL),
+   fCheckOverlaps(true)
+{
+  fNumOfCrystal = Crystal_Num;  //also need to change the fNDet in CloverEventAction.cc
+  fLogicCrystal = new G4LogicalVolume* [fNumOfCrystal];
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+CloverDetectorConstruction::~CloverDetectorConstruction()
+{
+  delete [] fLogicCrystal;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4VPhysicalVolume* CloverDetectorConstruction::Construct()
+{
+  // Define materials
+  DefineMaterials();
+  
+  // Define volumes
+  return DefineVolumes();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverDetectorConstruction::DefineMaterials()
+{ 
+  // Lead material defined using NIST Manager
+  auto nistManager = G4NistManager::Instance();
+  
+  // Vacuum
+  new G4Material("Galactic", 1., 1.01*g/mole, universe_mean_density, kStateGas, 2.73*kelvin, 3.e-18*pascal);
+
+  // Air
+  nistManager->FindOrBuildMaterial("G4_AIR");
+
+  // Iron
+  nistManager->FindOrBuildMaterial("G4_Fe");
+
+  // Aluminium
+  nistManager->FindOrBuildMaterial("G4_Al");
+
+  // Ge
+  fCrystalMaterial = nistManager->FindOrBuildMaterial("G4_Ge");
+
+  // Print materials
+  G4cout << *(G4Material::GetMaterialTable()) << G4endl;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4VPhysicalVolume* CloverDetectorConstruction::DefineVolumes()
+{
+  // Geometry parameters
+  G4double crystalLength = 8.*cm;
+  G4double crystalRadius = 25.*mm;
+  G4double crystalZPos = Radius; //14.*cm + crystalLength/2.;
+  G4double cutXY = 46.0 * mm;
+
+  G4double worldSizeZ  = 2*(crystalLength + crystalZPos);
+  G4double worldSizeXY = worldSizeZ;
+/*
+  //pipe
+  G4double pipeOuterRadius     = 100/2.* mm;
+  G4double pipeWallThickness   =  5 * mm;
+  G4double pipeLength          = 70* mm;
+  G4double pipeZPos            = 0.*cm;
+*/
+
+  //spherical chamber
+  G4double chamberOuterRadius     = 150/2.* mm;
+  G4double chamberWallThickness   =  5 * mm;
+  G4double chamberInnerRadius      = chamberOuterRadius - chamberWallThickness;
+  G4double chamberZPos            = 0.*cm;
+
+
+  //Clover casing
+  G4double caseXYInner = 2*(cutXY + 1.0) * mm;
+  G4double caseZInner  = crystalLength/2 + 4 * mm + crystalLength + 15*mm;
+  G4double caseXYWallThickness = 3.0 * mm;
+  G4double caseZWallThickness = 3.0 * mm;
+  
+  // Get materials
+  auto backgroundMaterial = G4Material::GetMaterial("G4_AIR");
+  //auto backgroundMaterial = G4Material::GetMaterial("Galactic");
+   
+  // World
+  auto worldS = new G4Box("World",  worldSizeXY/2, worldSizeXY/2, worldSizeZ); 
+                         
+  auto worldLV = new G4LogicalVolume( worldS,  backgroundMaterial, "World"); 
+                                   
+  auto worldPV = new G4PVPlacement( 0,                // no rotation
+                                    G4ThreeVector(),  // at (0,0,0)
+                                    worldLV,          // its logical volume                         
+                                    "World",          // its name
+                                    0,                // its mother  volume
+                                    false,            // no boolean operation
+                                    0,                // copy number
+                                    fCheckOverlaps);  // checking overlaps
+
+  worldLV->SetVisAttributes (G4VisAttributes::GetInvisible());
+
+  // vaccum Pipe
+  
+ // G4Tubs * pipe = new G4Tubs("pipe", pipeOuterRadius - pipeWallThickness , pipeOuterRadius, pipeLength, 0, 360*degree);
+//  G4Tubs * pipe0 = new G4Tubs("pipe0", 0 , pipeOuterRadius, pipeLength, 0, 360*degree);
+ 
+  G4Sphere *chamber = new G4Sphere("chamber",chamberInnerRadius,chamberOuterRadius,0,2*pi,0,pi);
+
+  G4RotationMatrix * rotX = new G4RotationMatrix(); rotX->rotateX(90*degree);
+  G4RotationMatrix * rotY = new G4RotationMatrix(); rotY->rotateY(90*degree);
+  
+  //G4SubtractionSolid * pipe1 = new G4SubtractionSolid("Pipe1", pipe, pipe0, rotY, G4ThreeVector());
+
+ // G4LogicalVolume * pipe1LV = new G4LogicalVolume( pipe1, G4Material::GetMaterial("G4_Fe"), "Pipe");
+ // G4LogicalVolume * pipe2LV = new G4LogicalVolume( pipe, G4Material::GetMaterial("G4_Fe"), "Pipe2");
+  G4LogicalVolume * chamberLV = new G4LogicalVolume(chamber, G4Material::GetMaterial("G4_Fe"),"targetchamber");
+  
+  new G4PVPlacement( 0,                // no rotation
+                     G4ThreeVector(0, 0, chamberZPos),  // at (0,0,0)
+                     chamberLV,          // its logical volume                         
+                     "targetchamber",          // its name
+                     worldLV,          // its mother  volume
+                     false,            // no boolean operation
+                     0,                // copy number
+                     fCheckOverlaps);  // checking overlaps
+  /*                   
+  new G4PVPlacement( rotY,                // no rotation
+                     G4ThreeVector(0, 0, pipeZPos),  // at (0,0,0)
+                     pipe2LV,          // its logical volume                         
+                     "Pipe2",          // its name
+                     worldLV,          // its mother  volume
+                     false,            // no boolean operation
+                     0,                // copy number
+                     fCheckOverlaps);  // checking overlaps
+  */
+ // pipe1LV->SetVisAttributes(new G4VisAttributes(G4Colour(1.0,1.0,1.0)));
+ // pipe2LV->SetVisAttributes(new G4VisAttributes(G4Colour(1.0,1.0,1.0)));
+    chamberLV->SetVisAttributes(new G4VisAttributes(G4Colour(1.0,1.0,1.0)));
+  
+
+  // Crystals
+  G4VisAttributes * crystalVisAtt= new G4VisAttributes(G4Colour(0.5,0.5,1.0));
+  crystalVisAtt->SetVisibility(true);
+
+//define crystal shapes
+
+  G4Tubs * base = new G4Tubs("base",  0, crystalRadius, crystalLength/2., 0, 360*degree);
+  G4Box * cut = new G4Box("cut", cutXY/2. , cutXY/2. , crystalLength/2. * 1.2);
+  G4IntersectionSolid * crystalS = new G4IntersectionSolid("HPGe", base, cut);
+
+  G4Box* case1 = new G4Box("case1", caseXYInner/2., caseXYInner/2., caseZInner  );
+  G4Box* case2 = new G4Box("case2", (caseXYInner + caseXYWallThickness)/2. , (caseXYInner + caseXYWallThickness)/2., caseZInner + caseZWallThickness*2.);
+  G4SubtractionSolid * casing = new G4SubtractionSolid("casing", case2, case1);
+  
+
+
+  double dTheta,dPsi,zmod; 
+
+  for( G4int i = 0 ; i < fNumOfCrystal ; i++){ 
+    G4String name = "HPGe"+ std::to_string(i);
+    G4cout << " crystal name : " << name << G4endl;
+    G4double phi = 360./4. * degree;
+
+    G4double rho = 0;
+    rho = (cutXY/2  + 0.3* mm) / sin(phi/2.);
+
+    int iClover  = i/4; //How to get each group of crystals.
+    int jCrystal = i%4; //How to renumber each group as 0,1,2,3. 
+    
+    fLogicCrystal[i] = new G4LogicalVolume(crystalS, fCrystalMaterial, "CrystalLV");
+
+
+    auto caseLV =  new G4LogicalVolume * [fNumOfCrystal/4];
+   // auto caseLV = new G4LogicalVolume(casing, G4Material::GetMaterial("G4_Al"), "Case")[i/4+1];
+
+    fLogicCrystal[i]->SetVisAttributes(crystalVisAtt);
+
+        
+    if(int j = iClover){    
+        
+        dTheta = Clover_Location[j][1];
+
+        dPsi = Clover_Location[j][2];
+
+        zmod = Clover_Location[j][3];
+
+        }
+        
+
+    G4ThreeVector pos = G4ThreeVector(rho * cos(jCrystal* phi + phi/2.), rho * sin(jCrystal * phi + phi/2.), crystalZPos+zmod);
+    
+//    if(iClover != DetNum){return 0;}
+
+
+    G4ThreeVector case_pos = G4ThreeVector(0,0,100*mm+crystalZPos + zmod); //The crystalZpos is the radius of the detectors. 
+ 
+    //Set rotation matrix settings for clovers.
+    
+    G4RotationMatrix * rotmat = new G4RotationMatrix(); 
+    G4RotationMatrix * crystalmat= new G4RotationMatrix();
+
+    //x = sin(theta)cos(phi) || y = sin(theta)sin(phi) || z = cos(theta) 
+   
+  
+    for(int j = 0; j <= iClover;j++){
+        if(j == iClover){
+  //  rotmat->rotateY(yangc *degree);
+  
+    rotmat->rotateZ(-dPsi*degree);
+    rotmat->rotateX(-dTheta*degree);
+   
+  //  crystalmat->rotateY(yang*degree);
+    crystalmat->rotateZ(-dPsi*degree);
+    crystalmat->rotateX(-dTheta*degree);
+   
+    //rotate position vector based on iClover.
+
+  //  pos.rotateY(dPhi*degree);
+    pos.rotateX(dTheta*degree);
+    pos.rotateZ(dPsi *degree);
+
+    //set rotation for the detector case vector; 
+
+  //  case_pos.rotateY(dPhi*degree); 
+    case_pos.rotateX(dTheta *degree); 
+    case_pos.rotateZ(dPsi*degree); 
+        }
+    }
+    new G4PVPlacement( crystalmat,     // orientation matrix
+                       pos,              // its position vector
+                       fLogicCrystal[i],        // its logical volume                         
+                       name,             // its name
+                       worldLV,          // its mother  volume
+                       false,            // no boolean operation
+                       i,                // copy number
+                       fCheckOverlaps);               // checking overlaps
+     
+    
+  if( jCrystal  == 0 ) {  
+      caseLV[iClover]  = new G4LogicalVolume(casing, G4Material::GetMaterial("G4_Al"), "Case");
+  new G4PVPlacement( rotmat ,                // case orientation
+                     case_pos,          // case location
+                     caseLV[i/4],           // its logical volume                         
+                     "Casing",         // its name
+                     worldLV,          // its mother  volume
+                     false,            // no boolean operation
+                     iClover,                // copy number
+                     fCheckOverlaps);  // checking overlaps
+  
+  caseLV[iClover]->SetVisAttributes(new G4VisAttributes(G4Colour(1.0,1.0,0.0)));
+    }
+       //Set a i+1%4 for each 
+     //
+
+  }
+
+  // Al casing
+  
+  
+  
+  // Always return the physical World
+  return worldPV;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverDetectorConstruction::ConstructSDandField()
+{
+
+  G4cout << "********************* CloverDetectorConstruction::ConstructSDandField() " << G4endl;
+  G4SDManager::GetSDMpointer()->SetVerboseLevel(1);
+
+  // Sensitive detectors
+  CloverCrystalSD * crystalSD = new CloverCrystalSD("Crystal", "CrystalHitsCollection", fNumOfCrystal);
+  G4SDManager::GetSDMpointer()->AddNewDetector(crystalSD);
+
+  //Set crystalSD to all logical volumn under the same name of "CrystalLV"
+  SetSensitiveDetector("CrystalLV",crystalSD, true);
+
+  // Create global magnetic field messenger.
+  // Uniform magnetic field is then created automatically if
+  // the field value is not zero.
+  G4ThreeVector fieldValue;
+  fMagFieldMessenger = new G4GlobalMagFieldMessenger(fieldValue);
+  fMagFieldMessenger->SetVerboseLevel(1);
+  
+  // Register the field messenger for deleting
+  G4AutoDelete::Register(fMagFieldMessenger);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

CloverDetectorConstruction.hh

@@ -0,0 +1,45 @@
+#ifndef CloverDetectorConstruction_h
+#define CloverDetectorConstruction_h 1
+
+#include "G4VUserDetectorConstruction.hh"
+#include "globals.hh"
+
+class G4VPhysicalVolume;
+class G4LogicalVolume;
+class G4Material;
+class G4GlobalMagFieldMessenger;
+
+
+class CloverDetectorConstruction : public G4VUserDetectorConstruction
+{
+  public:
+    CloverDetectorConstruction();
+    virtual ~CloverDetectorConstruction();
+
+  public:
+    virtual G4VPhysicalVolume* Construct();
+    virtual void ConstructSDandField();
+     
+  private:
+    // methods
+    void DefineMaterials();
+    G4VPhysicalVolume* DefineVolumes();
+
+    // data members
+    G4int   fNumOfCrystal;     // number of Crystal
+
+    G4LogicalVolume**  fLogicCrystal;    // pointer to the logical Crystal
+
+    G4Material*        fCrystalMaterial; // pointer to the crystal material
+    
+    // other things
+    static G4ThreadLocal G4GlobalMagFieldMessenger*  fMagFieldMessenger;  // magnetic field messenger
+
+    G4bool  fCheckOverlaps; // option to activate checking of volumes overlaps
+
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
+

CloverEventAction.cc

@@ -0,0 +1,169 @@
+
+#include "CloverEventAction.hh"
+#include "CloverCrystalSD.hh"
+#include "CloverCrystalHit.hh"
+#include "CloverAnalysis.hh"
+#include "Detector_Parameters.hh"
+
+
+#include "G4RunManager.hh"
+#include "G4Event.hh"
+#include "G4SDManager.hh"
+#include "G4HCofThisEvent.hh"
+#include "G4UnitsTable.hh"
+#include "G4SystemOfUnits.hh"
+
+//#include <cmath.h>
+#include <stdlib.h>
+
+#include "Randomize.hh"
+#include <iomanip>
+#include <vector>
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverEventAction::CloverEventAction()
+ : G4UserEventAction(),
+   fCrystalHCID(-1),
+   fNDet(Crystal_Num) // number of crystal
+{
+  fdEList.clear();
+  fdLList.clear();
+  for(int i = 0 ; i < fNDet; i++){
+    fdEList.push_back(0.);
+    fdLList.push_back(0.);
+  }
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverEventAction::~CloverEventAction()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverCrystalHitsCollection* CloverEventAction::GetHitsCollection(G4int hcID, const G4Event* event) const
+{
+  auto hitsCollection  = static_cast<CloverCrystalHitsCollection*>(event->GetHCofThisEvent()->GetHC(hcID));
+  
+  if ( ! hitsCollection ) {
+    G4ExceptionDescription msg;
+    msg << "Cannot access hitsCollection ID " << hcID; 
+    G4Exception("CloverEventAction::GetHitsCollection()",
+      "MyCode0003", FatalException, msg);
+  }         
+
+  return hitsCollection;
+}    
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverEventAction::PrintEventStatistics( G4double absoEdep, G4double absoTrackLength) const
+{
+  // print event statistics
+  G4cout
+     << "CloverEventAction::PrintEventStatistics" << G4endl
+     << "   Absorber: total energy: " 
+     << std::setw(7) << G4BestUnit(absoEdep, "Energy")
+     << "       total track length: " 
+     << std::setw(7) << G4BestUnit(absoTrackLength, "Length")
+     << G4endl;
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverEventAction::BeginOfEventAction(const G4Event* /*event*/)
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverEventAction::EndOfEventAction(const G4Event* event)
+{
+
+  G4TrajectoryContainer* trajectoryContainer = event->GetTrajectoryContainer();
+  G4int n_trajectories = 0;
+  if (trajectoryContainer) n_trajectories = trajectoryContainer->entries();
+
+  // periodic printing
+  //G4cout << "######################################### CloverEventAction::EndOfEventAction" << G4endl;
+  //G4int eventID = event->GetEventID();
+  //if ( eventID < 100 || eventID % 1000 == 0) {
+  //  G4cout << ">>>>>>>> Event: " << eventID  << G4endl;
+  //  if ( trajectoryContainer ) {
+  //    G4cout << "    " << n_trajectories
+  //           << " trajectories stored in this event." << G4endl;
+  //  }
+  //  G4VHitsCollection* hc = event->GetHCofThisEvent()->GetHC(0);
+  //  G4cout << "    "  
+  //         << hc->GetSize() << " hits stored in this event" << G4endl;
+  //}
+
+  G4PrimaryVertex * pv = event->GetPrimaryVertex();
+  G4PrimaryParticle * pp = pv->GetPrimary();
+
+  G4double beamEnergy = pp->GetKineticEnergy() * 1000; // to keV
+
+  G4ThreeVector pmomt = pp->GetMomentumDirection();
+  G4double beamTheta = pmomt.theta() / degree;
+  G4double beamPhi = pmomt.phi() / degree;
+  
+  //G4cout << "Beam Energy : " << beamEnergy << G4endl;
+  //       << "Beam theta : " << beamTheta << G4endl
+  //       << "Beam phi   : " << beamPhi << G4endl;
+    
+  // Get hits collections IDs (only once)
+  if ( fCrystalHCID == -1 ) {
+    fCrystalHCID = G4SDManager::GetSDMpointer()->GetCollectionID("CrystalHitsCollection");
+  }
+
+  // Get hits collections
+  auto crystalHC = GetHitsCollection(fCrystalHCID, event); //this is G4VHitsCollection
+
+  for( int i = 0 ; i < fNDet ; i++){
+    CloverCrystalHit * crystalHit = (*crystalHC)[i];  //this is CloverCrystalHit :: G4VHit
+
+    //add detector resolution of 1%
+
+    G4double edep = crystalHit->GetEdep();
+
+ 
+     G4double resol = G4RandGauss::shoot(1, 0.001);
+  
+    //Energy Resolution of the detectors is Energy dependent. 
+/*
+    double Fano = 0.112;
+    double electron_hole_work_function_GE = .003; //3 ev in HPGe. 
+                                                //
+                                                   
+    double energy_width = 2*std::pow((2* 0.693 * Fano * edep * electron_hole_work_function_GE),0.5);
+
+    double electron_width = 0.001;
+
+    double e_res =  std::pow( std::pow(energy_width,2) + std::pow(electron_width,2),0.5);
+
+ 
+    G4double resol = G4RandGauss::shoot(1,e_res);
+  */  
+    fdEList[i] = edep * resol * 1000;  // to keV
+    fdLList[i] = crystalHit->GetStepLength();
+
+  }
+
+  
+  // get analysis manager
+  auto analysisManager = G4AnalysisManager::Instance();
+ 
+  // fill ntuple
+  analysisManager->FillNtupleIColumn(0, n_trajectories);
+
+  analysisManager->FillNtupleDColumn(3, beamEnergy);
+  analysisManager->FillNtupleDColumn(4, beamTheta);
+  analysisManager->FillNtupleDColumn(5, beamPhi);
+  
+  analysisManager->AddNtupleRow();
+
+}  
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

CloverEventAction.hh

@@ -0,0 +1,54 @@
+
+#ifndef CloverEventAction_h
+#define CloverEventAction_h 1
+
+#include "G4UserEventAction.hh"
+
+#include "CloverCrystalHit.hh"
+
+#include "globals.hh"
+
+#include <vector>
+
+/// Event action class
+///
+/// In EndOfEventAction(), it prints the accumulated quantities of the energy 
+/// deposit and track lengths of charged particles in Crystal 
+/// stored in the hits collections.
+
+class CloverEventAction : public G4UserEventAction
+{
+public:
+  CloverEventAction();
+  virtual ~CloverEventAction();
+
+  virtual void  BeginOfEventAction(const G4Event* event);
+  virtual void    EndOfEventAction(const G4Event* event);
+
+  std::vector<G4double>& GetdEList()  {return fdEList;};
+  std::vector<G4double>& GetStepLengthList()  {return fdLList;};
+
+
+private:
+  // methods
+  CloverCrystalHitsCollection* GetHitsCollection(G4int hcID, const G4Event* event) const;
+  void PrintEventStatistics(G4double absoEdep, G4double absoTrackLength) const;
+  
+  // data members
+  G4int  fCrystalHCID; // Hit collection ID
+
+  G4int fNDet;
+
+  std::vector<G4double> fdEList;   // dE of each crystal
+  std::vector<G4double> fdLList;   // step length of each crystal
+
+  //std::vector<G4double> fLastPost; //last tracking position
+  
+  
+};
+                     
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
+
+    

CloverPhysicsList.cc

@@ -0,0 +1,126 @@
+
+#include "CloverPhysicsList.hh"
+
+#include "G4ParticleDefinition.hh"
+#include "G4ProcessManager.hh"
+#include "G4ParticleTypes.hh"
+#include "G4ParticleTable.hh"
+
+#include "G4SystemOfUnits.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverPhysicsList::CloverPhysicsList()
+: G4VUserPhysicsList()
+{
+  defaultCutValue = 1.0*mm;
+  SetVerboseLevel(1);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverPhysicsList::~CloverPhysicsList()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverPhysicsList::ConstructParticle()
+{
+  // In this method, static member functions should be called
+  // for all particles which you want to use.
+  // This ensures that objects of these particle types will be
+  // created in the program. 
+  ConstructBosons();
+  ConstructLeptons();
+  G4GenericIon::GenericIonDefinition();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverPhysicsList::ConstructBosons()
+{
+  G4Gamma::GammaDefinition();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverPhysicsList::ConstructLeptons()
+{
+  G4Electron::ElectronDefinition();
+  G4Positron::PositronDefinition();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverPhysicsList::ConstructProcess()
+{
+  AddTransportation();
+  ConstructEM();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "G4PhysicsListHelper.hh"
+
+#include "G4ComptonScattering.hh"
+#include "G4GammaConversion.hh"
+#include "G4PhotoElectricEffect.hh"
+
+#include "G4eMultipleScattering.hh"
+
+#include "G4eIonisation.hh"
+#include "G4eBremsstrahlung.hh"
+#include "G4eplusAnnihilation.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverPhysicsList::ConstructEM()
+{
+  G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
+  
+  auto particleIterator=GetParticleIterator();
+  particleIterator->reset();
+  while( (*particleIterator)() ){
+    G4ParticleDefinition* particle = particleIterator->value();
+    G4String particleName = particle->GetParticleName();
+     
+    if (particleName == "gamma") {
+          
+      ph->RegisterProcess(new G4PhotoElectricEffect, particle);
+      ph->RegisterProcess(new G4ComptonScattering,   particle);
+      ph->RegisterProcess(new G4GammaConversion,     particle);
+      
+    } else if (particleName == "e-") {
+
+      ph->RegisterProcess(new G4eMultipleScattering, particle);
+      ph->RegisterProcess(new G4eIonisation,         particle);
+      ph->RegisterProcess(new G4eBremsstrahlung,     particle);      
+
+    } else if (particleName == "e+") {
+
+      ph->RegisterProcess(new G4eMultipleScattering, particle);
+      ph->RegisterProcess(new G4eIonisation,         particle);
+      ph->RegisterProcess(new G4eBremsstrahlung,     particle);
+      ph->RegisterProcess(new G4eplusAnnihilation,   particle);      
+    }
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverPhysicsList::SetCuts()
+{
+  if (verboseLevel >0){
+    G4cout << "PhysicsList::SetCuts:";
+    G4cout << "CutLength : " << defaultCutValue/mm << " (mm)" << G4endl;
+  }
+     
+  // set cut values for gamma at first and for e- second and next for e+,
+  // because some processes for e+/e- need cut values for gamma 
+  SetCutValue(defaultCutValue, "gamma");
+  SetCutValue(defaultCutValue, "e-");
+  SetCutValue(defaultCutValue, "e+");
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+

CloverPhysicsList.hh

@@ -0,0 +1,38 @@
+
+#ifndef CloverPhysicsList_h
+#define CloverPhysicsList_h 1
+
+#include "G4VUserPhysicsList.hh"
+#include "globals.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+class CloverPhysicsList: public G4VUserPhysicsList
+{
+  public:
+    CloverPhysicsList();
+   ~CloverPhysicsList();
+
+  protected:
+    // Construct particle and physics
+    virtual void ConstructParticle();
+    virtual void ConstructProcess();
+ 
+    virtual void SetCuts();
+    
+  protected:
+    // these methods Construct particles 
+    void ConstructBosons();
+    void ConstructLeptons();
+
+  protected:
+  // these methods Construct physics processes and register them
+    void ConstructGeneral();
+    void ConstructEM();
+    
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
+

CloverPrimaryGeneratorAction.cc

@@ -0,0 +1,183 @@
+
+#include "CloverPrimaryGeneratorAction.hh"
+#include "CloverEventAction.hh"
+
+#include "G4RunManager.hh"
+#include "G4LogicalVolumeStore.hh"
+#include "G4LogicalVolume.hh"
+#include "G4Box.hh"
+#include "G4Event.hh"
+#include "G4ParticleGun.hh"
+#include "G4ParticleTable.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4GenericMessenger.hh"
+#include "Randomize.hh"
+
+#include <vector>
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverPrimaryGeneratorAction::CloverPrimaryGeneratorAction()
+ : G4VUserPrimaryGeneratorAction(),
+   fMessenger(nullptr),
+   fParticleGun(nullptr),
+   fAngle(180 * degree)
+   //,fEventAction(eventAction)
+{
+  G4int nofParticles = 1;
+  fParticleGun = new G4ParticleGun(nofParticles);
+
+  // default particle kinematic
+  //
+  auto particleDefinition = G4ParticleTable::GetParticleTable()->FindParticle("gamma");
+  fParticleGun->SetParticleDefinition(particleDefinition);
+  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
+  fParticleGun->SetParticleEnergy(6.13*MeV);
+  fParticleGun->SetParticlePosition(G4ThreeVector(0., 0., 0));
+
+
+  //Cs137
+ //fEnergyList = {{.663,100}};
+ 
+  //Co60
+
+  fEnergyList = {{1.173,99.973},
+                {1.332,99.986}};
+
+  //N16 iso
+ /*
+  fEnergyList = {{0.120, 100}, //energy, branch
+                {1.775, 0.114},
+                {1.955, 0.036},
+                {2.742, 0.777},
+                {8.872, 0.072},
+                {7.117, 5},
+                {6.917, 0.036},
+                {6.130, 68.0+0.777}};
+   */ 
+   //Eu152 
+ /*
+   fEnergyList = {{.12178,28.58},
+                  {.14801,0.037},
+                  {.21257,0.019},
+                  {.24469,7.583},
+                  {.29594,0.477},
+                  {.32913,0.128},
+                  {.34428,26.50},
+                  {.36778,0.861},
+                  {.41112,2.234},
+                  {.44396,0.327},
+                  {.44397,2.821},
+                  {.48868,0.419},
+                  {.50347,0.149},
+                  {.56399,0.489},
+                  {.56643,0.129},
+                  {.58627,0.459},
+                  {.65648,0.145},
+                  {.67468,0.172},
+                  {.67862,0.471},
+                  {.68867,0.857},
+                  {.71935,0.321},
+                  {.76490,0.215},
+                  {.77891,12.942},
+                  {.81045,0.320},
+                  {.86738,4.245},
+                  {.91933,0.427},
+                  {.92632,0.278},
+                  {.96339,0.135},
+                  {.96408,14.605},
+                  {1.00527,0.646},
+                  {1.0840,0.246},
+                  {1.08586,10.207},
+                  {1.089737,1.727},
+                  {1.109174,0.186},
+                  {1.11207,13.644},
+                  {1.21294,1.422},
+                  {1.249938,0.188},
+                  {1.292778,0.105},
+                  {1.29914,1.623},
+                  {1.408006,21.005},
+                  {1.457643,0.502},
+                  {1.528103,0.281}};
+
+*/
+
+  //accumulate the branch
+  G4int numEnergy = (G4int) fEnergyList.size();
+  for( int i = 1 ; i < numEnergy ; i++){
+    fEnergyList[i][1] += fEnergyList[i-1][1]; 
+  }
+  //renormalized the branching ratio
+  for( int i = 0 ; i < numEnergy ; i++){
+    fEnergyList[i][1] = fEnergyList[i][1]/fEnergyList[numEnergy-1][1];
+  }
+
+  //Define command for this class
+  DefineCommands();
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverPrimaryGeneratorAction::~CloverPrimaryGeneratorAction()
+{
+  delete fParticleGun;
+  delete fMessenger;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
+{
+
+  G4double cosTheta = 1.- (1.-cos(fAngle ))*G4UniformRand();
+
+  G4double phi = 360*degree*G4UniformRand();
+
+  G4double sinTheta = std::sqrt(1. - cosTheta*cosTheta);
+  G4double ux = sinTheta*std::cos(phi),
+           uy = sinTheta*std::sin(phi),
+           uz = cosTheta;
+  
+  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(ux,uy,uz));
+
+  fParticleGun->GeneratePrimaryVertex(anEvent);
+
+
+  //generate a random number form 0 to 1
+  G4double pop = G4UniformRand();
+  
+  G4double beamEnergy = 6.13;
+  
+  G4int numEnergy = (G4int) fEnergyList.size();
+  for( int i = 0; i < numEnergy  ; i++){
+    if( i == 0 &&   pop <= fEnergyList[i][1] ) beamEnergy = fEnergyList[0][0] ;
+    if( i > 0 &&  fEnergyList[i-1][1] < pop && pop <= fEnergyList[i][1] ) beamEnergy = fEnergyList[i][0] ;
+  }
+  
+  fParticleGun->SetParticleEnergy(beamEnergy);
+
+  //G4double energy = fParticleGun->GetParticleEnergy();
+ 
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+
+void CloverPrimaryGeneratorAction::DefineCommands()
+{
+  // Define /Clover/generator command directory using generic messenger class
+  fMessenger = new G4GenericMessenger(this, "/Clover/generator/", "Primary generator control");
+
+  // OpeningAngle command
+  auto& openingAngleCmd = fMessenger->DeclarePropertyWithUnit("OpeningAngle", "deg", fAngle, "Opening Angle of beam in deg");
+  openingAngleCmd.SetParameterName("t", true);
+  openingAngleCmd.SetRange("t>=0.");
+  openingAngleCmd.SetDefaultValue("180.");
+
+  
+
+}
+

CloverPrimaryGeneratorAction.hh

@@ -0,0 +1,42 @@
+
+#ifndef CloverPrimaryGeneratorAction_h
+#define CloverPrimaryGeneratorAction_h 1
+
+#include "G4VUserPrimaryGeneratorAction.hh"
+#include "globals.hh"
+
+#include <vector>
+
+class G4ParticleGun;
+class G4Event;
+class G4GenericMessenger;
+
+/// The primary generator action class with particle gum.
+
+class CloverPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
+{
+public:
+
+  CloverPrimaryGeneratorAction();    
+  virtual ~CloverPrimaryGeneratorAction();
+
+  virtual void GeneratePrimaries(G4Event* event);
+  
+  // set methods
+  void SetRandomFlag(G4bool value);
+
+private:
+  void DefineCommands();
+  G4GenericMessenger* fMessenger;
+
+  G4ParticleGun*  fParticleGun; // G4 particle gun
+
+  G4double fAngle ; 
+
+  std::vector<std::vector<G4double>> fEnergyList;
+
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif

CloverRunAction.cc

@@ -0,0 +1,87 @@
+
+#include "CloverRunAction.hh"
+#include "CloverAnalysis.hh"
+#include "CloverEventAction.hh"
+
+#include "G4Run.hh"
+#include "G4RunManager.hh"
+#include "G4UnitsTable.hh"
+#include "G4SystemOfUnits.hh"
+
+#include <vector>
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverRunAction::CloverRunAction(CloverEventAction * eventAction)
+//CloverRunAction::CloverRunAction()
+ : G4UserRunAction()
+ ,fEventAction(eventAction)
+{ 
+  // set printing event number per each event
+  G4RunManager::GetRunManager()->SetPrintProgress(1);     
+
+  // Create analysis manager
+  // The choice of analysis technology is done via selectin of a namespace
+  // in CloverAnalysis.hh
+  auto analysisManager = G4AnalysisManager::Instance();
+  G4cout << "Using " << analysisManager->GetType() << G4endl;
+
+  // Create directories 
+  analysisManager->SetVerboseLevel(1);
+  analysisManager->SetNtupleMerging(true);
+  
+  // Creating ntuple
+  //
+  analysisManager->CreateNtuple("Clover", "Edep and TrackL");
+  analysisManager->CreateNtupleIColumn("nTraj");
+  analysisManager->CreateNtupleDColumn("energy", fEventAction->GetdEList() );
+  analysisManager->CreateNtupleDColumn("stepLength", fEventAction->GetStepLengthList() );
+
+  // event inital angle, energy
+  analysisManager->CreateNtupleDColumn("bEnergy");
+  analysisManager->CreateNtupleDColumn("theta");
+  analysisManager->CreateNtupleDColumn("phi");
+  
+  analysisManager->FinishNtuple();
+   
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+CloverRunAction::~CloverRunAction()
+{
+  delete G4AnalysisManager::Instance();  
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverRunAction::BeginOfRunAction(const G4Run* /*run*/)
+{ 
+  //inform the runManager to save random number seed
+  //G4RunManager::GetRunManager()->SetRandomNumberStore(true);
+  
+  // Get analysis manager
+  auto analysisManager = G4AnalysisManager::Instance();
+
+  // Open an output file
+  //G4String fileName = "Clover.root";
+  //analysisManager->OpenFile(fileName);
+  analysisManager->OpenFile();
+ 
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void CloverRunAction::EndOfRunAction(const G4Run* /*run*/)
+{
+  
+  auto analysisManager = G4AnalysisManager::Instance();
+
+  // save histograms & ntuple
+  //
+  analysisManager->Write();
+  analysisManager->CloseFile();
+   
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

CloverRunAction.hh

@@ -0,0 +1,29 @@
+
+#ifndef CloverRunAction_h
+#define CloverRunAction_h 1
+
+#include "G4UserRunAction.hh"
+#include "globals.hh"
+
+class CloverEventAction;
+
+class G4Run;
+
+
+class CloverRunAction : public G4UserRunAction
+{
+  public:
+    CloverRunAction(CloverEventAction * eventAction);
+    virtual ~CloverRunAction();
+
+    virtual void BeginOfRunAction(const G4Run*);
+    virtual void   EndOfRunAction(const G4Run*);
+
+  private:
+    CloverEventAction * fEventAction;
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
+

Detector_Parameters.hh

@@ -0,0 +1,42 @@
+//Global Variables. 
+
+#ifndef __iamauniqueid_h__
+#define __iamauniqueid_h__
+#include <stdio.h>
+#include <cstdio>
+#include <cstdlib>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+
+   //location data array for each clover needs to be set here;
+   // clover number , theta, phi, relative z-distance to the Radius
+   //
+  
+  const double Clover_Location[16][4] = {{0,131.75,326.96,0.},
+                                         {1,150.,243.57,0.},
+                                         {2,90.0,257.14,0.},
+                                         {3,90.0,205.71,0.},
+                                         {4,48.25,326.68,0.},
+                                         {5,48.25,252.64,0.},
+                                         {6,131.75,33.31,0.},
+                                         {7,150.,116.42,0.},
+                                         {8,90.,102.85,0.},
+                                         {9,90.,154.28,0.},
+                                         {10,48.25,33.31,0.},
+                                         {11,90.,308.57,0.},
+                                         {12,90.,51.42,0.},
+                                         {13,48.25,180.,0.},
+                                         {14,129.,180.,0.},
+                                         {15,48.25,107.35,0.}};//,
+  
+
+
+
+  const int Crystal_Num = 64; // has to be divisible by 4. 
+                         
+  const double Radius = 200.; //Z POSITION -- radius away from target, in milimeters.  
+
+
+
+#endif

analysis.C

@@ -0,0 +1,38 @@
+{
+
+  TFile * file = new TFile("000Everything.root");
+
+  TTree * tree = (TTree *) file->FindObjectAny("Clover");
+
+  TH1F * h0 = new TH1F("h0", "Crystal-0; Energy[keV]; count / 1 keV", 7000,  1000, 8000); h0->SetLineColor(2);
+  TH1F * h1 = new TH1F("h1", "Crystal-1; Energy[keV]; count / 1 keV", 7000,  1000, 8000); h1->SetLineColor(4);
+  TH1F * h2 = new TH1F("h2", "Crystal-2; Energy[keV]; count / 1 keV", 7000,  1000, 8000); h2->SetLineColor(6);
+  TH1F * h3 = new TH1F("h3", "Crystal-3; Energy[keV]; count / 1 keV", 7000,  1000, 8000); h3->SetLineColor(7);
+
+  TCanvas * c1 = new TCanvas("c1", "Clover plots", 1200, 400);
+  //c1->Divide(1,4);
+
+  /*c1->cd(1);*/ tree->Draw("energy[0]>>h0", "energy[0]>0");
+  /*c1->cd(2);*/ tree->Draw("energy[1]>>h1", "energy[1]>0");
+  /*c1->cd(3);*/ tree->Draw("energy[2]>>h2", "energy[2]>0");
+  /*c1->cd(4);*/ tree->Draw("energy[3]>>h3", "energy[3]>0");
+
+  //find maximum height
+  double yMax = 0;
+
+  if( h0->GetMaximum() > yMax) yMax = h0->GetMaximum();
+  if( h1->GetMaximum() > yMax) yMax = h1->GetMaximum();
+  if( h2->GetMaximum() > yMax) yMax = h2->GetMaximum();
+  if( h3->GetMaximum() > yMax) yMax = h3->GetMaximum();
+
+  h0->SetMaximum( yMax * 1.2);
+  h1->SetMaximum( yMax * 1.2);
+  h2->SetMaximum( yMax * 1.2);
+  h3->SetMaximum( yMax * 1.2);
+
+  h0->Draw("");
+  h1->Draw("same");
+  h2->Draw("same");
+  h3->Draw("same");
+
+}

clover.cc

@@ -0,0 +1,80 @@
+
+#include "CloverDetectorConstruction.hh"
+#include "CloverActionInitialization.hh"
+#include "CloverPhysicsList.hh"
+
+#include "G4RunManagerFactory.hh"
+
+#include "G4UImanager.hh"
+#include "FTFP_BERT.hh"
+
+#include "Randomize.hh"
+
+#include "G4VisExecutive.hh"
+#include "G4UIExecutive.hh"
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+int main(int argc,char** argv)
+{
+  
+  // Detect interactive mode (if no macro provided) and define UI session
+  //
+  G4UIExecutive* ui = 0;
+  if ( argc == 1 ) {
+    ui = new G4UIExecutive(argc, argv);
+  }
+
+  // Optionally: choose a different Random engine...
+  // G4Random::setTheEngine(new CLHEP::MTwistEngine);
+  
+  // Construct the default run manager
+  auto* runManager = G4RunManagerFactory::CreateRunManager(G4RunManagerType::Default);
+
+  // Create the world and set the B-field
+  runManager->SetUserInitialization(new CloverDetectorConstruction());
+
+  // use build in physics
+  G4VModularPhysicsList* physicsList = new FTFP_BERT;
+  runManager->SetUserInitialization(physicsList);
+
+  //or use custom gamma physics
+  //runManager->SetUserInitialization(new CloverPhysicsList);
+
+  //Action Initialization  
+  runManager->SetUserInitialization(new CloverActionInitialization());
+  
+  // Initialize visualization
+  auto visManager = new G4VisExecutive;
+  // G4VisExecutive can take a verbosity argument - see /vis/verbose guidance.
+  // G4VisManager* visManager = new G4VisExecutive("Quiet");
+  visManager->Initialize();
+
+  // Get the pointer to the User Interface manager
+  auto UImanager = G4UImanager::GetUIpointer();
+
+  // Process macro or start UI session
+  if ( ! ui ) { 
+    // batch mode
+    G4String command = "/control/execute ";
+    G4String fileName = argv[1];
+    UImanager->ApplyCommand(command+fileName);
+  }
+  else { 
+    // interactive mode
+    UImanager->ApplyCommand("/control/execute init_vis.mac");
+    ui->SessionStart();
+    delete ui;
+  }
+
+  // Job termination
+  // Free the store: user actions, physics_list and detector_description are
+  // owned and deleted by the run manager, so they should not be deleted 
+  // in the main() program !
+
+  delete visManager;
+  delete runManager;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....

init_vis.mac

@@ -0,0 +1,16 @@
+# Macro file for the initialization of clover
+# in interactive session
+
+
+/control/saveHistory
+
+# Change the default number of threads (in multi-threaded mode)
+#/run/numberOfThreads 4
+
+
+# Initialize kernel
+/run/initialize
+
+
+# Visualization setting
+/control/execute vis.mac

my.in

@@ -0,0 +1,43 @@
+#
+# type "help" in the interactive shell for command
+#
+# endit vis.mac for visualization control
+
+# Initialize kernel, 
+/run/initialize
+/control/saveHistory
+
+
+/control/verbose 2
+/run/verbose 2
+/tracking/verbose 0
+/event/verbose 0
+/process/verbose 0
+
+#------ disable visualization
+/vis/disable
+/tracking/storeTrajectory 0
+
+#------ Sensitive detetcor verbose control
+/hits/verbose 0
+
+#-------override B-field
+/globalField/verbose 1
+/globalField/setValue 0.0 0.0 0.0 tesla
+
+#-------override particle, 
+#/gun/particle gamma
+#/gun/energy 6.13 MeV
+
+#-------Override Opening angle
+/Clover/generator/OpeningAngle  10 deg
+
+#-------Override output root file
+/analysis/setFileName 000SingleCrystal.root
+
+/run/printProgress 100000
+#/run/beamOn 15000000
+/run/beamOn 50000
+#/run/beamOn 20
+
+

vis.mac

@@ -0,0 +1,58 @@
+
+# Use this open statement to create a .wrl file suitable for
+# viewing in a VRML viewer:
+#/vis/open VRML2FILE
+
+# Use this open statement to create an OpenGL view:
+/vis/open OGL 600x600-0+0
+
+# Draw geometry:
+/vis/drawVolume
+
+
+/vis/viewer/set/auxiliaryEdge true
+/vis/viewer/set/lineSegmentsPerCircle 100
+
+
+# Draw smooth trajectories at end of event, showing trajectory points
+# as markers 2 pixels wide:
+/vis/scene/add/trajectories smooth
+/vis/modeling/trajectories/create/drawByCharge
+/vis/modeling/trajectories/drawByCharge-0/default/setDrawStepPts true
+/vis/modeling/trajectories/drawByCharge-0/default/setStepPtsSize 2
+
+# Draw hits at end of event:
+/vis/scene/add/hits
+
+# To draw only gammas:
+#/vis/filtering/trajectories/create/particleFilter
+#/vis/filtering/trajectories/particleFilter-0/add gamma
+
+# To invert the above, drawing all particles except gammas,
+# keep the above two lines but also add:
+#/vis/filtering/trajectories/particleFilter-0/invert true
+
+
+# Many other options are available with /vis/modeling and /vis/filtering.
+# For example, to select colour by particle ID:
+#/vis/modeling/trajectories/create/drawByParticleID
+#/vis/modeling/trajectories/drawByParticleID-0/default/setDrawStepPts true
+# To select or override default colours (note: e+ is blue by default):
+#/vis/modeling/trajectories/list
+#/vis/modeling/trajectories/drawByParticleID-0/set e+ yellow
+
+
+# To superimpose all of the events from a given run:
+/vis/scene/endOfEventAction accumulate
+
+
+# Re-establish auto refreshing and verbosity:
+/vis/viewer/set/autoRefresh true
+/vis/verbose warnings
+
+
+# Draw coordinate axes:
+/vis/scene/add/axes 0 0 0 0.1 m
+
+# output only detector
+#/vis/viewer/flush