BGO Fully operational

.gitignore

@@ -0,0 +1,4 @@
+*.d
+*.so
+*.swp
+*.swo

Analyzer.C

@@ -30,6 +30,8 @@
 TH1F * hEnergy[64];
 TH1F * hAddback[16];
 
+//if you need to make the energy histograms then just copy how the addback histograms are made.
+
 void Analyzer::Begin(TTree * /*tree*/)
 {
    // The Begin() function is called at the start of the query.
@@ -39,10 +41,10 @@ void Analyzer::Begin(TTree * /*tree*/)
    TString option = GetOption();
 
    for( int i = 0; i < 64; i++){
-        hEnergy[i] = new TH1F(Form("hE%02d", i), Form("hE%02d", i), 4000, 0, 4000);
+        hEnergy[i] = new TH1F(Form("hE%02d", i), Form("hE%02d", i), 8000, 0, 8000);
     
         if( i < 16  ){
-            hAddback[i] = new TH1F(Form("AE%02d",i), Form("AE%02d",i), 4000,0,4000);
+            hAddback[i] = new TH1F(Form("AE%02d",i), Form("AE%02d",i), 8000,0,8000);
         }
 
    }
@@ -79,19 +81,21 @@ Bool_t Analyzer::Process(Long64_t entry)
    //
    // The return value is currently not used.
    b_energy->GetEntry(entry);
+   b_bgoenergy->GetEntry(entry);
 
    double addback[16] ={0};
    for(int i =0; i < 64; i++){
-     if((*energy)[i] <= 0.) { continue ; }
+     if((*energy)[i] <= 0.) { continue ; }   //Energy threshold Addback selector.
      ///   printf("%2d, %f \n", i, (*energy)[i]);
      hEnergy[i]->Fill((*energy)[i]);
 
      int iclover = i/4;
-
+    
+     if( entry < 1000 && (*bgoenergy)[iclover] > 10.) printf("%d , %f, %d, %f \n", i, (*energy)[i], iclover,(*bgoenergy)[iclover] );    
+     if((*bgoenergy)[iclover] <= 0){
      addback[iclover] += (*energy)[i];
-  //   if( entry < 10 ) printf("%d , %f, %d, %f \n", i, (*energy)[i], iclover, addback[iclover]);  
+    }
    }
-   
    for( int i = 0; i < 16; i++){
      if(addback[i] <= 0.){continue;}
      hAddback[i]->Fill(addback[i]);
@@ -129,5 +133,5 @@ void Analyzer::Terminate()
    // a query. It always runs on the client, it can be used to present
    // the results graphically or save the results to file.
 
-    hAddback[14]->Draw();
+    hAddback[0]->Draw();
 }

Analyzer.h

@@ -26,6 +26,8 @@ public :
    Int_t           nTraj;
    vector<double>  *energy;
    vector<double>  *stepLength;
+   vector<double>  *bgoenergy;
+   vector<double>  *bgostepLength;
    Double_t        bEnergy;
    Double_t        theta;
    Double_t        phi;
@@ -34,6 +36,8 @@ public :
    TBranch        *b_nTraj;   //!
    TBranch        *b_energy;   //!
    TBranch        *b_stepLength;   //!
+   TBranch        *b_bgoenergy;   //!
+   TBranch        *b_bgostepLength;   //!
    TBranch        *b_bEnergy;   //!
    TBranch        *b_theta;   //!
    TBranch        *b_phi;   //!
@@ -73,6 +77,8 @@ void Analyzer::Init(TTree *tree)
    // Set object pointer
    energy = 0;
    stepLength = 0;
+   bgoenergy = 0;
+   bgostepLength = 0;
    // Set branch addresses and branch pointers
    if (!tree) return;
    fChain = tree;
@@ -81,6 +87,8 @@ void Analyzer::Init(TTree *tree)
    fChain->SetBranchAddress("nTraj", &nTraj, &b_nTraj);
    fChain->SetBranchAddress("energy", &energy, &b_energy);
    fChain->SetBranchAddress("stepLength", &stepLength, &b_stepLength);
+   fChain->SetBranchAddress("bgoE", &bgoenergy, &b_bgoenergy);
+   fChain->SetBranchAddress("bgostep", &bgostepLength, &b_bgostepLength);
    fChain->SetBranchAddress("bEnergy", &bEnergy, &b_bEnergy);
    fChain->SetBranchAddress("theta", &theta, &b_theta);
    fChain->SetBranchAddress("phi", &phi, &b_phi);

AnalyzerCloveronly.C

@@ -0,0 +1,133 @@
+#define Analyzer_cxx
+// The class definition in Analyzer.h has been generated automatically
+// by the ROOT utility TTree::MakeSelector(). This class is derived
+// from the ROOT class TSelector. For more information on the TSelector
+// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.
+
+// The following methods are defined in this file:
+//    Begin():        called every time a loop on the tree starts,
+//                    a convenient place to create your histograms.
+//    SlaveBegin():   called after Begin(), when on PROOF called only on the
+//                    slave servers.
+//    Process():      called for each event, in this function you decide what
+//                    to read and fill your histograms.
+//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
+//                    called only on the slave servers.
+//    Terminate():    called at the end of the loop on the tree,
+//                    a convenient place to draw/fit your histograms.
+//
+// To use this file, try the following session on your Tree T:
+//
+// root> T->Process("Analyzer.C")
+// root> T->Process("Analyzer.C","some options")
+// root> T->Process("Analyzer.C+")
+//
+
+#include "Analyzer.h"
+#include <TH2.h>
+#include <TStyle.h>
+
+TH1F * hEnergy[64];
+TH1F * hAddback[16];
+
+void Analyzer::Begin(TTree * /*tree*/)
+{
+   // The Begin() function is called at the start of the query.
+   // When running with PROOF Begin() is only called on the client.
+   // The tree argument is deprecated (on PROOF 0 is passed).
+
+   TString option = GetOption();
+
+   for( int i = 0; i < 64; i++){
+        hEnergy[i] = new TH1F(Form("hE%02d", i), Form("hE%02d", i), 4000, 0, 4000);
+    
+        if( i < 16  ){
+            hAddback[i] = new TH1F(Form("AE%02d",i), Form("AE%02d",i), 4000,0,4000);
+        }
+
+   }
+   
+}
+
+void Analyzer::SlaveBegin(TTree * /*tree*/)
+{
+   // The SlaveBegin() function is called after the Begin() function.
+   // When running with PROOF SlaveBegin() is called on each slave server.
+   // The tree argument is deprecated (on PROOF 0 is passed).
+
+   TString option = GetOption();
+
+}
+
+Bool_t Analyzer::Process(Long64_t entry)
+{
+   // The Process() function is called for each entry in the tree (or possibly
+   // keyed object in the case of PROOF) to be processed. The entry argument
+   // specifies which entry in the currently loaded tree is to be processed.
+   // It can be passed to either Analyzer::GetEntry() or TBranch::GetEntry()
+   // to read either all or the required parts of the data. When processing
+   // keyed objects with PROOF, the object is already loaded and is available
+   // via the fObject pointer.
+   //
+   // This function should contain the "body" of the analysis. It can contain
+   // simple or elaborate selection criteria, run algorithms on the data
+   // of the event and typically fill histograms.
+   //
+   // The processing can be stopped by calling Abort().
+   //
+   // Use fStatus to set the return value of TTree::Process().
+   //
+   // The return value is currently not used.
+   b_energy->GetEntry(entry);
+
+   double addback[16] ={0};
+   for(int i =0; i < 64; i++){
+     if((*energy)[i] <= 0.) { continue ; }
+     ///   printf("%2d, %f \n", i, (*energy)[i]);
+     hEnergy[i]->Fill((*energy)[i]);
+
+     int iclover = i/4;
+
+     addback[iclover] += (*energy)[i];
+  //   if( entry < 10 ) printf("%d , %f, %d, %f \n", i, (*energy)[i], iclover, addback[iclover]);  
+   }
+   
+   for( int i = 0; i < 16; i++){
+     if(addback[i] <= 0.){continue;}
+     hAddback[i]->Fill(addback[i]);
+
+//     if( entry < 10 ) printf("=--==%d, %f \n", i, addback[i]); 
+
+   }
+
+   return kTRUE;
+}
+
+void Analyzer::SlaveTerminate()
+{
+   // The SlaveTerminate() function is called after all entries or objects
+   // have been processed. When running with PROOF SlaveTerminate() is called
+   // on each slave server.
+
+}
+
+void Analyzer::Terminate()
+{
+   TFile * out = new TFile("CloverAddback.root","recreate");
+   out->cd();
+
+   for(int j = 0; j< 64; j++){
+       
+        hEnergy[j]->Write(); 
+
+   }
+   for(int k = 0; k< 16; k++){
+        hAddback[k]->Write();
+   }
+    out->Close();
+   // The Terminate() function is the last function to be called during
+   // a query. It always runs on the client, it can be used to present
+   // the results graphically or save the results to file.
+
+    hAddback[14]->Draw();
+}

AnalyzerCloveronly.h

@@ -0,0 +1,100 @@
+//////////////////////////////////////////////////////////
+// This class has been automatically generated on
+// Mon Oct 17 14:56:21 2022 by ROOT version 6.26/06
+// from TTree Clover/Edep and TrackL
+// found on file: AllCrystal_Co60.root
+//////////////////////////////////////////////////////////
+
+#ifndef Analyzer_h
+#define Analyzer_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TSelector.h>
+
+// Header file for the classes stored in the TTree if any.
+#include "vector"
+
+class Analyzer : public TSelector {
+public :
+   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
+
+// Fixed size dimensions of array or collections stored in the TTree if any.
+
+   // Declaration of leaf types
+   Int_t           nTraj;
+   vector<double>  *energy;
+   vector<double>  *stepLength;
+   Double_t        bEnergy;
+   Double_t        theta;
+   Double_t        phi;
+
+   // List of branches
+   TBranch        *b_nTraj;   //!
+   TBranch        *b_energy;   //!
+   TBranch        *b_stepLength;   //!
+   TBranch        *b_bEnergy;   //!
+   TBranch        *b_theta;   //!
+   TBranch        *b_phi;   //!
+
+   Analyzer(TTree * /*tree*/ =0) : fChain(0) { }
+   virtual ~Analyzer() { }
+   virtual Int_t   Version() const { return 2; }
+   virtual void    Begin(TTree *tree);
+   virtual void    SlaveBegin(TTree *tree);
+   virtual void    Init(TTree *tree);
+   virtual Bool_t  Notify();
+   virtual Bool_t  Process(Long64_t entry);
+   virtual Int_t   GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
+   virtual void    SetOption(const char *option) { fOption = option; }
+   virtual void    SetObject(TObject *obj) { fObject = obj; }
+   virtual void    SetInputList(TList *input) { fInput = input; }
+   virtual TList  *GetOutputList() const { return fOutput; }
+   virtual void    SlaveTerminate();
+   virtual void    Terminate();
+
+   ClassDef(Analyzer,0);
+};
+
+#endif
+
+#ifdef Analyzer_cxx
+void Analyzer::Init(TTree *tree)
+{
+   // The Init() function is called when the selector needs to initialize
+   // a new tree or chain. Typically here the branch addresses and branch
+   // pointers of the tree will be set.
+   // It is normally not necessary to make changes to the generated
+   // code, but the routine can be extended by the user if needed.
+   // Init() will be called many times when running on PROOF
+   // (once per file to be processed).
+
+   // Set object pointer
+   energy = 0;
+   stepLength = 0;
+   // Set branch addresses and branch pointers
+   if (!tree) return;
+   fChain = tree;
+   fChain->SetMakeClass(1);
+
+   fChain->SetBranchAddress("nTraj", &nTraj, &b_nTraj);
+   fChain->SetBranchAddress("energy", &energy, &b_energy);
+   fChain->SetBranchAddress("stepLength", &stepLength, &b_stepLength);
+   fChain->SetBranchAddress("bEnergy", &bEnergy, &b_bEnergy);
+   fChain->SetBranchAddress("theta", &theta, &b_theta);
+   fChain->SetBranchAddress("phi", &phi, &b_phi);
+}
+
+Bool_t Analyzer::Notify()
+{
+   // The Notify() function is called when a new file is opened. This
+   // can be either for a new TTree in a TChain or when when a new TTree
+   // is started when using PROOF. It is normally not necessary to make changes
+   // to the generated code, but the routine can be extended by the
+   // user if needed. The return value is currently not used.
+
+   return kTRUE;
+}
+
+#endif // #ifdef Analyzer_cxx

Analyzer_C.d

@@ -1,89 +0,0 @@
-
-# DO NOT DELETE
-
-./Analyzer_C.so: Analyzer.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TROOT.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TDirectory.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TNamed.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TObject.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/Rtypes.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/RtypesCore.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/ROOT/RConfig.hxx
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/RVersion.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/RConfigure.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/DllImport.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/strtok.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/strlcpy.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/snprintf.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TGenericClassInfo.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TSchemaHelper.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TIsAProxy.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TVirtualIsAProxy.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TStorage.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TVersionCheck.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/RVersion.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TString.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TMathBase.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/ROOT/RStringView.hxx
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/ROOT/TypeTraits.hxx
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TClass.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TDictionary.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/ESTLType.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TObjArray.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TSeqCollection.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TCollection.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TIterator.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TVirtualRWMutex.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TVirtualMutex.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/ROOT/RRangeCast.hxx
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/ROOT/RSpan.hxx
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/ROOT/span.hxx
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TUUID.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TList.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TBuffer.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TDataType.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/Bytes.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/Byteswap.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TChain.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TTree.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/Compression.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/ROOT/TIOFeatures.hxx
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TArrayD.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TArray.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TArrayI.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TAttFill.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TAttLine.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TAttMarker.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TVirtualTreePlayer.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TBranch.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TBranchCacheInfo.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TBits.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TFile.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TDirectoryFile.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TDatime.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TUrl.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/ROOT/RConcurrentHashColl.hxx
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/ROOT/TRWSpinLock.hxx
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/ROOT/TSpinMutex.hxx
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TSelector.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TSelectorList.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/THashList.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TH2.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TH1.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TAxis.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TAttAxis.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TArrayC.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TArrayS.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TArrayF.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/Foption.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/ROOT/EExecutionPolicy.hxx
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TVectorFfwd.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TVectorDfwd.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TFitResultPtr.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TMatrixFBasefwd.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TMatrixDBasefwd.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TStyle.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TAttText.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/TColor.h
-./Analyzer_C.so: /usr/local/Cellar/root/6.26.06_1/include/root/RVersion.h /usr/local/Cellar/root/6.26.06_1/include/root/RConfig.h /usr/local/Cellar/root/6.26.06_1/include/root/TClass.h /usr/local/Cellar/root/6.26.06_1/include/root/TDictAttributeMap.h /usr/local/Cellar/root/6.26.06_1/include/root/TInterpreter.h /usr/local/Cellar/root/6.26.06_1/include/root/TROOT.h /usr/local/Cellar/root/6.26.06_1/include/root/TBuffer.h /usr/local/Cellar/root/6.26.06_1/include/root/TMemberInspector.h /usr/local/Cellar/root/6.26.06_1/include/root/TError.h /usr/local/Cellar/root/6.26.06_1/include/root/RtypesImp.h /usr/local/Cellar/root/6.26.06_1/include/root/TIsAProxy.h /usr/local/Cellar/root/6.26.06_1/include/root/TFileMergeInfo.h /usr/local/Cellar/root/6.26.06_1/include/root/TCollectionProxyInfo.h /usr/local/bin/rootcling
-Analyzer_C__ROOTBUILDVERSION= 6.26/06

BGOHIT.cc

@@ -0,0 +1,81 @@
+
+#include "BGOHIT.hh"
+
+#include "G4UnitsTable.hh"
+#include "G4VVisManager.hh"
+#include "G4Circle.hh"
+#include "G4Colour.hh"
+#include "G4VisAttributes.hh"
+
+#include <iomanip>
+
+G4ThreadLocal G4Allocator<BGOHit>* BGOHitAllocator = 0;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+BGOHit::BGOHit()
+ : G4VHit(),
+   fTrackID(-1),
+   fBGOID(-1),
+   fEdep(0.),
+   fPos(G4ThreeVector()),
+   fStepLength(0.)
+{
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+BGOHit::~BGOHit() {}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+BGOHit::BGOHit(const BGOHit& right)
+  : G4VHit()
+{
+  fTrackID   = right.fTrackID;
+  fBGOID  = right.fBGOID;
+  fEdep      = right.fEdep;
+  fPos       = right.fPos;
+
+  fStepLength = right.fStepLength;
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+const BGOHit& BGOHit::operator=(const BGOHit& right)
+{
+  fTrackID   = right.fTrackID;
+  fBGOID  = right.fBGOID;
+  fEdep      = right.fEdep;
+  fPos       = right.fPos;
+
+  fStepLength = right.fStepLength;
+
+  return *this;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4bool BGOHit::operator==(const BGOHit& right) const
+{
+  return ( this == &right ) ? true : false;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void BGOHit::Print()
+{
+  G4cout
+     
+     << "  trackID: " << fTrackID << " Crystal: " << fBGOID 
+     << "  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......

BGOHIT.hh

@@ -0,0 +1,87 @@
+
+#ifndef BGOHit_h
+#define BGOHit_h 
+
+#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 BGOHit : public G4VHit
+{
+  public:
+    BGOHit();
+    BGOHit(const BGOHit&);
+    virtual ~BGOHit();
+
+    // operators
+    const BGOHit& operator=(const BGOHit&);
+    G4bool operator==(const BGOHit&) 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 SetBGOID  (G4int id)         { fBGOID = id; };
+    void SetEdep       (G4double de)      { fEdep += de; };
+    void SetPos        (G4ThreeVector xyz){ fPos = xyz; };
+    
+    void SetStepLength (G4double sl)      { fStepLength += sl; };
+
+    // Get methods
+    G4int GetTrackID() const     { return fTrackID; };
+    G4int GetBGOID() const   { return fBGOID; };
+    G4double GetEdep() const     { return fEdep; };
+    G4ThreeVector GetPos() const { return fPos; };
+
+    G4double GetStepLength () const  { return fStepLength ; };
+
+  private:
+
+    G4int         fTrackID;
+    G4int         fBGOID;
+    G4double      fEdep;
+    G4ThreeVector fPos;
+    G4double      fStepLength;
+
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+using BGOHitsCollection = G4THitsCollection<BGOHit>;
+
+extern G4ThreadLocal G4Allocator<BGOHit>* BGOHitAllocator;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+inline void* BGOHit::operator new(size_t)
+{
+  if (!BGOHitAllocator) BGOHitAllocator = new G4Allocator<BGOHit>;
+  
+  return (void *) BGOHitAllocator->MallocSingle();
+}
+
+inline void BGOHit::operator delete(void *hit)
+{
+  if (!BGOHitAllocator) BGOHitAllocator = new G4Allocator<BGOHit>;
+  
+  BGOHitAllocator->FreeSingle((BGOHit*) hit);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif

BGOSD.cc

@@ -0,0 +1,114 @@
+
+//#include "CloverCrystalSD.hh"
+#include "BGOSD.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)
+BGOSD::BGOSD(const G4String& name, const G4String& hitsCollectionName, const G4int nCrystal)
+ : G4VSensitiveDetector(name),
+   fHitsCollection(nullptr),
+   fNDet(nCrystal)
+{
+  collectionName.insert(hitsCollectionName);
+}
+
+///....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+//CloverCrystalSD::~CloverCrystalSD() 
+BGOSD::~BGOSD() 
+{ 
+}
+
+///....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+//void CloverCrystalSD::Initialize(G4HCofThisEvent* hce)
+void BGOSD::Initialize(G4HCofThisEvent* hce)
+{
+  // Create hits collection
+  //fHitsCollection = new CloverCrystalHitsCollection(SensitiveDetectorName, collectionName[0]); 
+  fHitsCollection = new BGOHitsCollection(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 BGOHit());
+  }
+
+  //G4cout << "######### size of fHitCollection : " <<  fHitsCollection->GetSize() << G4endl;
+
+}
+
+///....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+//G4bool CloverCrystalSD::ProcessHits(G4Step* step, G4TouchableHistory* /*history*/)
+G4bool BGOSD::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 BGOID = step->GetPreStepPoint()->GetTouchableHandle() ->GetCopyNumber();
+  //----------- save hit in each crystal
+  BGOHit * hit = (*fHitsCollection)[BGOID];
+
+  hit->SetTrackID  (step->GetTrack()->GetTrackID());
+  hit->SetEdep(edep);
+  hit->SetPos (step->GetPostStepPoint()->GetPosition());
+  hit->SetStepLength( stepLength);
+  hit->SetBGOID(BGOID);
+
+
+  //---------- Save indivual hit
+ // CloverCrystalHit* newHit = new CloverCrystalHit();
+  BGOHit* newHit = new BGOHit();
+  
+  newHit->SetTrackID  (step->GetTrack()->GetTrackID());
+  newHit->SetBGOID(BGOID);
+  newHit->SetEdep(edep);
+  newHit->SetPos (step->GetPostStepPoint()->GetPosition());
+  newHit->SetStepLength( stepLength);
+  
+  fHitsCollection->insert( newHit );
+
+
+  return true;
+}
+
+///....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+//void CloverCrystalSD::EndOfEvent(G4HCofThisEvent*)
+void BGOSD::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......

BGOSD.hh

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

CMakeLists.txt

@@ -49,6 +49,8 @@ set(CLOVER_SCRIPTS
   vis.mac
   my.in
   analysis.C
+  Analyzer.C
+  Analyzer.h
 )
 
 foreach(_script ${CLOVER_SCRIPTS})

CloverDetectorConstruction.cc

@@ -1,6 +1,7 @@
 
 #include "CloverDetectorConstruction.hh"
 #include "CloverCrystalSD.hh"
+#include "BGOSD.hh"
 #include "Detector_Parameters.hh"
 
 #include "G4Material.hh"
@@ -43,17 +44,19 @@ G4GlobalMagFieldMessenger* CloverDetectorConstruction::fMagFieldMessenger = 0;
 CloverDetectorConstruction::CloverDetectorConstruction()
  : G4VUserDetectorConstruction(),
    fNumOfCrystal(0),
-   fLogicCrystal(NULL), fCrystalMaterial(NULL),
+   fLogicCrystal(NULL), fLogicBGO(NULL), fCrystalMaterial(NULL),
    fCheckOverlaps(true)
 {
   fNumOfCrystal = Crystal_Num;  //also need to change the fNDet in CloverEventAction.cc
   fLogicCrystal = new G4LogicalVolume* [fNumOfCrystal];
+  fLogicBGO = new G4LogicalVolume* [fNumOfCrystal/4];
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 CloverDetectorConstruction::~CloverDetectorConstruction()
 {
   delete [] fLogicCrystal;
+  delete [] fLogicBGO;
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -86,6 +89,24 @@ void CloverDetectorConstruction::DefineMaterials()
   // Aluminium
   nistManager->FindOrBuildMaterial("G4_Al");
 
+
+
+  //nistManager->FindOrBuildMaterial("G4_Bi4-Ge3-O12");
+  
+    auto a = 208.98*g/mole;
+    G4Element* bi = new G4Element("Bismuth","Bi",83.,a);
+    auto b = 72.64*g/mole;
+    G4Element* ge = new G4Element("Germanium","Ge",32.,b);
+    auto c = 16.0*g/mole;
+    G4Element* ox = new G4Element("Oxygen","ox",16.,c);
+    auto density = 7.13*g/cm3;
+    G4Material* bgocrystalmaterial = new G4Material("bgocrystalmaterial",density,3);
+
+    bgocrystalmaterial->AddElement(bi, 4);
+    bgocrystalmaterial->AddElement(ge, 3);
+    bgocrystalmaterial->AddElement(ox, 12);
+  
+
   // Ge
   fCrystalMaterial = nistManager->FindOrBuildMaterial("G4_Ge");
 
@@ -212,16 +233,25 @@ G4VPhysicalVolume* CloverDetectorConstruction::DefineVolumes()
 
 
   G4VSolid* bgo_u = new G4UnionSolid("bgo2 + bgo3",bgo2,bgo3,bgo3_m,bgo3_v);
+  
+  G4Trd* bgoactive = new G4Trd("bgoactive", (bgoXYInner+ bgoXYWallThickness)/2.  - bgoanglecut1 - 22*mm ,(bgoXYInner+ bgoXYWallThickness)/2. - bgoanglecut1 + 12*mm
+                                , (bgoXYInner + bgoXYWallThickness)/2. - bgoanglecut1 - 22*mm, (bgoXYInner + bgoXYWallThickness)/2. - bgoanglecut1 +12*mm, zbgo1);
   //G4Box* bgo2 = new G4Box("bgo2", (bgoXYInner+ bgoXYWallThickness)/2., (bgoXYInner + bgoXYWallThickness)/2., bgoZInner + bgoZWallThickness*2.);
 
+  G4Trd* bgo33 = new G4Trd("bgo33", (bgoXYInner+ bgoXYWallThickness)/2.  - bgoanglecut0 + bgoanglecut2 - 10*mm,(bgoXYInner+ bgoXYWallThickness)/2. - bgoanglecut1
+                                , (bgoXYInner + bgoXYWallThickness)/2. - bgoanglecut0 + bgoanglecut2 - 10*mm, (bgoXYInner + bgoXYWallThickness)/2.- bgoanglecut1, zbgo1);
  
-  G4Trd* bgowindow = new G4Trd("bgowindow", (bgoXYInner+ bgoXYWallThickness)/2.  - bgoanglecut0 ,(bgoXYInner+ bgoXYWallThickness)/2. - bgoanglecut0 +bgoanglecut2
+  G4Trd* bgowindow1 = new G4Trd("bgowindow1", (bgoXYInner+ bgoXYWallThickness)/2.  - bgoanglecut0 ,(bgoXYInner+ bgoXYWallThickness)/2. - bgoanglecut0 +bgoanglecut2
                                 , (bgoXYInner + bgoXYWallThickness)/2. - bgoanglecut0, (bgoXYInner + bgoXYWallThickness)/2. -bgoanglecut0 +bgoanglecut2, zbgo2/2);
 
+  G4Box* bgo11 = new G4Box("bgo11", bgoXYInner/2. - 10*mm, bgoXYInner/2.-10*mm, bgoZInner +zbgo1);
+  G4SubtractionSolid * bgowindow = new G4SubtractionSolid("bgowindow",bgowindow1,bgo11);
 
+  G4SubtractionSolid * bgoSH= new G4SubtractionSolid("bgoSH", bgo_u, bgo1);
 
-  G4SubtractionSolid * bgo = new G4SubtractionSolid("bgo", bgo_u, bgo1);
+  G4SubtractionSolid * bgo = new G4SubtractionSolid("bgo", bgoactive, bgo33);
  
+  //G4SubtractionSolid * bgoSH= new G4SubtractionSolid("bgoSH", bgoSH0, bgo);
   // creating the subtraction solid that will form the actual bgo sheild.
 
   
@@ -243,7 +273,7 @@ G4VPhysicalVolume* CloverDetectorConstruction::DefineVolumes()
 
 
     auto caseLV =  new G4LogicalVolume * [fNumOfCrystal/4];
-    auto bgoLV  =  new G4LogicalVolume * [fNumOfCrystal/4];
+    auto bgoshieldLV  =  new G4LogicalVolume * [fNumOfCrystal/4];
     auto bgowindowLV  =  new G4LogicalVolume * [fNumOfCrystal/4];
 
     fLogicCrystal[i]->SetVisAttributes(crystalVisAtt);
@@ -268,13 +298,15 @@ G4VPhysicalVolume* CloverDetectorConstruction::DefineVolumes()
     G4ThreeVector case_pos = G4ThreeVector(0,0,88*mm+crystalZPos); //The crystalZpos is the radius of the detector,
                                                                    //bgo uses an identical rotation matrix as the case, 
                                                                    //you can also add in zmod if you want to shift the shielding. 
-    G4ThreeVector bgo_pos = G4ThreeVector(0,0,115*mm+crystalZPos);
+    G4ThreeVector bgoS_pos = G4ThreeVector(0,0,115*mm+crystalZPos);
+    G4ThreeVector bgoA_pos = G4ThreeVector(0,0,115*mm+crystalZPos - 80*mm);
     G4ThreeVector bgowindow_pos = G4ThreeVector(0,0,-63*mm + crystalZPos);
                                                                           
     //Set rotation matrix settings for clovers.
     
     G4RotationMatrix * bgowindowmat = new G4RotationMatrix(); 
-    G4RotationMatrix * bgomat = new G4RotationMatrix(); 
+    G4RotationMatrix * bgoSmat = new G4RotationMatrix(); 
+    G4RotationMatrix * bgoAmat = new G4RotationMatrix(); 
     G4RotationMatrix * rotmat = new G4RotationMatrix(); 
     G4RotationMatrix * crystalmat= new G4RotationMatrix();
 
@@ -288,9 +320,12 @@ G4VPhysicalVolume* CloverDetectorConstruction::DefineVolumes()
     bgowindowmat->rotateZ(-dPsi*degree);
     bgowindowmat->rotateX(-dTheta*degree);
 
-    bgomat->rotateZ(-dPsi*degree);
+    bgoSmat->rotateZ(-dPsi*degree);
-    bgomat->rotateX(-dTheta*degree);
+    bgoSmat->rotateX(-dTheta*degree);
    
+    bgoAmat->rotateZ(-dPsi*degree);
+    bgoAmat->rotateX(-dTheta*degree);
+    
     rotmat->rotateZ(-dPsi*degree);
     rotmat->rotateX(-dTheta*degree);
    
@@ -305,9 +340,12 @@ G4VPhysicalVolume* CloverDetectorConstruction::DefineVolumes()
     case_pos.rotateX(dTheta *degree); 
     case_pos.rotateZ(dPsi*degree); 
     
-    bgo_pos.rotateX(dTheta *degree); 
+    bgoA_pos.rotateX(dTheta *degree); 
-    bgo_pos.rotateZ(dPsi*degree); 
+    bgoA_pos.rotateZ(dPsi*degree); 
     
+    bgoS_pos.rotateX(dTheta *degree); 
+    bgoS_pos.rotateZ(dPsi*degree); 
+   
     bgowindow_pos.rotateX(dTheta *degree); 
     bgowindow_pos.rotateZ(dPsi*degree); 
         }
@@ -324,8 +362,10 @@ G4VPhysicalVolume* CloverDetectorConstruction::DefineVolumes()
     
   if( jCrystal  == 0 ) {  
       caseLV[iClover]  = new G4LogicalVolume(casing, G4Material::GetMaterial("G4_Al"), "Case");
-      bgoLV[iClover]   = new G4LogicalVolume(bgo   , G4Material::GetMaterial("G4_Al"), "BGO ");
+      fLogicBGO[iClover]   = new G4LogicalVolume(bgo   , G4Material::GetMaterial("G4_Ge"), "BGOLV");
+      //fLogicBGO[iClover]   = new G4LogicalVolume(bgo   , bgocrystalmaterial, "BGOLV");
       bgowindowLV[iClover]   = new G4LogicalVolume(bgowindow, G4Material::GetMaterial("G4_Al"), "BGO WINDOW ");
+      bgoshieldLV[iClover]   = new G4LogicalVolume(bgoSH, G4Material::GetMaterial("G4_Al"), "BGO SHIELD ");
      
       new G4PVPlacement( rotmat ,                // case orientation
                      case_pos,          // case location
@@ -336,10 +376,19 @@ G4VPhysicalVolume* CloverDetectorConstruction::DefineVolumes()
                      iClover,                // copy number
                      fCheckOverlaps);  // checking overlaps
   
-      new G4PVPlacement( bgomat ,      // bgo orientation
+      new G4PVPlacement( bgoAmat ,      // bgo orientation
-                     bgo_pos,          // bgo location
+                     bgoA_pos,          // bgo location
-                     bgoLV[i/4],      // its logical volume                         
+                     fLogicBGO[i/4],      // its logical volume                         
-                     "Bgo",         // its name
+                     "BGO",         // its name
+                     worldLV,          // its mother  volume
+                     false,            // no boolean operation
+                     iClover,                // copy number
+                     fCheckOverlaps);  // checking overlaps
+      
+      new G4PVPlacement( bgoSmat ,      // bgo orientation
+                     bgoS_pos,          // bgo location
+                     bgoshieldLV[i/4],      // its logical volume                         
+                     "BGO SHIELD",         // its name
                      worldLV,          // its mother  volume
                      false,            // no boolean operation
                      iClover,                // copy number
@@ -355,8 +404,9 @@ G4VPhysicalVolume* CloverDetectorConstruction::DefineVolumes()
                      fCheckOverlaps);  // checking overlaps
                                        
   caseLV[iClover]->SetVisAttributes(new G4VisAttributes(G4Colour(1.0,1.0,0.0)));
-  bgoLV[iClover]->SetVisAttributes(new G4VisAttributes(G4Colour(0.,1.0,0.0)));
+  fLogicBGO[iClover]->SetVisAttributes(new G4VisAttributes(G4Colour(1.0,0.0,0.0)));
   bgowindowLV[iClover]->SetVisAttributes(new G4VisAttributes(G4Colour(0.0,1.0,1.0)));
+  bgoshieldLV[iClover]->SetVisAttributes(new G4VisAttributes(G4Colour(0.0,1.0,0.0)));
     }
        //Set a i+1%4 for each 
      //
@@ -383,8 +433,11 @@ void CloverDetectorConstruction::ConstructSDandField()
   CloverCrystalSD * crystalSD = new CloverCrystalSD("Crystal", "CrystalHitsCollection", fNumOfCrystal);
   G4SDManager::GetSDMpointer()->AddNewDetector(crystalSD);
 
+  BGOSD * bgoSD = new BGOSD("BGO", "BGOHitsCollection", fNumOfCrystal/4);
+  G4SDManager::GetSDMpointer()->AddNewDetector(bgoSD);
   //Set crystalSD to all logical volumn under the same name of "CrystalLV"
   SetSensitiveDetector("CrystalLV",crystalSD, true);
+  SetSensitiveDetector("BGOLV",bgoSD, true);
 
   // Create global magnetic field messenger.
   // Uniform magnetic field is then created automatically if

CloverDetectorConstruction.hh

@@ -29,6 +29,8 @@ class CloverDetectorConstruction : public G4VUserDetectorConstruction
     G4int   fNumOfCrystal;     // number of Crystal
 
     G4LogicalVolume**  fLogicCrystal;    // pointer to the logical Crystal
+    
+    G4LogicalVolume**  fLogicBGO;
 
     G4Material*        fCrystalMaterial; // pointer to the crystal material
     

CloverEventAction.cc

@@ -5,6 +5,9 @@
 #include "CloverAnalysis.hh"
 #include "Detector_Parameters.hh"
 
+#include "BGOSD.hh"
+#include "BGOHIT.hh"
+
 
 #include "G4RunManager.hh"
 #include "G4Event.hh"
@@ -25,13 +28,19 @@
 CloverEventAction::CloverEventAction()
  : G4UserEventAction(),
    fCrystalHCID(-1),
-   fNDet(Crystal_Num) // number of crystal
+   fBGOHCID(-1),
+   fNBGODet(Crystal_Num/4), // number of bgo
+   fNDet(Crystal_Num) // number of crystal   
 {
   fdEList.clear();
   fdLList.clear();
+  fBEList.clear();
+  fBLList.clear();
   for(int i = 0 ; i < fNDet; i++){
     fdEList.push_back(0.);
     fdLList.push_back(0.);
+    fBEList.push_back(0.);
+    fBLList.push_back(0.);
   }
 
 }
@@ -57,6 +66,19 @@ CloverCrystalHitsCollection* CloverEventAction::GetHitsCollection(G4int hcID, co
   return hitsCollection;
 }    
 
+BGOHitsCollection* CloverEventAction::GetHitsBGOCollection(G4int hcID, const G4Event* event) const
+{
+  auto hitsCollection  = static_cast<BGOHitsCollection*>(event->GetHCofThisEvent()->GetHC(hcID));
+  
+  if ( ! hitsCollection ) {
+    G4ExceptionDescription msg;
+    msg << "Cannot access hitsCollection ID " << hcID; 
+    G4Exception("CloverEventAction::GetHitsCollection() BGOHIT",
+      "MyCode0003", FatalException, msg);
+  }         
+
+  return hitsCollection;
+}    
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 void CloverEventAction::PrintEventStatistics( G4double absoEdep, G4double absoTrackLength) const
@@ -118,6 +140,9 @@ void CloverEventAction::EndOfEventAction(const G4Event* event)
     fCrystalHCID = G4SDManager::GetSDMpointer()->GetCollectionID("CrystalHitsCollection");
   }
 
+  if ( fBGOHCID == -1 ) {
+    fBGOHCID = G4SDManager::GetSDMpointer()->GetCollectionID("BGOHitsCollection");
+  }
   // Get hits collections
   auto crystalHC = GetHitsCollection(fCrystalHCID, event); //this is G4VHitsCollection
 
@@ -148,9 +173,41 @@ void CloverEventAction::EndOfEventAction(const G4Event* event)
   */  
     fdEList[i] = edep * resol * 1000;  // to keV
     fdLList[i] = crystalHit->GetStepLength();
-
+  //  std::cout<< "crystalhere "<< i << ", " << fdEList[i] << std::endl;
   }
 
+  auto BGOHC = GetHitsBGOCollection(fBGOHCID, event); //this is G4VHitsCollection
+
+  for( int i = 0 ; i < fNBGODet ; i++){
+    BGOHit * bgoHit = (*BGOHC)[i];  //this is CloverCrystalHit :: G4VHit
+
+    //add detector resolution of 1%
+
+    G4double edep = bgoHit->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);
+  */  
+    fBEList[i] = edep * resol * 1000;  // to keV
+    fBLList[i] = bgoHit->GetStepLength();
+
+   // std::cout<< "bgohere "<< i << ", " << fBEList[i] << std::endl;
+  }
   
   // get analysis manager
   auto analysisManager = G4AnalysisManager::Instance();
@@ -158,9 +215,9 @@ void CloverEventAction::EndOfEventAction(const G4Event* event)
   // fill ntuple
   analysisManager->FillNtupleIColumn(0, n_trajectories);
 
-  analysisManager->FillNtupleDColumn(3, beamEnergy);
+  analysisManager->FillNtupleDColumn(5, beamEnergy);
-  analysisManager->FillNtupleDColumn(4, beamTheta);
+  analysisManager->FillNtupleDColumn(6, beamTheta);
-  analysisManager->FillNtupleDColumn(5, beamPhi);
+  analysisManager->FillNtupleDColumn(7, beamPhi);
   
   analysisManager->AddNtupleRow();
 

CloverEventAction.hh

@@ -5,6 +5,7 @@
 #include "G4UserEventAction.hh"
 
 #include "CloverCrystalHit.hh"
+#include "BGOHIT.hh"
 
 #include "globals.hh"
 
@@ -28,20 +29,28 @@ public:
   std::vector<G4double>& GetdEList()  {return fdEList;};
   std::vector<G4double>& GetStepLengthList()  {return fdLList;};
 
+  std::vector<G4double>& GetBEList000()  {return fBEList;};
+  std::vector<G4double>& GetBStepLengthList()  {return fBLList;};
 
 private:
   // methods
   CloverCrystalHitsCollection* GetHitsCollection(G4int hcID, const G4Event* event) const;
+  BGOHitsCollection*           GetHitsBGOCollection(G4int hcID, const G4Event* event) const;
   void PrintEventStatistics(G4double absoEdep, G4double absoTrackLength) const;
   
   // data members
   G4int  fCrystalHCID; // Hit collection ID
-
+  G4int  fBGOHCID; // Hit collection ID
+//clover
   G4int fNDet;
+//bgo  
+  G4int fNBGODet;
 
   std::vector<G4double> fdEList;   // dE of each crystal
   std::vector<G4double> fdLList;   // step length of each crystal
 
+  std::vector<G4double> fBEList;   // dE of each crystal
+  std::vector<G4double> fBLList;   // step length of each crystal
   //std::vector<G4double> fLastPost; //last tracking position
   
   

CloverPrimaryGeneratorAction.cc

@@ -43,9 +43,11 @@ CloverPrimaryGeneratorAction::CloverPrimaryGeneratorAction()
  
   //Co60
 
-  fEnergyList = {{1.173,99.973},
+  //fEnergyList = {{1.173,99.973},
-                {1.332,99.986}};
+  //              {1.332,99.986}};
 
+  fEnergyList = {{4.173,99.973},
+                {4.332,99.986}};
   //N16 iso
  /*
   fEnergyList = {{0.120, 100}, //energy, branch

CloverRunAction.cc

@@ -37,6 +37,10 @@ CloverRunAction::CloverRunAction(CloverEventAction * eventAction)
   analysisManager->CreateNtupleDColumn("energy", fEventAction->GetdEList() );
   analysisManager->CreateNtupleDColumn("stepLength", fEventAction->GetStepLengthList() );
 
+  //BGO
+  analysisManager->CreateNtupleDColumn("bgoE", fEventAction->GetBEList000() );
+  analysisManager->CreateNtupleDColumn("bgostep", fEventAction->GetBStepLengthList() );
+
   // event inital angle, energy
   analysisManager->CreateNtupleDColumn("bEnergy");
   analysisManager->CreateNtupleDColumn("theta");

Detector_Parameters.hh

@@ -17,7 +17,7 @@ const double Radius = 200.; //Z POSITION-- radius away from target, in milimeter
                             //
 
 const double Clover_Location[Crystal_Num/4][4] = {{0,131.75,326.96,0.},
-  
+ 
     {1,150.,243.57,0.},
     {2,90.0,257.14,0.},
     {3,90.0,205.71,0.},
@@ -33,7 +33,7 @@ const double Clover_Location[Crystal_Num/4][4] = {{0,131.75,326.96,0.},
     {13,48.25,180.,0.},
     {14,129.,180.,0.},
     {15,48.25,107.35,0.}
-   
+  
     };//,