Coordinates of hottest cell and energy distribution

Dear experts below is the part of the code in which i am studying energy depositions in crystal cells. I need to see coordinates of the hottest crystal cell and its energy distribution and also distributions for several near by cells.
here cells are bins Please help in this regard.

Screenshot from 2023-05-19 23-32-46


   `

Please read tips for efficient and successful posting and posting code

___`  map<string,TH3*> h3;

    float xmin = -50;
    float xmax = 50;
    int nx = 100;
    float ymin = -50;
    float ymax = 50;
    int ny = 100;
    float zmin = 66000;
    float zmax = 66300;
    int nz = 300;


    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();
    int nev = 0;
    vector<int> nphoton;
    float xxx = 0;
    float yyy = 0;
    float zzz = 0;
    float emax = 0;
    double binx= 0;
    double biny= 0;
    double binz= 0;
    while (myReader.Next())
    {
        nev++;
        cout << "Event: " << nev << endl;
        nphoton.push_back(0);
        float tot_e = 0;

        for (int i = 0; i < energy.GetSize(); i++)
        {
            nphoton.back()++;
            tot_e += energy[i];
            h3["energy"]->Fill(x[i],y[i],abs(z[i]) - zmin,energy[i]);
            //h3["energy"]->Fill(x[i]+rndm->Gaus(0,1) ,y[i] + rndm->Gaus(0,1) , abs(z[i]) - zmin,energy[i] );
            xxx = max(abs(x[i]), xxx);
            yyy = max(abs(y[i]), yyy);
            zzz = max(abs(z[i]), zzz);
             
   
        }
    }

    cout << "Max-X: " << xxx << endl;
    cout << "Max-Y: " << yyy << endl;
    cout << "Max-Z: " << zzz << endl;
    //cout << "Energy: " << energy << "GeV" <<endl;
    // Fill Histograms
    // *************************************************************************

    // *************************************************************************
    // Draw Histograms
    for (map<string,TH3*>::const_iterator itr = h3.begin(); itr != h3.end(); itr++)
    {
        //gStyle->SetOptStat(0);

        string hname = itr->first;
        TH3* h = (TH3*)itr->second->Clone();
        h->Scale(1.0/float(h->Integral()));

       
        TH2* h2yx = (TH2*)h->Project3D("yx")->Clone();
        TH2* h2zx = (TH2*)h->Project3D("zx")->Clone();
        TH2* h2zy = (TH2*)h->Project3D("zy")->Clone();
        TH1* h1x = (TH2*)h->Project3D("x")->Clone();
        TH1* h1y = (TH2*)h->Project3D("y")->Clone();
        TH1* h1z = (TH2*)h->Project3D("z")->Clone();
        TH1* h1E = (TH2*)h->Project3D("energy")->Clone();
        TH3* h3xyz = (TH3*)h->Project3D("xyz")->Clone();
        
        c[0]->cd();
        h2yx->GetXaxis()->SetTitle("x (mm)");
        h2yx->GetYaxis()->SetTitle("y (mm)");
        gPad->SetLogx(0); gPad->SetLogy(0); gPad->SetLogz(0);
        h2yx->Draw("colz");
        cout << "Bin Max XY:" << h2yx->GetMaximumBin() <<endl;
        cout << "Energy XY" << h2yx->GetMaximum() <<endl;
        
        //ss.str(""); ss << "test2D_yx_" << hname << ".png";
        //c->SaveAs(TString(ss.str()));
        
        c[1]->cd();
        h2zx->GetXaxis()->SetTitle("x (mm)");
        h2zx->GetYaxis()->SetTitle("z (mm)");
        gPad->SetLogx(0); gPad->SetLogy(0); gPad->SetLogz(0);
        h2zx->Draw("colz");
      
        c[2]->cd();
        h2zy->GetXaxis()->SetTitle("y (mm)");
        h2zy->GetYaxis()->SetTitle("z (mm)");
        gPad->SetLogx(0); gPad->SetLogy(0); gPad->SetLogz(0);
        h2zy->Draw("colz");
        
        c[3]->cd();
        h1x->GetXaxis()->SetTitle("x (mm)");
        gPad->SetLogx(0); gPad->SetLogy(1); gPad->SetLogz(0);
        h1x->Draw("colz");
      
        c[4]->cd();
        h1y->GetXaxis()->SetTitle("y (mm)");
        gPad->SetLogx(0); gPad->SetLogy(1); gPad->SetLogz(0);
        h1y->Draw("colz");
       
        c[5]->cd();
        h1z->Draw("colz");
        h1z->GetXaxis()->SetTitle("z (mm)");
        h1z->GetXaxis()->SetRangeUser(0., 200);
        h1z->SetTitle("Energy distribution");
        h1z->GetYaxis()->SetTitle("#Delta E/#Delta X (GeV/mm)");
        gPad->SetLogx(0); gPad->SetLogy(0); gPad->SetLogz(0);
      
       
         c[6]->cd();
        h3xyz->Draw("LEGO2");
        h3xyz->SetTitle("");
        h3xyz->GetXaxis()->SetTitle("X (mm)");
        h3xyz->GetYaxis()->SetTitle("Y (mm)");
        h3xyz->GetZaxis()->SetTitle("Z (mm)");
       

ROOT Version: Not Provided
Platform: Not Provided
Compiler: Not Provided


I think you need to be more specific. What is not working in your code? i.e. what are you trying to do, what are you getting and what exactly do you want to get (and in which form: histogram/s? a value? a list/table of values?..)

Thanks for your response. Actually I have code which is taking Root file and reading Energy X Y and Z positions for Fibers. The coordinates of the hottest (highest energy fiber) is X Y Z(-2.5, -4.5, 34.5) (Already printed in the code).

I would like to have 1D Energy vs Z(mm) distribution for this hottest fiber and neighboring 8 fibers. Please help to add such distributions in the following code.

#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 fiber_test()
{
    // *************************************************************************
    // Rootfile and Tree
    TFile* f = new TFile("lumidet_FBS.edm4hep.root");
    f->cd();

    //TTree* t = (TTree*)f->Get("events");
    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");
    
    
     const int ngraph = 6; 
 TCanvas * c[ngraph]; 
 for (int i =0; i<ngraph; ++i){
 c[i] = new TCanvas(Form("c%d",i),Form("c%d",i),1200,1000);
 c[i]->SetMargin(0.09, 0.1 ,0.1,0.06);
 }
 
   // TCanvas *can = new TCanvas("Energy in cells", "Energy in cells", 800,800);
    //can->Divide(3,2);
 
    // Rootfile and Tree
    // *************************************************************************
    // *************************************************************************
    // Histograms
    map<string,TH3*> h3;

    float xmin = -50;
    float xmax = 50;
    int nx = 100;
    float ymin = -5;
    float ymax = 40;
    int ny = 45;
    float zmin = 66150;
    float zmax = 66350;
    int nz = 200;

    
    h3["energy"] = new TH3F("energy", "Total Energy Deposit;x;y;z", nx, xmin, xmax, ny, ymin, ymax, nz, 0, zmax - zmin);
    // Create a 1D histogram for the hottest bin
 
    
    // TRandom3 *rndm = new TRandom3();
    int nev = 0;
    vector<int> nphoton;
    float xxx = 0;
    float yyy = 0;
    float zzz = 0;
    float emax = 0;
    double hottest_x = 0.0;
    double hottest_y = 0.0;
    double hottest_z = 0.0;
   float  hottest_energy = 0.0;
  
     
    while (myReader.Next())
    {
        nev++;
        cout << "Event: " << nev << endl;
        nphoton.push_back(0);
        float tot_e = 0;

        for (int i = 0; i < energy.GetSize(); i++)
        {
            nphoton.back()++;
            tot_e += energy[i];
            float z_shifted = abs(z[i]) - zmin;
            h3["energy"]->Fill(x[i],y[i],abs(z[i]) - zmin,energy[i]);
            //h3["energy"]->Fill(x[i]+rndm->Gaus(0,1) ,y[i] + rndm->Gaus(0,1) , abs(z[i]) - zmin,energy[i] );
            xxx = max(abs(x[i]), xxx);
            yyy = max(abs(y[i]), yyy);
            zzz = max(abs(z[i]), zzz);
             
    	    //cout << "- Pos [" << i << "]:" << x[i] << " " << y[i] << " " << z[i] << " :" << z_shifted << endl;
           
        }
    }

    cout << "Max-X: " << xxx << endl;
    cout << "Max-Y: " << yyy << endl;
    cout << "Max-Z: " << zzz << endl;

    // Fill Histograms
    // *************************************************************************

    // *************************************************************************
    // Draw Histograms
    for (map<string,TH3*>::const_iterator itr = h3.begin(); itr != h3.end(); itr++)
    {
        //gStyle->SetOptStat(0);

        string hname = itr->first;
        TH3* h = (TH3*)itr->second->Clone();
        h->Scale(1.0/float(h->Integral()));

       
        TH2* h2yx = (TH2*)h->Project3D("yx")->Clone();
        TH2* h2zx = (TH2*)h->Project3D("zx")->Clone();
        TH2* h2zy = (TH2*)h->Project3D("zy")->Clone();
        TH1* h1x = (TH2*)h->Project3D("x")->Clone();
        TH1* h1y = (TH2*)h->Project3D("y")->Clone();
        TH1* h1z = (TH2*)h->Project3D("z")->Clone();
        
        
        c[0]->cd();
        h2yx->GetXaxis()->SetTitle("x (mm)");
        h2yx->GetYaxis()->SetTitle("y (mm)");
        gPad->SetLogx(0); gPad->SetLogy(0); gPad->SetLogz(0);
        h2yx->Draw("colz");
                //ss.str(""); ss << "test2D_yx_" << hname << ".png";
        //c->SaveAs(TString(ss.str()));
        
        c[1]->cd();
        h2zx->GetXaxis()->SetTitle("x (mm)");
        h2zx->GetYaxis()->SetTitle("z (mm)");
        gPad->SetLogx(0); gPad->SetLogy(0); gPad->SetLogz(0);
        h2zx->Draw("colz");
      
        c[2]->cd();
        h2zy->GetXaxis()->SetTitle("y (mm)");
        h2zy->GetYaxis()->SetTitle("z (mm)");
        gPad->SetLogx(0); gPad->SetLogy(0); gPad->SetLogz(0);
        h2zy->Draw("colz");
        
        c[3]->cd();
        h1x->GetXaxis()->SetTitle("x (mm)");
        gPad->SetLogx(0); gPad->SetLogy(1); gPad->SetLogz(0);
        h1x->Draw("colz");
      
        c[4]->cd();
        h1y->GetXaxis()->SetTitle("y (mm)");
        gPad->SetLogx(0); gPad->SetLogy(1); gPad->SetLogz(0);
        h1y->Draw("colz");
       
        c[5]->cd();
        h1z->Draw("colz");
        h1z->GetXaxis()->SetTitle("z (mm)");
        h1z->GetXaxis()->SetRangeUser(0., 200);
        h1z->SetTitle("Energy distribution");
        h1z->GetYaxis()->SetTitle("E(GeV)");
        gPad->SetLogx(0); gPad->SetLogy(0); gPad->SetLogz(0);
        
        
     
    }
    // Draw Histograms
    // *************************************************************************
  // Find the hottest bin with the highest energy
    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++) {
                float energy = h3["energy"]->GetBinContent(ix, iy, iz);
                if (energy > hottest_energy) {
                    hottest_x = h3["energy"]->GetXaxis()->GetBinCenter(ix);
                    hottest_y = h3["energy"]->GetYaxis()->GetBinCenter(iy);
                    hottest_z = h3["energy"]->GetZaxis()->GetBinCenter(iz);
                    hottest_energy = energy;
                }
            }
        }
    }

    cout << "Hottest Fiber coordinates: X Y Z(" << hottest_x << ", " << hottest_y << ", " << hottest_z << ")" << endl;
    //cout << "Hottest cell total energy: " << hottest_energy << endl


}

Once you’ve identified the “hottest cell”, you need to identify its fiber and its 8 neighboring fibers and then loop over the tree entries again (and fill the 1+8 new 1D histograms with their energy distributions).
Note: you will need something like this: myReader.Restart(); while (myReader.Next()) { ... }

@Wile_E_Coyote Thanks a lot for your response. I am not very much Expert it will be great if you kindly add code for that It will be very helpful for me

This requires detailed knowledge of your detector, so you need to talk to your colleagues and/or supervisor about it.

@Wile_E_Coyote Can You help in the Following I need to fill hced with Maximum bin in X Y and Z from h3[energy]. So hced will contain the highest energy. How can I do that?

map<string,TH3*> h3;

float xmin = -50;
float xmax = 50;
int nx = 100;
float ymin = -50;
float ymax = 50;
int ny = 100;
float zmin = 66000;
float zmax = 66300;
int nz = 300;


TH1F* hced = new TH1F("hced", "Hottest cell total E. distr.", 1000, 0.0, 1);

h3["energy"] = new TH3F("energy", "Total Energy Deposit;x;y;z", nx, xmin, xmax, ny, ymin, ymax, nz, 0, zmax - zmin);
// Create a 1D histogram for the hottest bin


//TRandom3* rndm = new TRandom3(); // outside the loop


int nev = 0;
vector<int> nphoton;
float xxx = 0;
float yyy = 0;
float zzz = 0;
float emax = 0;
double hottest_x = 0.0;
double hottest_y = 0.0;
double hottest_z = 0.0;
float hottest_energy = 0.0;
float hottest_cell_energy = 0.0;
double maxBinContent = 0.0;
int maxBinX = 0, maxBinY = 0, maxBinZ = 0;
while (myReader.Next())
{
    nev++;
    cout << "Event: " << nev << endl;
    nphoton.push_back(0);
    float tot_e = 0;
    
    for (int i = 0; i < energy.GetSize(); i++)
    {
        nphoton.back()++;
        tot_e += energy[i];
        float z_shifted = abs(z[i]) - zmin;
        h3["energy"]->Fill(x[i],y[i],abs(z[i]) - zmin,energy[i]);
        //h3["energy"]->Fill(x[i]+rndm->Gaus(0,10) ,y[i] + rndm->Gaus(0,10) , abs(z[i]) - zmin,energy[i] );
        xxx = max(abs(x[i]), xxx);
        yyy = max(abs(y[i]), yyy);
        zzz = max(abs(z[i]), zzz);
         
	   
        // Find maximum bin


       hced->Fill(energy[i]);
    	    
    
    }
}

cout << "Max-X: " << xxx << endl;
cout << "Max-Y: " << yyy << endl;
cout << "Max-Z: " << zzz << endl;


// *************************************************************************

// *************************************************************************
// Draw Histograms
for (map<string,TH3*>::const_iterator itr = h3.begin(); itr != h3.end(); itr++)
{
    //gStyle->SetOptStat(0);

    string hname = itr->first;
    TH3* h = (TH3*)itr->second->Clone();
    h->Scale(1.0/float(h->Integral()));

   
    TH2* h2yx = (TH2*)h->Project3D("yx")->Clone();
   
    
    c[0]->cd();
    h2yx->GetXaxis()->SetTitle("x (mm)");
    h2yx->GetYaxis()->SetTitle("y (mm)");
    h2yx->GetZaxis()->SetTitle("Energy (GeV)");
    gPad->SetLogx(0); gPad->SetLogy(0); gPad->SetLogz(0);
    h2yx->Draw("colz");
    
    c[1]->cd();
    hced->GetXaxis()->SetTitle("Energy (GeV)");
    hced->GetYaxis()->SetTitle("Counts");
    hced->GetXaxis()->SetRangeUser(maxBinX - 1, maxBinX + 1);
    hced->GetYaxis()->SetRangeUser(maxBinY - 1, maxBinY + 1);
    hced->GetZaxis()->SetRangeUser(maxBinZ - 1, maxBinZ + 1);
    hced->SetLineColor(kRed);
    hced->Draw();

Your code does not seems complete.

@couet here is the complete code What I dont Know is how to fill hced histogram using the maxbin of h3[energy] histogram in X Y and Z. So this hced will be the histogram with maximum bin content which is energy here:


#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 test()
{
    // *************************************************************************
    // Rootfile and Tree
    TFile* f = new TFile("lumidet_PbWO4_1K.edm4hep.root");
    f->cd();

    //TTree* t = (TTree*)f->Get("events");
    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");
    
    
     const int ngraph = 2; 
 TCanvas * c[ngraph]; 
 for (int i =0; i<ngraph; ++i){
 c[i] = new TCanvas(Form("c%d",i),Form("c%d",i),1200,1000);
 c[i]->SetMargin(0.09, 0.1 ,0.1,0.06);
 }
 
   // TCanvas *can = new TCanvas("Energy in cells", "Energy in cells", 800,800);
  
    // Rootfile and Tree
    // *************************************************************************
    // *************************************************************************
    // Histograms
    map<string,TH3*> h3;

    float xmin = -50;
    float xmax = 50;
    int nx = 100;
    float ymin = -50;
    float ymax = 50;
    int ny = 100;
    float zmin = 66000;
    float zmax = 66300;
    int nz = 300;


    TH1F* hced = new TH1F("hced", "Hottest cell total E. distr.", 1000, 0.0, 1);
    
    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(); // outside the loop
  
    
    int nev = 0;
    vector<int> nphoton;
    float xxx = 0;
    float yyy = 0;
    float zzz = 0;
    float emax = 0;
    double hottest_x = 0.0;
    double hottest_y = 0.0;
    double hottest_z = 0.0;
    float hottest_energy = 0.0;
    float hottest_cell_energy = 0.0;
    double maxBinContent = 0.0;
    int maxBinX = 0, maxBinY = 0, maxBinZ = 0;
    double maxValue = -1.0;
    while (myReader.Next())
    {
        nev++;
        cout << "Event: " << nev << endl;
        nphoton.push_back(0);
        float tot_e = 0;
        
        for (int i = 0; i < energy.GetSize(); i++)
        {
            nphoton.back()++;
            tot_e += energy[i];
            float z_shifted = abs(z[i]) - zmin;
            h3["energy"]->Fill(x[i],y[i],abs(z[i]) - zmin,energy[i]);
            //h3["energy"]->Fill(x[i]+rndm->Gaus(0,10) ,y[i] + rndm->Gaus(0,10) , abs(z[i]) - zmin,energy[i] );
            xxx = max(abs(x[i]), xxx);
            yyy = max(abs(y[i]), yyy);
            zzz = max(abs(z[i]), zzz);
             
    	   
            // Find maximum bin
           
Double_t maxValue = -1.0;

for (int binX = 1; binX <= h3["energy"]->GetNbinsX(); binX++) {
    for (int binY = 1; binY <= h3["energy"]->GetNbinsY(); binY++) {
        for (int binZ = 1; binZ <= h3["energy"]->GetNbinsZ(); binZ++) {
            double binContent = h3["energy"]->GetBinContent(binX, binY, binZ);
            if (binContent > maxValue) {
                maxValue = binContent;
                maxBinX = binX;
                maxBinY = binY;
                maxBinZ = binZ;
         
         
         
            }
        }
    }
}
          
           hced->Fill(energy[i]);
   	    	    
   	    
        }
    }

    cout << "Max-X: " << xxx << endl;
    cout << "Max-Y: " << yyy << endl;
    cout << "Max-Z: " << zzz << endl;
   
  
    // *************************************************************************

    // *************************************************************************
    // Draw Histograms
    for (map<string,TH3*>::const_iterator itr = h3.begin(); itr != h3.end(); itr++)
    {
        //gStyle->SetOptStat(0);

        string hname = itr->first;
        TH3* h = (TH3*)itr->second->Clone();
        h->Scale(1.0/float(h->Integral()));

       
        TH2* h2yx = (TH2*)h->Project3D("yx")->Clone();
       
        
        c[0]->cd();
        h2yx->GetXaxis()->SetTitle("x (mm)");
        h2yx->GetYaxis()->SetTitle("y (mm)");
        h2yx->GetZaxis()->SetTitle("Energy (GeV)");
        gPad->SetLogx(0); gPad->SetLogy(0); gPad->SetLogz(0);
        h2yx->Draw("colz");
        
        c[1]->cd();
        hced->GetXaxis()->SetTitle("Energy (GeV)");
        hced->GetYaxis()->SetTitle("Counts");
        hced->GetXaxis()->SetRangeUser(maxBinX - 1, maxBinX + 1);
        hced->GetYaxis()->SetRangeUser(maxBinY - 1, maxBinY + 1);
        hced->GetZaxis()->SetRangeUser(maxBinZ - 1, maxBinZ + 1);
        hced->SetLineColor(kRed);
        hced->Draw();

        
  
     
    }
   
    // *************************************************************************
  
  //  cout << "Hottest cell coordinates: X Y Z(" << hottest_x << ", " << hottest_y << ", " << hottest_z << ")" << endl;
}    

lumidet_PbWO4_1K.edm4hep.root is missing

Sorry @couet Please download the file from the following link as it is larger than 3MB

for each event this loop takes ages…:

         for (int binX = 1; binX <= h3["energy"]->GetNbinsX(); binX++) {
            for (int binY = 1; binY <= h3["energy"]->GetNbinsY(); binY++) {
               for (int binZ = 1; binZ <= h3["energy"]->GetNbinsZ(); binZ++) {
                  double binContent = h3["energy"]->GetBinContent(binX, binY, binZ);
                  if (binContent > maxValue) {
                     maxValue = binContent;
                     maxBinX = binX;
                     maxBinY = binY;
                     maxBinZ = binZ;
                  }
               }
            }
         }

Ok Please suggest an alternative will be very greatful to fill for the hottest cell/bin then I will extend it to neighbours

I am not sure I have an alternative … it is running right now … just super long.

Ok Please suggest where to fill hced histogram

Or If you can send the complete code If You have made some changes

What sense does it make to search for the “hottest cell/bin” after every single “Fill” call?

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.