Dear Experts
I am beginner in Analysis. I am plotting the Energy distribution for the hottest cell/bin and its nearest neighbours. (4 Nearest Neighbours) How can I print the coordinates of this hottest bin with highest energy and coordinates of the 4 nearest neighbours of this cell/bin My Code is as Follows Please suggest in the code
Thanks
#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>
stringstream ss;
void celldepC()
{
// *************************************************************************
// Rootfile and Tree
TFile* f = new TFile("lumi1.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 = 66200 + 0.5;
int nz = 200 + 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);
TRandom3 *rndm = new TRandom3();
// Histograms
// *************************************************************************
// *************************************************************************
// Event Loop
// The first for loop to fill 3D histogram
// int nev_max = 10000;
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];
}
}*/
// Event Loop
// *************************************************************************
// *************************************************************************
map<vector<int>, TH1*> h1max; // Added h1max definition
// *************************************************************************
// Plot Energy spectrum of the max bin
// 27 X 1D histograms around and of the bin with the max energy deposit
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};
TH1* hist = new TH1F(TString(ss.str()), "", 100, 0., 0.01); // Moved the TH1 creation before assigning it to h1max
hist->GetXaxis()->SetTitle("Energy (GeV)");
hist->GetYaxis()->SetTitle("Counts");
h1max[v] = hist; // Assign the TH1 to h1max
}
}
}
int nevts_nohits = 0;
myReader.Restart();
TRandom3 rnd;
rnd.SetSeed(0); // Set the seed for reproducibility
while (myReader.Next()) {
// Identify Max Energy Bin
// *************************************************************************
float e_max = -1;
float x_max = -999;
float y_max = -999;
float z_max = -999;
if (energy.GetSize() == 0)
nevts_nohits++;
//std::cout << "Empty Event: " << nev << std::endl;
for (int i = 0; i < energy.GetSize(); ++i) {
if (energy[i] > e_max) {
e_max = energy[i];
x_max = x[i];
y_max = y[i];
z_max = abs(z[i]) - zmin;
}
}
cout << "Entries: "<< energy.GetSize() <<endl;
nev++;
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;
// Apply Gaussian smearing with 1 cm resolution
//float smeared_energy = rnd.Gaus(energy[i], 0.01);
//cout << ix << endl;
vector<int> v = {ix, iy, iz};
if (h1max.count(v) > 0) {
h1max[v]->Fill(energy[i]); // Use the correct histogram from h1max to fill
}
}
}
cout << "Total number of events: " << nev << ", without hits: " << nevts_nohits << endl;
gStyle->SetOptStat(111111);
// 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("");
ss.str(""); ss << "test_" << ix << "_" << iy << "_" << iz << ".png";
c->SaveAs(TString(ss.str()));
c->Close();
}
}