Dear Root Experts
I am trying to fill my all histograms for 1000 Entries but I can see very less number of Entries in each histogram. Could You Please check what could be the Problem.
My code is as follows:
#include <iostream>
#include <iomanip>
#include <sstream>
#include <fstream>
#include <string>
#include <vector>
#include <map>
#include <stdlib.h>
#include "TROOT.h"
#include "TTree.h"
#include "TChain.h"
#include "TFile.h"
#include "TString.h"
#include "TSystem.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TColor.h"
#include "Riostream.h"
#include "TDatime.h"
#include "TMath.h"
#include "TStyle.h"
#include "TLorentzVector.h"
#include <TTreeReader.h>
#include <TTreeReaderValue.h>
#include <TTreeReaderArray.h>
using namespace std;
stringstream ss;
void hotcell()
{
// *************************************************************************
// Rootfile and Tree
TFile* f = new TFile("lumi_1K.edm4hep.root");
f->cd();
TTreeReader myReader("events", f); // name of tree and file
TTreeReaderArray<Float_t> energy(myReader, "LumiDirectPCALHits.energy");
TTreeReaderArray<Float_t> x(myReader, "LumiDirectPCALHits.position.x");
TTreeReaderArray<Float_t> y(myReader, "LumiDirectPCALHits.position.y");
TTreeReaderArray<Float_t> z(myReader, "LumiDirectPCALHits.position.z");
// Rootfile and Tree
// *************************************************************************
// *************************************************************************
// Histograms
map<string, TH1*> h1;
map<string, TH2*> h2;
map<string, TH3*> h3;
float xmin = -50 - 0.5;
float xmax = 50 + 0.5;
int nx = 100 + 1;
float ymin = -50 - 0.5;
float ymax = 50 + 0.5;
int ny = 100 + 1;
float zmin = 66000 - 0.5;
float zmax = 66300 + 0.5;
int nz = 300 + 1;
// 3D histogram x, y, z. total energy deposit as T axis
h3["energy"] = new TH3F("energy", "Total Energy Deposit;x;y;z", nx, xmin, xmax, ny, ymin, ymax, nz, 0, zmax - zmin);
// Histograms
// *************************************************************************
// *************************************************************************
// Event Loop
// The first for loop to fill 3D histogram
int nev_max = 1000; // Set the maximum number of events to process to 1000
int nev = 0;
while (myReader.Next())
{
if (nev >= nev_max) break;
nev++;
if (nev % 100 == 0) cout << "Event: " << nev << endl;
float tot_e = 0;
for (int i = 0; i < energy.GetSize(); i++)
{
tot_e += energy[i];
float z_shifted = abs(z[i]) - zmin;
h3["energy"]->Fill(x[i], y[i], z_shifted, energy[i]);
}
}
// Event Loop
// *************************************************************************
// *************************************************************************
// Identify Max Energy Bin
float e_max = 0;
float x_max = -999;
float y_max = -999;
float z_max = -999;
for (int ix = 1; ix <= h3["energy"]->GetNbinsX(); ix++)
{
for (int iy = 1; iy <= h3["energy"]->GetNbinsY(); iy++)
{
for (int iz = 1; iz <= h3["energy"]->GetNbinsZ(); iz++)
{
if (h3["energy"]->GetBinContent(ix, iy, iz) > e_max)
{
e_max = h3["energy"]->GetBinContent(ix, iy, iz);
x_max = h3["energy"]->GetXaxis()->GetBinCenter(ix);
y_max = h3["energy"]->GetYaxis()->GetBinCenter(iy);
z_max = h3["energy"]->GetZaxis()->GetBinCenter(iz);
}
}
}
}
cout << "max-E : " << e_max << endl;
cout << "max-X : " << x_max << endl;
cout << "max-Y : " << y_max << endl;
cout << "max-Z : " << z_max << endl;
// Identify Max Energy Bin
// *************************************************************************
// *************************************************************************
// Plot Energy spectrum of the max bin
// 27 X 1D histograms around and of the bin with the max energy deposit
map<vector<int>, TH1*> h1max;
for (int ix = -1; ix <= 1; ix++)
{
for (int iy = -1; iy <= 1; iy++)
{
for (int iz = -1; iz <= 1; iz++)
{
ss.str("");
ss << "hist_" << ix << "_" << iy << "_" << iz;
vector<int> v = {ix, iy, iz};
h1max[v] = new TH1F(TString(ss.str()), "", 100, 0, 0.1);
}
}
}
// *************************************************************************
// Loop again to fill h1max histograms
myReader.Restart();
nev = 0; // Reset the event counter to zero
int requiredEntries = 1000;
while (true)
{
if (nev >= requiredEntries) break;
if (!myReader.Next())
{
myReader.Restart();
continue;
}
nev++;
if (nev % 100 == 0) cout << "Event: " << nev << endl;
for (int i = 0; i < energy.GetSize(); i++)
{
int ix = x[i] - x_max;
int iy = y[i] - y_max;
float z_shifted = abs(z[i]) - zmin;
int iz = z_shifted - z_max;
if (abs(ix) > 1 || abs(iy) > 1 || abs(iz) > 1)
continue;
vector<int> v = {ix, iy, iz};
h1max[v]->Fill(energy[i]);
}
}
// Loop again
// *************************************************************************
for (map<vector<int>, TH1*>::const_iterator itr = h1max.begin(); itr != h1max.end(); itr++)
{
int ix = itr->first[0];
int iy = itr->first[1];
int iz = itr->first[2];
TCanvas* c = new TCanvas("c", "", 1000, 1000);
c->cd();
itr->second->Draw("E");
ss.str("");
ss << "test_" << ix << "_" << iy << "_" << iz << ".png";
c->SaveAs(TString(ss.str()));
c->Close();
}
}