Dear Experts
I am studying the Energy deposition in crystals with my bin/cell size = 1mm. But now want to study the deposition in circular fiber size of 1mm. How can I add cut because otherwise there is a chance that I can also deposit around the fibers as for crystal case bins are square. 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>
stringstream ss;
void crystal_dep()
{
// *************************************************************************
// Rootfile and Tree
TFile* f = new TFile("lumi_crystal_segmented.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");
// Rootfile and Tree
// *************************************************************************
// *************************************************************************
// Histograms
map<string,TH3*> h3;
float xmin = -100;
float xmax = 100;
int nx = 200;
float ymin = -100;
float ymax = 100;
int ny = 200;
float zmin = 66000;
float zmax = 66500;
int nz = 500;
h3["energy"] = new TH3F("energy", "Total Energy Deposit;x;y;z", nx, xmin, xmax, ny, ymin, ymax, nz, 0, zmax - zmin);
/*
map<string,vector<TH1*>> h1;
for (int i = 0; i < h3["hits"]->GetNbinsX()*h3["hits"]->GetNbinsY()*h3["hits"]->GetNbinsZ(); i++)
{
ss.str(""); ss << "h1_cell_energy_" << i;
h1["energy"].push_back(new TH1F(TString(ss.str()), "", 100, 0, 0.001));
ss.str(""); ss << "Cell " << i << ";E_{#gamma} (GeV); Entries";
h1["energy"][i]->SetTitle(TString(ss.str()));
}
*/
// Histograms
// *************************************************************************
// *************************************************************************
// Fill Histograms
int nev = 0;
vector<int> nphoton;
float xxx = 0;
float yyy = 0;
float zzz = 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]);
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;
// 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()));
TCanvas* c = new TCanvas("c", "", 1000, 1000);
/*
c->cd();
h->Draw("BOX2Z");
h->Draw();
gPad->SetLogx(0);
gPad->SetLogy(0);
gPad->SetLogz(0);
ss.str(""); ss << "test_" << hname << ".png";
c->SaveAs(TString(ss.str()));
*/
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->cd();
h2yx->Draw("colz");
gPad->SetLogx(0); gPad->SetLogy(0); gPad->SetLogz(0);
ss.str(""); ss << "test2D_yx_" << hname << ".png";
c->SaveAs(TString(ss.str()));
h2zx->Draw("colz");
gPad->SetLogx(0); gPad->SetLogy(0); gPad->SetLogz(0);
ss.str(""); ss << "test2D_zx_" << hname << ".png";
c->SaveAs(TString(ss.str()));
h2zy->Draw("colz");
gPad->SetLogx(0); gPad->SetLogy(0); gPad->SetLogz(0);
ss.str(""); ss << "test2D_zy_" << hname << ".png";
c->SaveAs(TString(ss.str()));
h1x->Draw("colz");
gPad->SetLogx(0); gPad->SetLogy(1); gPad->SetLogz(0);
ss.str(""); ss << "test1D_x_" << hname << ".png";
c->SaveAs(TString(ss.str()));
h1y->Draw("colz");
gPad->SetLogx(0); gPad->SetLogy(1); gPad->SetLogz(0);
ss.str(""); ss << "test1D_y_" << hname << ".png";
c->SaveAs(TString(ss.str()));
h1z->Draw("colz");
gPad->SetLogx(0); gPad->SetLogy(1); gPad->SetLogz(0);
ss.str(""); ss << "test1D_z_" << hname << ".png";
c->SaveAs(TString(ss.str()));
c->Close();
}
// Draw Histograms
// *************************************************************************
}
Please read tips for efficient and successful posting and posting code
ROOT Version: Not Provided
Platform: Not Provided
Compiler: Not Provided