diff --git a/ecal/ecal_config.yml b/ecal/ecal_config.yml
index 69ffb803f79caa6f8b0b1352328e9444d577dd7b..a4bb8ff633ac519b6ae9c064e6250e104939c521 100644
--- a/ecal/ecal_config.yml
+++ b/ecal/ecal_config.yml
@@ -3,7 +3,7 @@ ecal_1_emcal_electrons:
tags:
- silicon
needs: ["get_data"]
- timeout: 12 hours 30 minutes
+ timeout: 24 hours
artifacts:
expire_in: 20 weeks
paths:
diff --git a/ecal/emcal_electrons.sh b/ecal/emcal_electrons.sh
index 3173c35d8e1a779ad80be9eb0aaec0d11a4df289..29ec061704e40659f1436f5d54c86a3cee993108 100644
--- a/ecal/emcal_electrons.sh
+++ b/ecal/emcal_electrons.sh
@@ -8,6 +8,14 @@ if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
export JUGGLER_N_EVENTS=100
fi
+if [[ ! -n "${E_start}" ]] ; then
+ export E_start=0.0
+fi
+
+if [[ ! -n "${E_end}" ]] ; then
+ export E_end=30.0
+fi
+
export JUGGLER_FILE_NAME_TAG="emcal_uniform_electrons"
export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
@@ -28,7 +36,8 @@ popd
# generate the input events
# note datasets is now only used to develop datasets.
#git clone https://eicweb.phy.anl.gov/EIC/datasets.git datasets
-root -b -q "ecal/scripts/emcal_electrons.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+root -b -q "ecal/scripts/emcal_electrons.cxx(${JUGGLER_N_EVENTS}, ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+root -b -q "ecal/scripts/emcal_electrons_reader.cxx(${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
pushd ${JUGGLER_DETECTOR}
ls -l
@@ -46,7 +55,12 @@ popd
pwd
mkdir -p results
-root -b -q "ecal/scripts/makeplot.C(\"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\")"
+root -b -q "ecal/scripts/makeplot.C(${E_start}, ${E_end}, \"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\", \"results/rec_${JUGGLER_FILE_NAME_TAG}.txt\")"
+root -b -q "ecal/scripts/makeplot_input.C(\"${JUGGLER_DETECTOR}/${JUGGLER_SIM_FILE}\", \"results/sim_${JUGGLER_FILE_NAME_TAG}.txt\")"
+
+paste results/sim_${JUGGLER_FILE_NAME_TAG}.txt results/rec_${JUGGLER_FILE_NAME_TAG}.txt > results/eng_${JUGGLER_FILE_NAME_TAG}.txt
+root -b -q "ecal/scripts/read_eng.C(\"results/eng_${JUGGLER_FILE_NAME_TAG}.root\", \"results/eng_${JUGGLER_FILE_NAME_TAG}.txt\")"
+root -b -q "ecal/scripts/cal_eng_res.C(\"results/eng_${JUGGLER_FILE_NAME_TAG}.root\")"
#mkdir -p sim_output
#cp "${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}" sim_output/.
diff --git a/ecal/scripts/cal_eng_res.C b/ecal/scripts/cal_eng_res.C
new file mode 100644
index 0000000000000000000000000000000000000000..3aa4220a65877bd1932ad67371550e50fb386136
--- /dev/null
+++ b/ecal/scripts/cal_eng_res.C
@@ -0,0 +1,84 @@
+int cal_eng_res(const char* input_fname = "sim_output/read_eng_output.root")
+{
+ // Setting Figures
+ gROOT->SetStyle("Plain");
+ gStyle->SetLineWidth(3);
+ gStyle->SetOptStat("nem");
+ gStyle->SetOptFit(1);
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetPadGridX(1);
+ gStyle->SetPadGridY(1);
+ gStyle->SetPadLeftMargin(0.14);
+
+ // Read ROOT File
+ TFile *f = new TFile(input_fname,"r");
+ TTree *tEngRes = (TTree *)f->Get("tEngRes");
+ // Variables from ROOT tree
+ Int_t ievent;
+ Double_t tot_clust_eng; // total cluster energy [GeV]
+ Int_t ncluster; // number of clusters per event
+ Double_t mc_eng; // thrown energy [GeV]
+ Double_t sim_tru_eng; // total truth simulated energy [GeV]
+ Double_t sim_eng; // total simulated energy [GeV]
+ // Read Branches
+ tEngRes->SetBranchAddress("ievent", &ievent);
+ tEngRes->SetBranchAddress("tot_clust_eng", &tot_clust_eng);
+ tEngRes->SetBranchAddress("ncluster", &ncluster);
+ tEngRes->SetBranchAddress("mc_eng", &mc_eng);
+ tEngRes->SetBranchAddress("sim_tru_eng", &sim_tru_eng);
+ tEngRes->SetBranchAddress("sim_eng", &sim_eng);
+
+ // Setting for Canvas
+ TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
+ TCanvas *c2 = new TCanvas("c2", "c2", 700, 500);
+
+ // Declare Histrograms
+ TH1D *h1 = new TH1D("hEnergyRes","Energy Resolution", 100,-0.3,0.3);
+ TH1D *h2 = new TH1D("hEnergyResCUT","Energy Resolution with CUT", 100,-0.3,0.3);
+
+ // Variables that used in calculations
+ Double_t diff_energy;
+ Double_t eng_res;
+
+ Int_t nentries = tEngRes->GetEntries();
+
+ for(Int_t ievent = 0; ievent < nentries; ievent++)
+ {
+ tEngRes->GetEntry(ievent);
+
+ if(ncluster == 1)
+ {
+ diff_energy = tot_clust_eng - mc_eng;
+ eng_res = (diff_energy / mc_eng);
+ cout << tot_clust_eng << "\t" << mc_eng << "\t" << eng_res << endl;
+ h1->Fill(eng_res, 1.0);
+
+ if(tot_clust_eng > 0.5)
+ h2->Fill(eng_res, 1.0);
+ }
+ }
+
+ // Drawing and Saving Figures
+ c1->cd();
+ h1->SetLineColor(kBlue);
+ h1->SetLineWidth(2);
+ h1->GetXaxis()->SetTitle("#DeltaE/E");
+ h1->GetYaxis()->SetTitle("events");
+ h1->GetYaxis()->SetTitleOffset(1.4);
+ h1->Fit("gaus");
+ h1->Draw();
+ c1->SaveAs("results/hEngRes.png");
+
+ c2->cd();
+ h2->SetLineColor(kBlue);
+ h2->SetLineWidth(2);
+ h2->GetXaxis()->SetTitle("#DeltaE/E");
+ h2->GetYaxis()->SetTitle("events");
+ h2->GetYaxis()->SetTitleOffset(1.4);
+ h2->Fit("gaus");
+ h2->Draw();
+ c2->SaveAs("results/hEngRes_CUT.png");
+
+ return 0;
+}
diff --git a/ecal/scripts/emcal_electrons.cxx b/ecal/scripts/emcal_electrons.cxx
index 02ac3fe0ec441971a5a48add99858e4f45b2a90f..d17f60925eec621fdab4bcf129c76f0f8eb5f2ca 100644
--- a/ecal/scripts/emcal_electrons.cxx
+++ b/ecal/scripts/emcal_electrons.cxx
@@ -19,7 +19,7 @@
using namespace HepMC3;
-void emcal_electrons(int n_events = 1e3, const char* out_fname = "./data/emcal_electron_0GeVto30GeV_100kEvt.hepmc")
+void emcal_electrons(int n_events = 1e3, double e_start = 1.0, double e_end = 1.0, const char* out_fname = "./data/emcal_electron_0GeVto30GeV_100kEvt.hepmc")
{
WriterAscii hepmc_output(out_fname);
int events_parsed = 0;
@@ -40,7 +40,7 @@ void emcal_electrons(int n_events = 1e3, const char* out_fname = "./data/emcal_e
FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
// Define momentum
- Double_t p = r1->Uniform(0.0, 30.0);
+ Double_t p = r1->Uniform(e_start, e_end);
Double_t px;
Double_t py;
Double_t pz;
diff --git a/ecal/scripts/emcal_electrons_reader.cxx b/ecal/scripts/emcal_electrons_reader.cxx
index 5537e3b8783be44e464f714157bb1621a5835c73..f916d23bc0c1fcd601ba1660a016f346b6d544e1 100644
--- a/ecal/scripts/emcal_electrons_reader.cxx
+++ b/ecal/scripts/emcal_electrons_reader.cxx
@@ -14,7 +14,7 @@
using namespace HepMC3;
-void emcal_electrons_reader(){
+void emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const char* in_fname = "./data/emcal_electron_0GeVto30GeV_100kEvt.hepmc"){
// Setting for graphs
gROOT->SetStyle("Plain");
@@ -27,18 +27,18 @@ void emcal_electrons_reader(){
gStyle->SetPadLeftMargin(0.14);
gStyle->SetPadRightMargin(0.14);
- ReaderAscii hepmc_input("./data/emcal_electron_0GeVto30GeV_100kEvt.hepmc");
+ ReaderAscii hepmc_input(in_fname);
int events_parsed = 0;
GenEvent evt(Units::GEV, Units::MM);
// Histograms
- TH1F* h_electrons_energy = new TH1F("h_electron_energy", "electron energy;E [GeV];Events", 100,-0.5,30.5);
+ TH1F* h_electrons_energy = new TH1F("h_electron_energy", "electron energy;E [GeV];Events", 100,e_start-0.5,e_end+0.5);
TH1F* h_electrons_eta = new TH1F("h_electron_eta", "electron #eta;#eta;Events", 100,-10.0,10.0);
TH1F* h_electrons_theta = new TH1F("h_electron_theta", "electron #theta;#theta [degree];Events", 100,-0.5,180.5);
TH1F* h_electrons_phi = new TH1F("h_electron_phi", "electron #phi;#phi [degree];Events", 100,-180.5,180.5);
- TH2F* h_electrons_pzpt = new TH2F("h_electrons_pzpt", "electron pt vs pz;pt [GeV];pz [GeV]", 100,-0.5,30.5,100,-30.5,30.5);
- TH2F* h_electrons_pxpy = new TH2F("h_electrons_pxpy", "electron px vs py;px [GeV];py [GeV]", 100,-30.5,30.5,100,-30.5,30.5);
- TH3F* h_electrons_p = new TH3F("h_electron_p", "electron p;px [GeV];py [GeV];pz [GeV]", 100,-30.5,30.5,100,-30.5,30.5,100,-30.5,30.5);
+ TH2F* h_electrons_pzpt = new TH2F("h_electrons_pzpt", "electron pt vs pz;pt [GeV];pz [GeV]", 100,e_start-0.5,e_end+0.5,100,-e_end-0.5,e_end+0.5);
+ TH2F* h_electrons_pxpy = new TH2F("h_electrons_pxpy", "electron px vs py;px [GeV];py [GeV]", 100,-e_end-0.5,e_end+0.5,100,-e_end-0.5,e_end+0.5);
+ TH3F* h_electrons_p = new TH3F("h_electron_p", "electron p;px [GeV];py [GeV];pz [GeV]", 100,-e_end-0.5,e_end+0.5,100,-e_end-0.5,e_end+0.5,100,-e_end-0.5,e_end+0.5);
while(!hepmc_input.failed()) {
// Read event from input file
@@ -69,24 +69,24 @@ void emcal_electrons_reader(){
h_electrons_energy->SetLineWidth(2);
h_electrons_energy->SetLineColor(kBlue);
h_electrons_energy->DrawClone();
- c->SaveAs("results/emcal_electrons_energy_0GeVto30GeV.png");
- c->SaveAs("results/emcal_electrons_energy_0GeVto30GeV.pdf");
+ c->SaveAs("results/input_electrons_energy_0GeVto30GeV.png");
+ c->SaveAs("results/input_electrons_energy_0GeVto30GeV.pdf");
TCanvas* c1 = new TCanvas("c1","c1", 500, 500);
h_electrons_eta->GetYaxis()->SetTitleOffset(1.9);
h_electrons_eta->SetLineWidth(2);
h_electrons_eta->SetLineColor(kBlue);
h_electrons_eta->DrawClone();
- c1->SaveAs("results/emcal_electrons_eta_0GeVto30GeV.png");
- c1->SaveAs("results/emcal_electrons_eta_0GeVto30GeV.pdf");
+ c1->SaveAs("results/input_electrons_eta_0GeVto30GeV.png");
+ c1->SaveAs("results/input_electrons_eta_0GeVto30GeV.pdf");
TCanvas* c2 = new TCanvas("c2","c2", 500, 500);
h_electrons_theta->GetYaxis()->SetTitleOffset(1.8);
h_electrons_theta->SetLineWidth(2);
h_electrons_theta->SetLineColor(kBlue);
h_electrons_theta->DrawClone();
- c2->SaveAs("results/emcal_electrons_theta_0GeVto30GeV.png");
- c2->SaveAs("results/emcal_electrons_theta_0GeVto30GeV.pdf");
+ c2->SaveAs("results/input_electrons_theta_0GeVto30GeV.png");
+ c2->SaveAs("results/input_electrons_theta_0GeVto30GeV.pdf");
TCanvas* c3 = new TCanvas("c3","c3", 500, 500);
h_electrons_phi->GetYaxis()->SetTitleOffset(1.8);
@@ -94,24 +94,24 @@ void emcal_electrons_reader(){
h_electrons_phi->GetYaxis()->SetRangeUser(0.0,h_electrons_phi->GetMaximum()+100.0);
h_electrons_phi->SetLineColor(kBlue);
h_electrons_phi->DrawClone();
- c3->SaveAs("results/emcal_electrons_phi_0GeVto30GeV.png");
- c3->SaveAs("results/emcal_electrons_phi_0GeVto30GeV.pdf");
+ c3->SaveAs("results/input_electrons_phi_0GeVto30GeV.png");
+ c3->SaveAs("results/input_electrons_phi_0GeVto30GeV.pdf");
TCanvas* c4 = new TCanvas("c4","c4", 500, 500);
h_electrons_pzpt->GetYaxis()->SetTitleOffset(1.4);
h_electrons_pzpt->SetLineWidth(2);
h_electrons_pzpt->SetLineColor(kBlue);
h_electrons_pzpt->DrawClone("COLZ");
- c4->SaveAs("results/emcal_electrons_pzpt_0GeVto30GeV.png");
- c4->SaveAs("results/emcal_electrons_pzpt_0GeVto30GeV.pdf");
+ c4->SaveAs("results/input_electrons_pzpt_0GeVto30GeV.png");
+ c4->SaveAs("results/input_electrons_pzpt_0GeVto30GeV.pdf");
TCanvas* c5 = new TCanvas("c5","c5", 500, 500);
h_electrons_pxpy->GetYaxis()->SetTitleOffset(1.4);
h_electrons_pxpy->SetLineWidth(2);
h_electrons_pxpy->SetLineColor(kBlue);
h_electrons_pxpy->DrawClone("COLZ");
- c5->SaveAs("results/emcal_electrons_pxpy_0GeVto30GeV.png");
- c5->SaveAs("results/emcal_electrons_pxpy_0GeVto30GeV.pdf");
+ c5->SaveAs("results/input_electrons_pxpy_0GeVto30GeV.png");
+ c5->SaveAs("results/input_electrons_pxpy_0GeVto30GeV.pdf");
TCanvas* c6 = new TCanvas("c6","c6", 500, 500);
h_electrons_p->GetYaxis()->SetTitleOffset(1.8);
@@ -120,7 +120,7 @@ void emcal_electrons_reader(){
h_electrons_p->SetLineWidth(2);
h_electrons_p->SetLineColor(kBlue);
h_electrons_p->DrawClone();
- c6->SaveAs("results/emcal_electrons_p_0GeVto30GeV.png");
- c6->SaveAs("results/emcal_electrons_p_0GeVto30GeV.pdf");
+ c6->SaveAs("results/input_electrons_p_0GeVto30GeV.png");
+ c6->SaveAs("results/input_electrons_p_0GeVto30GeV.pdf");
}
diff --git a/ecal/scripts/makeplot.C b/ecal/scripts/makeplot.C
index f188b1641ef54ac1ebf9f89e4be0821f78de55d9..e7b88e2611568dbcf8075f0216011964ddd79e67 100644
--- a/ecal/scripts/makeplot.C
+++ b/ecal/scripts/makeplot.C
@@ -2,7 +2,7 @@
// Read reconstruction ROOT output file
// Plot variables
////////////////////////////////////////
-int makeplot(const char* input_fname = "sim_output/sim_emcal_electrons_output.root")
+int makeplot(double e_start = 1.0, double e_end = 1.0, const char* input_fname = "sim_output/sim_emcal_electrons_output.root", const char *fout = "results/rec_1000k_info_1GeV.txt")
{
// Setting figures
gROOT->SetStyle("Plain");
@@ -15,6 +15,11 @@ int makeplot(const char* input_fname = "sim_output/sim_emcal_electrons_output.ro
gStyle->SetPadLeftMargin(0.14);
gStyle->SetPadRightMargin(0.14);
+ // Output ASCII file
+ FILE *fp = fopen(fout,"a");
+ if(!fp)
+ fprintf(stderr,"error: can't open %s\n", fout);
+
// Input ROOT file
TFile *f = new TFile(input_fname,"READ");
TTree *t = (TTree *)f->Get("events");
@@ -66,12 +71,12 @@ int makeplot(const char* input_fname = "sim_output/sim_emcal_electrons_output.ro
TCanvas *c7 = new TCanvas("c7", "c7", 500, 500);
TCanvas *c8 = new TCanvas("c8", "c8", 500, 500);
TCanvas *c9 = new TCanvas("c9", "c9", 500, 500);
-
+
// Declare histograms
TH1D *h1 = new TH1D("h1","Scattering Angle(#theta)", 100,130.0,180.0);
TH1D *h2 = new TH1D("h2","Pseudo-rapidity(#eta)", 100,-5.0,0.0);
- TH2D *h3 = new TH2D("h3","Cluster E vs Pseudo-rapidity", 100,-0.5,30.5,100,-5.0,0.0);
- TH1D *h4 = new TH1D("h4","Reconstructed energy per event", 100,-0.5,30.5);
+ TH2D *h3 = new TH2D("h3","Cluster E vs Pseudo-rapidity", 100,e_start-0.5,e_end+0.5,100,-5.0,0.0);
+ TH1D *h4 = new TH1D("h4","Reconstructed energy per event", 100,e_start-0.5,e_end+0.5);
TH1D *h5 = new TH1D("h5","Number of Clusters per event", 5,-0.5,4.5);
TH1D *h6 = new TH1D("h6","Scattering Angle(#theta) with CUT", 100,130.0,180.0);
TH1D *h7 = new TH1D("h7","Pseudo-rapidity(#eta) with CUT", 100,-5.0,0.0);
@@ -100,7 +105,7 @@ int makeplot(const char* input_fname = "sim_output/sim_emcal_electrons_output.ro
for (int ievent = 0; ievent < nentries; ievent++)
{
t->GetEntry(ievent);
-
+
Int_t ncluster = EcalClusters_;
Int_t nreconhits = RecoEcalHits_;
@@ -120,6 +125,8 @@ int makeplot(const char* input_fname = "sim_output/sim_emcal_electrons_output.ro
cluster_e = cluster_energy[icluster] / 1.e+3;
total_cluster_e += cluster_e;
}
+ // Writing into ASCII file
+ fprintf(fp,"%d\t%d\t%2.4f\n",ievent,ncluster,total_cluster_e);
// Select events with one cluster
if(ncluster == 1)
@@ -141,12 +148,12 @@ int makeplot(const char* input_fname = "sim_output/sim_emcal_electrons_output.ro
for(int ireconhit=0; ireconhit < nreconhits; ireconhit++)
h9->Fill(rec_x_pos[ireconhit],rec_y_pos[ireconhit], 1.0);
}
-
+
// Drawing and Saving figures
c1->cd();
h1->SetLineColor(kBlue);
h1->SetLineWidth(2);
- h1->GetYaxis()->SetRangeUser(0.0,h1->GetMaximum()+10.0);
+ h1->GetYaxis()->SetRangeUser(0.0,h1->GetMaximum()+50.0);
h1->GetXaxis()->SetTitle("#theta [degree]");
h1->GetYaxis()->SetTitle("events");
h1->GetYaxis()->SetTitleOffset(1.4);
@@ -157,7 +164,7 @@ int makeplot(const char* input_fname = "sim_output/sim_emcal_electrons_output.ro
c2->cd();
h2->SetLineColor(kBlue);
h2->SetLineWidth(2);
- h2->GetYaxis()->SetRangeUser(0.0,h2->GetMaximum()+10.0);
+ h2->GetYaxis()->SetRangeUser(0.0,h2->GetMaximum()+50.0);
h2->GetXaxis()->SetTitle("#eta");
h2->GetYaxis()->SetTitle("events");
h2->GetYaxis()->SetTitleOffset(1.4);
@@ -198,7 +205,7 @@ int makeplot(const char* input_fname = "sim_output/sim_emcal_electrons_output.ro
c6->cd();
h6->SetLineColor(kBlue);
h6->SetLineWidth(2);
- h6->GetYaxis()->SetRangeUser(0.0,h1->GetMaximum()+10.0);
+ h6->GetYaxis()->SetRangeUser(0.0,h1->GetMaximum()+50.0);
h6->GetXaxis()->SetTitle("#theta [degree]");
h6->GetYaxis()->SetTitle("events");
h6->GetYaxis()->SetTitleOffset(1.4);
@@ -209,7 +216,7 @@ int makeplot(const char* input_fname = "sim_output/sim_emcal_electrons_output.ro
c7->cd();
h7->SetLineColor(kBlue);
h7->SetLineWidth(2);
- h7->GetYaxis()->SetRangeUser(0.0,h2->GetMaximum()+10.0);
+ h7->GetYaxis()->SetRangeUser(0.0,h2->GetMaximum()+50.0);
h7->GetXaxis()->SetTitle("#eta");
h7->GetYaxis()->SetTitle("events");
h7->GetYaxis()->SetTitleOffset(1.4);
@@ -237,5 +244,6 @@ int makeplot(const char* input_fname = "sim_output/sim_emcal_electrons_output.ro
c9->SaveAs("results/electron_hit_pos_all_0GeVto30GeV.png");
c9->SaveAs("results/electron_hit_pos_all_0GeVto30GeV.pdf");
+ fclose(fp);
return 0;
}
diff --git a/ecal/scripts/makeplot_input.C b/ecal/scripts/makeplot_input.C
new file mode 100644
index 0000000000000000000000000000000000000000..8feab5156e0df7e77032057bf4da14f3bae22e55
--- /dev/null
+++ b/ecal/scripts/makeplot_input.C
@@ -0,0 +1,102 @@
+////////////////////////////////////////
+// Read simulation ROOT input file
+// Write ASCII file
+////////////////////////////////////////
+int makeplot_input(const char* input_fname = "sim_output/sim_emcal_electrons_input.root", const char *fout = "results/sim_1000k_info_1GeV.txt")
+{
+ // Setting figures
+ gROOT->SetStyle("Plain");
+ gStyle->SetLineWidth(3);
+ gStyle->SetOptStat("nem");
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetPadGridX(1);
+ gStyle->SetPadGridY(1);
+ gStyle->SetPadLeftMargin(0.14);
+ gStyle->SetPadRightMargin(0.14);
+
+ // Output ASCII file
+ FILE *fp = fopen(fout,"a");
+ if(!fp)
+ fprintf(stderr,"error: can't open %s\n", fout);
+
+ // Input ROOT file
+ TFile *f = new TFile(input_fname,"READ");
+ TTree *t = (TTree *)f->Get("events");
+
+ // Set Branch status and addressed
+ t->SetMakeClass(1);
+ t->SetBranchStatus("*", 0);
+
+ Int_t CrystalEcalHits_;
+ t->SetBranchStatus("CrystalEcalHits", 1);
+ t->SetBranchAddress("CrystalEcalHits", &CrystalEcalHits_);
+
+ Int_t mcparticles_;
+ t->SetBranchStatus("mcparticles", 1);
+ t->SetBranchAddress("mcparticles", &mcparticles_);
+
+ const Int_t kMaxCrystalEcalHits = 100000;
+ Double_t truth_deposit[kMaxCrystalEcalHits];
+ Double_t energyDeposit[kMaxCrystalEcalHits];
+ t->SetBranchStatus("CrystalEcalHits.truth.deposit",1);
+ t->SetBranchStatus("CrystalEcalHits.energyDeposit",1);
+ t->SetBranchAddress("CrystalEcalHits.truth.deposit",truth_deposit);
+ t->SetBranchAddress("CrystalEcalHits.energyDeposit",energyDeposit);
+
+ const Int_t kMaxmcparticles = 100000;
+ Double_t px[kMaxmcparticles];
+ Double_t py[kMaxmcparticles];
+ Double_t pz[kMaxmcparticles];
+ Double_t mass[kMaxmcparticles];
+ t->SetBranchStatus("mcparticles.psx",1);
+ t->SetBranchStatus("mcparticles.psy",1);
+ t->SetBranchStatus("mcparticles.psz",1);
+ t->SetBranchStatus("mcparticles.mass",1);
+ t->SetBranchAddress("mcparticles.psx",px);
+ t->SetBranchAddress("mcparticles.psy",py);
+ t->SetBranchAddress("mcparticles.psz",pz);
+ t->SetBranchAddress("mcparticles.mass",mass);
+
+ // Total number of entries
+ Int_t nentries = t->GetEntries();
+
+ // Variables are used in calculation
+ Double_t total_thr_e; // Thrown energy [GeV]
+ Double_t truth_sim_e; // truth simulated deposit energy [GeV]
+ Double_t total_truth_sim_e; // Add up truth simulated deposit energy [GeV]
+ Double_t sim_e; // simulated deposit energy [GeV]
+ Double_t total_sim_e; // Add up truth simulated deposit energy [GeV]
+ Double_t momentum2; // momentum squared
+ Double_t mass2; // mass squared
+
+ // Loop over event by event
+ for (int ievent = 0; ievent < nentries; ievent++)
+ {
+ t->GetEntry(ievent);
+
+ Int_t nCrystalEcalHits = CrystalEcalHits_;
+ Int_t nmcparticle = 2;
+
+ total_thr_e = 0.0;
+ total_truth_sim_e = 0.0;
+ total_sim_e = 0.0;
+
+ momentum2 = px[nmcparticle]*px[nmcparticle]+py[nmcparticle]*py[nmcparticle]+pz[nmcparticle]*pz[nmcparticle];
+ mass2 = mass[nmcparticle]*mass[nmcparticle];
+
+ total_thr_e = sqrt(momentum2 + mass2)/1.e+3;
+
+ // Loop over cluster by cluster
+ for (int isimhit=0; isimhit < nCrystalEcalHits; isimhit++)
+ {
+ truth_sim_e = truth_deposit[isimhit];
+ sim_e = energyDeposit[isimhit];
+ total_truth_sim_e += truth_sim_e;
+ total_sim_e += energyDeposit[isimhit]/1.e+3;
+ }
+ fprintf(fp,"%d\t%2.4f\t%2.4f\t%2.4f\n",ievent,total_thr_e,total_truth_sim_e,total_sim_e);
+ }
+ fclose(fp);
+ return 0;
+}
diff --git a/ecal/scripts/read_eng.C b/ecal/scripts/read_eng.C
new file mode 100644
index 0000000000000000000000000000000000000000..cedc16dad4d8837e7473bd0a48eecadeeab2e0af
--- /dev/null
+++ b/ecal/scripts/read_eng.C
@@ -0,0 +1,70 @@
+int read_eng(const char* output_fname = "sim_output/read_eng_output.root", const char *fname = "results/eng.txt")
+{
+ // Input ASCII File
+ FILE *fin = fopen(fname,"r");
+
+ if(!fin)
+ fprintf(stderr,"error: can't open %s\n", fname);
+
+ // Output ROOT File
+ TFile *rootFile = new TFile(output_fname,"recreate");
+ // ROOT tree that goes into this ROOT file
+ TTree *tEngRes = new TTree("tEngRes","Energy Resolution Information");
+ // Variables that go into ROOT file
+ Int_t ievent;
+ Double_t tot_clust_eng; // [GeV]
+ Int_t ncluster;
+ Double_t mc_eng; // [GeV]
+ Double_t sim_tru_eng; // [GeV]
+ Double_t sim_eng; // [GeV]
+ // Declare ROOT tree branches
+ tEngRes->Branch("ievent", &ievent, "ievent/I");
+ tEngRes->Branch("tot_clust_eng", &tot_clust_eng, "tot_clust_eng/D");
+ tEngRes->Branch("ncluster", &ncluster, "ncluster/I");
+ tEngRes->Branch("mc_eng", &mc_eng, "mc_eng/D");
+ tEngRes->Branch("sim_tru_eng", &sim_tru_eng, "sim_tru_eng/D");
+ tEngRes->Branch("sim_eng", &sim_eng, "sim_eng/D");
+
+ char inbuf[1024];
+ int nevents_readout = 0;
+ int linenum = 0;
+
+ TString s;
+ TObjArray *o = 0;
+
+ // Read each line
+ while(fgets(inbuf,1024,fin))
+ {
+ nevents_readout++;
+ s = inbuf;
+ o = s.Tokenize(" \t\n");
+ for(int icol = 0; icol < o->GetEntries(); icol++)
+ {
+ s = ((TObjString *)(*o)[icol])->GetString();
+ if(icol==0)
+ sscanf(s.Data(),"%d", &ievent);
+ if(icol==1)
+ sscanf(s.Data(),"%lf", &mc_eng);
+ if(icol==2)
+ sscanf(s.Data(),"%lf", &sim_tru_eng);
+ if(icol==3)
+ sscanf(s.Data(),"%lf", &sim_eng);
+ if(icol==5)
+ sscanf(s.Data(),"%d", &ncluster);
+ if(icol==6)
+ sscanf(s.Data(),"%lf", &tot_clust_eng);
+ }
+ tEngRes->Fill();
+ // Clean up
+ delete o;
+ }
+ // Done reading ASCII file. So finialize the ROOT tree and close the ROOT file
+ fprintf(stdout,"%d events read in!\n",nevents_readout);
+ // Finalize ROOT file
+ tEngRes->Write();
+ // Close ROOT file
+ rootFile->Close();
+ fclose(fin);
+
+ return 0;
+}