added BGOSD and BGOHit, compile fail.

.gitignore

@@ -0,0 +1,3 @@
+*.d
+*.so
+*.swp

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,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......

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),
+   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

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
+   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();

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
   
   

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");