From a1c098e50d07e8346b6c4e006012ce92f316008c Mon Sep 17 00:00:00 2001
From: Peter DeRosa <54421997+MagnusSnowleopard@users.noreply.github.com>
Date: Mon, 17 Oct 2022 18:02:56 -0400
Subject: [PATCH] added BGOSD and BGOHit, compile fail.

---
 .CloverDetectorConstruction.cc.swp | Bin 16384 -> 0 bytes
 .gitignore                         |   3 +
 BGOHIT.cc                          |  81 ++++++++++++++++++++++
 BGOHIT.hh                          |  87 +++++++++++++++++++++++
 BGOSD.cc                           | 106 +++++++++++++++++++++++++++++
 BGOSD.hh                           |  32 +++++++++
 CMakeLists.txt                     |   2 +
 CloverDetectorConstruction.cc      |  16 +++--
 CloverDetectorConstruction.hh      |   2 +
 CloverEventAction.cc               |  58 +++++++++++++++-
 CloverEventAction.hh               |  11 ++-
 CloverRunAction.cc                 |   4 ++
 12 files changed, 395 insertions(+), 7 deletions(-)
 delete mode 100644 .CloverDetectorConstruction.cc.swp
 create mode 100644 .gitignore
 create mode 100644 BGOHIT.cc
 create mode 100644 BGOHIT.hh
 create mode 100644 BGOSD.cc
 create mode 100644 BGOSD.hh

diff --git a/.CloverDetectorConstruction.cc.swp b/.CloverDetectorConstruction.cc.swp
deleted file mode 100644
index 8c3f81ce0e6637874b4278485d83483fd3cc84be..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001

literal 16384
zcmeI2UyK_^9mgjvg_g9TwDb`eu2k5@+4Ws~N$|yyF6U%9;yYi@J~zD(9qZk(y}I>|
zwSO*L11T+tM+B9Ecz~z!P*p8b9{8hHs1Jw-cmSTLiYhAg83a;k)rYp>JNw7>dhfzf
zaW80A`sCu7`OTm2{APCdx6>)LmsjX%WscywpO8nMyJ$T2=FK~AeTI-B^O-Bs|HEbW
zjo~b=Q}fEOzw23sr*yqQVP=pvydGAl*7#M&=DO*r8q*!W{AlZ3VA-a+RByM=H(RP<
z^Bv}H4RjC9$G_C2FtLUn89sNHxa0Y5VE7hy6vIe2BY#N+qyl>?a347_H(L@}56^s#
zKKQ`Ko~2}>R6r^q6_5%@1*8H}0jYpgKq~P6senJci#&;;-4PG(+wtEcQ-5#8*Tp#f
zJ5%`|#Q9g^^gmAJzaQtf<Md+OK+?{gar!6XzfVu)$I0@KR6r^q6_5%@1*8H}0jYpg
zKq?>=kP1izqyis80o^3z)37@)0s#2_pS=Iy!i)Hy!LPs%!42>Y@HManE`pQbv*0NB
z=P^S50A2(yfbRhpTmmKV5cm{$`#wVc1^xnF122Q0f@i@~zyKYv1{z=$ynint?}5LA
zpMW2Oi{LBZli(jmQ6G37JOjQDP6Gly_&Dsr3%~_KumPH29UK7!{PrF~egj?vKLR(v
zQ=kVJFu{Z10r1A%guD)Z2Cjo^fP*{05%BKE2>Ctu9e4#i1AYLW1{<IS>fl~*6uk3M
zLVgWi1}^~%^g$kc0UQH=y9@P!O>hP*fFt1bk3dK8C-5xT1#2J&9s!>R$G~3^V7Uo&
za1zXeGME8xVeW5&?}D#_%U~0H83@0qf!kh`5=vF|oYtVp2@MJZO0}}dgO1Ipp?hre
zG2OOXeaqN(nCH<dJ)x9nfer=>(<LuGe&C{;%Z|fbirNf!*YkC|&YYg#SLf&#(<%5p
zTPp0?SX4!o&J}V-;9q)N6je(3QQy?~^obaS<@BaT*4GCe=5DTWk4{`wT2cze)X<`$
zVRP%cX$7$b`X1#p0-Sh>bB-s8DJjxYMOLW-=t7m!`BC4;GmvdUbhk=p=S%9GlFWmu
zdc!_*Ekn1H4iHjM$;4a^tsy(ea7DRjHbg_wRi&p&NmHY)iGl|kgeW;FlBg6TD5Xar
zN)fd{Ci;ApySCZ1o?w{WDjogWOzAmJ7PSg*m~^SJLy5uX7t;P6`c0LVPwwf{Q9Yq(
zRT@n&Tyb1u;<U;FmXT?}j~#<ipHp9Vd(20{$wY%W@Ds|riHCbiY8-L5V@JYB=9P@&
zr$$!A0kGIvRh4Rt518-n(xL9^_#bmUa*b%^_F8zjOSjwH4hF11_dk<Wsnl54a#&Q-
z%jF4CRRtI4ywNB%8u%-Iit9rr5%SG&f(;({^voGr$Z6#j9i6pwyHe4ZpZsX)*no2R
zyeOoVnPb{QER1u_wH%)gVjEA)XLl*mC9#DYIGBwU-O+o@J%09r<(LgOJg$}$$<1ix
zEsf2<A;MZUCMr6R)uKJJT7;sOV+Df)8($bVx$JW1z;=tv>*H1}9PdMC*;ObrF(BR6
zakS2uDPkWKlx7P(b--;_q=94MaGuAu1~|L7WirRJ{M{nmZu-zs)4d{{Q%;>IY%_bu
za*EjE+40%ag`w_YLrElR>qdZ)hIb7;fR|6G1;w@U>R9tiWf@2P#B#_b=5?l<qiLgN
z=x8Pl#MVq#mzymbPZzmH$qu*71Ndn&_Cxekm-&I~;Qp9e7qQbssdOh?SYD(|ls3_`
zblJhkZ~yDHIXXRgc<a(+HZ~4=`3mjUdCxN9b+v$v<aWEp_5Cx^*TdODQXLcvK3ol9
zug6>6<?d=QXe8Ui6)_d56xhW>+(N+xu@mUF2hV2A6#F=C@J<gqdY{1#&q&mAVz)fM
zgWY^_bdvQNYQ{FyXU28|3(eMtUJfT+Sw2jk)H1iaQLpEVnLR5`56e$>I#<dU39Vpx
zwKqKGCJ$;OXTVX~97^mz1C`4Dd57akT3NGsN4Hn>-UW-<<_ex0@p#AF+uwV-nev#P
z(gGz~xz%^Ed(^pjyg7_}5_#Bi3_CCxJ+l*@R6dKNmPr!s6iqy%GC0zPWBH!a@6V)i
zlN)R}Cpa1%l}Hs?;x-T5jAh&M7Jc8fI(UNI!>U<}_r%aJRpC0v&W*XpM&nvv&*m>;
zgRX_EIn(Gqland8)?Q~r8;_KG<gDRV!v<IanY`2r40BWcZU!CS6}KtwG`Ve=8T)1N
zkmO+!_AP-)*w47-W-Dxp2V<6LrxkRvGkK2hA?3tvEmKRf&rc6DS*BZSqld96+h|e7
zyHB&WXKqr<<YJFn67T<8c*p)F5bytZl@G6<<3Ha5Pl9#u5O@pk@2`RzU>o$om%veQ
z54anMdVV45l9yCKDj*e*3P=T{0#X5~fK)&#AQg}b{4W$ZoKTon9>48ZD(kF=I4pB>
z*?3A4icN|{fh`fDp@SqwL%ZBXbQW<y@y#p@_GT+OezwW{Hgh}N!?&UxX6I(Zy0axm
zk-2@M4!N#qzOCDcm5K0N6zgH_FbYG*nQ;^~jKqZ<-Vr9$@~Fdo`UG>i0&9F8Ky&~j
z7$DNzb1()YszVfnt%FZ53=vq9iU2lZvEqZ7y-TfbWV{sS(q0(8nyz9bks_k+5TF@`
z(Ggfr-qUk4@%wzeJySHIZxYQ~k%p1LXaq1zKL0dpBE*RHsV!oy2yKOdWCy`|6W@u%
SP(6Iep^o?mUyUXzC;tXd<G<Pf

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..9d0f0fd
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,3 @@
+*.d
+*.so
+*.swp
diff --git a/BGOHIT.cc b/BGOHIT.cc
new file mode 100644
index 0000000..3936d1b
--- /dev/null
+++ b/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......
diff --git a/BGOHIT.hh b/BGOHIT.hh
new file mode 100644
index 0000000..7048981
--- /dev/null
+++ b/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
diff --git a/BGOSD.cc b/BGOSD.cc
new file mode 100644
index 0000000..2260fc8
--- /dev/null
+++ b/BGOSD.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......
diff --git a/BGOSD.hh b/BGOSD.hh
new file mode 100644
index 0000000..05d9935
--- /dev/null
+++ b/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
+
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0a66633..2913f2f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -49,6 +49,8 @@ set(CLOVER_SCRIPTS
   vis.mac
   my.in
   analysis.C
+  Analyzer.C
+  Analyzer.h
 )
 
 foreach(_script ${CLOVER_SCRIPTS})
diff --git a/CloverDetectorConstruction.cc b/CloverDetectorConstruction.cc
index 3bc0295..bb7a62a 100644
--- a/CloverDetectorConstruction.cc
+++ b/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),
+   fLogicBGO(NULL), fLogicCrystal(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 [] fLogicCrystal;
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -243,7 +246,7 @@ G4VPhysicalVolume* CloverDetectorConstruction::DefineVolumes()
 
 
     auto caseLV =  new G4LogicalVolume * [fNumOfCrystal/4];
-    auto bgoLV  =  new G4LogicalVolume * [fNumOfCrystal/4];
+//    fLogicBGO  =  new G4LogicalVolume * [fNumOfCrystal/4];
     auto bgowindowLV  =  new G4LogicalVolume * [fNumOfCrystal/4];
 
     fLogicCrystal[i]->SetVisAttributes(crystalVisAtt);
@@ -324,7 +327,7 @@ 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_Al"), "BGO ");
       bgowindowLV[iClover]   = new G4LogicalVolume(bgowindow, G4Material::GetMaterial("G4_Al"), "BGO WINDOW ");
      
       new G4PVPlacement( rotmat ,                // case orientation
@@ -338,7 +341,7 @@ G4VPhysicalVolume* CloverDetectorConstruction::DefineVolumes()
   
       new G4PVPlacement( bgomat ,      // bgo orientation
                      bgo_pos,          // bgo location
-                     bgoLV[i/4],      // its logical volume                         
+                     fLogicBGO[i/4],      // its logical volume                         
                      "Bgo",         // its name
                      worldLV,          // its mother  volume
                      false,            // no boolean operation
@@ -355,7 +358,7 @@ 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(0.,1.0,0.0)));
   bgowindowLV[iClover]->SetVisAttributes(new G4VisAttributes(G4Colour(0.0,1.0,1.0)));
     }
        //Set a i+1%4 for each 
@@ -383,8 +386,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
diff --git a/CloverDetectorConstruction.hh b/CloverDetectorConstruction.hh
index db39781..69cf223 100644
--- a/CloverDetectorConstruction.hh
+++ b/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
     
diff --git a/CloverEventAction.cc b/CloverEventAction.cc
index 6fe240d..053a3d1 100644
--- a/CloverEventAction.cc
+++ b/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
+   fNDet(Crystal_Num), // number of crystal   
+   fBGOHCID(-1),
+   fNBGODet(Crystal_Num/4) // number of bgo
 {
   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
 
@@ -151,6 +176,37 @@ void CloverEventAction::EndOfEventAction(const G4Event* event)
 
   }
 
+  auto BGOHC = GetHitsBGOCollection(fBGOHCID, event); //this is G4VHitsCollection
+
+  for( int i = 0 ; i < fNDet ; 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();
+
+  }
   
   // get analysis manager
   auto analysisManager = G4AnalysisManager::Instance();
diff --git a/CloverEventAction.hh b/CloverEventAction.hh
index 5a0a3d8..214ed3e 100644
--- a/CloverEventAction.hh
+++ b/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
   
   
diff --git a/CloverRunAction.cc b/CloverRunAction.cc
index 7c0291a..ef16531 100644
--- a/CloverRunAction.cc
+++ b/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");