Skip to content
Snippets Groups Projects
Commit fba70eca authored by Jihee Kim's avatar Jihee Kim
Browse files

Included mcparticles2 and CrystalEcalHits2 branches info

parent c3aa6ddf
No related branches found
No related tags found
1 merge request!15Acceptance
This commit is part of merge request !15. Comments created here will be created in the context of that merge request.
Hide whitespace changes
Inline Side-by-side
ecal/scripts/rec_emcal_electrons_reader.C
+ 39
− 48
View file @ fba70eca
......@@ -2,7 +2,6 @@
// Read reconstruction ROOT output file
// Plot variables
////////////////////////////////////////
//int rec_emcal_electrons_reader(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")
int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const char* input_fname = "sim_output/sim_emcal_electrons_output.root")
{
// Setting figures
......@@ -16,11 +15,6 @@ int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const c
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");
......@@ -29,37 +23,37 @@ int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const c
t->SetMakeClass(1);
t->SetBranchStatus("*", 0);
Int_t mcparticles2_;
//t->SetBranchStatus("mcparticles2", 1);
//t->SetBranchAddress("mcparticles2", &mcparticles2_);
//Int_t CrystalEcalHits2_;
//t->SetBranchStatus("CrystalEcalHits2", 1);
//t->SetBranchAddress("CrystalEcalHits2", &CrystalEcalHits2_);
t->SetBranchStatus("mcparticles2", 1);
t->SetBranchAddress("mcparticles2", &mcparticles2_);
Int_t CrystalEcalHits2_;
t->SetBranchStatus("CrystalEcalHits2", 1);
t->SetBranchAddress("CrystalEcalHits2", &CrystalEcalHits2_);
Int_t RecoEcalHits_;
t->SetBranchStatus("RecoEcalHits", 1);
t->SetBranchAddress("RecoEcalHits", &RecoEcalHits_);
Int_t EcalClusters_;
t->SetBranchStatus("EcalClusters", 1);
t->SetBranchAddress("EcalClusters", &EcalClusters_);
//const Int_t kMaxmcparticles2 = 100000;
//Double_t px[kMaxmcparticles2];
//Double_t py[kMaxmcparticles2];
//Double_t pz[kMaxmcparticles2];
//Double_t mass[kMaxmcparticles2];
//t->SetBranchStatus("mcparticles2.psx",1);
//t->SetBranchStatus("mcparticles2.psy",1);
//t->SetBranchStatus("mcparticles2.psz",1);
//t->SetBranchStatus("mcparticles2.mass",1);
//t->SetBranchAddress("mcparticles2.psx",px);
//t->SetBranchAddress("mcparticles2.psy",py);
//t->SetBranchAddress("mcparticles2.psz",pz);
//t->SetBranchAddress("mcparticles2.mass",mass);
//const Int_t kMaxCrystalEcalHits2 = 100000;
//Double_t truth_deposit[kMaxCrystalEcalHits2];
//Double_t energyDeposit[kMaxCrystalEcalHits2];
//t->SetBranchStatus("CrystalEcalHits2.truth.deposit",1);
//t->SetBranchStatus("CrystalEcalHits2.energyDeposit",1);
//t->SetBranchAddress("CrystalEcalHits2.truth.deposit",truth_deposit);
//t->SetBranchAddress("CrystalEcalHits2.energyDeposit",energyDeposit);
const Int_t kMaxmcparticles2 = 100000;
Double_t px[kMaxmcparticles2];
Double_t py[kMaxmcparticles2];
Double_t pz[kMaxmcparticles2];
Double_t mass[kMaxmcparticles2];
t->SetBranchStatus("mcparticles2.psx",1);
t->SetBranchStatus("mcparticles2.psy",1);
t->SetBranchStatus("mcparticles2.psz",1);
t->SetBranchStatus("mcparticles2.mass",1);
t->SetBranchAddress("mcparticles2.psx",px);
t->SetBranchAddress("mcparticles2.psy",py);
t->SetBranchAddress("mcparticles2.psz",pz);
t->SetBranchAddress("mcparticles2.mass",mass);
const Int_t kMaxCrystalEcalHits2 = 100000;
Double_t truth_deposit[kMaxCrystalEcalHits2];
Double_t energyDeposit[kMaxCrystalEcalHits2];
t->SetBranchStatus("CrystalEcalHits2.truth.deposit",1);
t->SetBranchStatus("CrystalEcalHits2.energyDeposit",1);
t->SetBranchAddress("CrystalEcalHits2.truth.deposit",truth_deposit);
t->SetBranchAddress("CrystalEcalHits2.energyDeposit",energyDeposit);
const Int_t kMaxRecoEcalHits = 100000;
Double_t rec_x_pos[kMaxRecoEcalHits];
Double_t rec_y_pos[kMaxRecoEcalHits];
......@@ -147,8 +141,8 @@ int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const c
// Read event by event
t->GetEntry(ievent);
// Read number of hits/clusters
//Int_t nmcparticle = 2;
//Int_t nCrystalEcalHits = CrystalEcalHits2_;
Int_t nmcparticle = 2;
Int_t nCrystalEcalHits = CrystalEcalHits2_;
Int_t nreconhits = RecoEcalHits_;
Int_t ncluster = EcalClusters_;
// Initialize total energy variables
......@@ -157,18 +151,18 @@ int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const c
total_sim_e = 0.0;
total_cluster_e = 0.0;
// Thrown energy, momentum, and mass
//momentum2 = px[nmcparticle]*px[nmcparticle]+py[nmcparticle]*py[nmcparticle]+pz[nmcparticle]*pz[nmcparticle];
//momentum = TMath::Sqrt(momentum2);
//mass2 = mass[nmcparticle]*mass[nmcparticle];
//total_thr_e = sqrt(momentum2 + mass2)/1.e+3;
//h12->Fill(momentum, 1.0);
momentum2 = px[nmcparticle]*px[nmcparticle]+py[nmcparticle]*py[nmcparticle]+pz[nmcparticle]*pz[nmcparticle];
momentum = TMath::Sqrt(momentum2);
mass2 = mass[nmcparticle]*mass[nmcparticle];
total_thr_e = sqrt(momentum2 + mass2)/1.e+3;
h12->Fill(momentum, 1.0);
// Loop over simulated hit by simulated hit
//for (int isimhit=0; isimhit < nCrystalEcalHits; isimhit++)
//{
// total_truth_sim_e += truth_deposit[isimhit];
// total_sim_e += energyDeposit[isimhit]/1.e+3;
//}
for (int isimhit=0; isimhit < nCrystalEcalHits; isimhit++)
{
total_truth_sim_e += truth_deposit[isimhit];
total_sim_e += energyDeposit[isimhit]/1.e+3;
}
// Loop over reconstructed hit by reconstructed hit
for(int ireconhit=0; ireconhit < nreconhits; ireconhit++)
......@@ -188,8 +182,6 @@ int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const c
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)
......@@ -212,8 +204,8 @@ int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const c
h11->Fill(eng_res, 1.0);
}
//if(total_cluster_e > 0.9*total_thr_e)
// h13->Fill(momentum, 1.0);
if(total_cluster_e > 0.9*total_thr_e)
h13->Fill(momentum, 1.0);
}
}
// Drawing and Saving figures
......@@ -366,6 +358,5 @@ int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const c
c14->SaveAs("results/hAcceptance.png");
c14->SaveAs("results/hAcceptance.pdf");
//fclose(fp);
return 0;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment