Incorrect NTuple filling

kitu.utkarsh · November 25, 2023, 6:54am

Hello Experts,
I am trying to fill a Ntuple with total energy deposited in a given detector. I tried modifying example B5/B4c for this.
I gave a condition to sore the energy if its greater than around 0.1 MeV, however still I am getting a large peak at 0 upon plotting the histogram. Manually inspecting the TTree values also reveal there are many 0 entries.
Can someone please share what I am messing here ?

Here are the relevant parts of application

Creating ntuples

B1RunAction::B1RunAction() {
  auto analysisManager = G4AnalysisManager::Instance();
  
  analysisManager->SetVerboseLevel(1);
  analysisManager->SetNtupleMerging(true);

  analysisManager->CreateNtuple("B4", "Edep");
  analysisManager->CreateNtupleDColumn("Eabs");
  analysisManager->FinishNtuple();
}

Opening root files

void B1RunAction::BeginOfRunAction(const G4Run*)
{ 
  auto analysisManager = G4AnalysisManager::Instance();
  
  G4String fileName = "B4";
  analysisManager->OpenFile(fileName);
}

Filling NTuples (User Eventaction class)

void B1EventAction::EndOfEventAction(const G4Event* event)
{   
    // get analysis manager
  auto analysisManager = G4AnalysisManager::Instance();
  
  //Calorimeter Hits
  array<G4int, kDim> totalCalHit = {{  0 }}; 
  array<G4double, kDim> totalCalEdep = {{ 0. }}; 
  
    for (G4int iDet = 0; iDet < kDim; ++iDet) {
	 auto hc = GetHitsCollection(fCalHCID[iDet], event); //kdim =1
	 if ( ! hc ) return;
	 
	 totalCalHit[iDet] = 0;
	 totalCalEdep[iDet] = 0.;

    for (unsigned long i = 0; i < hc->GetSize(); ++i){
      G4double edep = 0.;
      if (iDet == 0) {					// idet = 0 means calorimeter detector
        auto hit = static_cast<B2TrackerHit*>(hc->GetHit(i));
        edep = hit->GetEdep();
      } 
    /*else {                                                           //only used when more than one detector type
        auto hit = static_cast<B2TrackerHit*>(hc->GetHit(i));
        edep = hit->GetEdep();
      }*/
      if ( edep > 0. ) {
        totalCalHit[iDet]++;
        totalCalEdep[iDet] += edep;
      }
      fCalEdep[iDet][i] = edep;
    }
     if (totalCalEdep[iDet] > 100*keV)
     {
	     G4cout << "Fill Energy" << totalCalEdep[iDet] << G4endl;
	     analysisManager->FillNtupleDColumn(0, totalCalEdep[iDet]); }//totoal enrgy in one type of detector
    }
  
    analysisManager->AddNtupleRow();

I think that the condition if (totalCalEdep[iDet] > 100*keV) should only add energies above 100keV, this is also also evident from the G4cout dump. But when I see the ntuple it has entries less than 100keV as well.

The co60 spectrum shows a high peak at 0 as well.

Thanks in advance !

Edit 1: Just went through this post and found a similar kind of issue. I am not sure but my best guess is that perhaps the Ntuples are getting initialised by 0 values
If i just comment analysisManager->FillNtupleDColumn(0, totalCalEdep[iDet]);
The spectrum looks like