Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
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;
}