-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathanalyzeFF_eicrecon.C
356 lines (255 loc) · 14 KB
/
analyzeFF_eicrecon.C
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
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
//-------------------------
//
// Simple analysis code to analyze EPIC simulation output for Roman Pots, OMD, B0, and ZDC.
//
// NOTES:
//
// 1) B0 works out of the box for charged particles - the tracks are reconstrucuted using ACTS with truth
// seeding and stored as ReconstructedChargedParticles.
//
// 2) RP + OMD are using transport matrix reconstruction with a central-orbit matrix, and assuming decoupled
// x and y motion (which is true for central orbits, we only steer in x).
// This means that as XL increasingly moves away from 1, the matrix-induced smearing is increased. This
// will be fixed in the future, but is not analysis breaking.
//
// For the OMD, it's a bit worse because the transport for off-momentum protons through the magnets is MUCH
// messier. So if you are doing tagged-DIS, you may need to watch out for particle at large angles (theta > 1mrad),
// and for XL != 0.5 (which is the assumed OMD "central orbit" for calculating the matrix).
//
// 3) For ZDC, no full reco is currently available, so only acceptance studies with poorly calibrated energy information
// is available.
//
//-------------------------------------------------
// Please email [email protected] with any issues.
//-------------------------------------------------
//
// Author: Alex Jentsch
//
// Date of last author update: 8/21/2023
//
//------------------------
using namespace std;
void analyzeFF_eicrecon(){
TString fileList = "./inputFileList_particleGun.list";
TString outputName = "ePIC_fullReco_LOCAL_COORD_RP_Output_";
TString date = "8_21_2023_";
TString run = "run_0";
cout << "Input FileList: " << fileList << endl;
TString fileType_ROOT = ".root";
TString outputFileName = outputName + date + run + fileType_ROOT;
string fileName;
TFile * inputRootFile;
TTree * rootTree;
cout << "Output file: " << outputFileName << endl;
ifstream fileListStream;
fileListStream.open(fileList);
if(!fileListStream) { cout << "NO_LIST_FILE " << fileList << endl; return;}
//--------------------------------------------------------------------------
//histograms -- only a few for now
//MC information
TH1D* h_eta_MC = new TH1D("h_eta",";Pseudorapidity, #eta",100,0.0,15.0);
TH1D* h_px_MC = new TH1D("px_MC", ";p_{x} [GeV/c]", 100, -10.0, 10.0);
TH1D* h_py_MC = new TH1D("py_MC", ";p_{y} [GeV/c]", 100, -10.0, 10.0);
TH1D* h_pt_MC = new TH1D("pt_MC", ";p_{t} [GeV/c]", 100, 0.0, 2.0);
TH1D* h_pz_MC = new TH1D("pz_MC", ";p_{z} [GeV/c]", 100, 0.0, 320.0);
//Roman pots
TH1D* h_px_RomanPots = new TH1D("px_RomanPots", ";p_{x} [GeV/c]", 100, -10.0, 10.0);
TH1D* h_py_RomanPots = new TH1D("py_RomanPots", ";p_{y} [GeV/c]", 100, -10.0, 10.0);
TH1D* h_pt_RomanPots = new TH1D("pt_RomanPots", ";p_{t} [GeV/c]", 100, 0.0, 2.0);
TH1D* h_pz_RomanPots = new TH1D("pz_RomanPots", ";p_{z} [GeV/c]", 100, 0.0, 320.0);
TH2D* h_rp_occupancy_map = new TH2D("Roman_pots_occupancy_map", "hit y [mm];hit x [mm]", 100, -150, 150, 100, -70, 70);
//OMD
TH1D* h_px_OMD = new TH1D("px_OMD", ";p_{x} [GeV/c]", 100, -10.0, 10.0);
TH1D* h_py_OMD = new TH1D("py_OMD", ";p_{y} [GeV/c]", 100, -10.0, 10.0);
TH1D* h_pt_OMD = new TH1D("pt_OMD", ";p_{t} [GeV/c]", 100, 0.0, 2.0);
TH1D* h_pz_OMD = new TH1D("pz_OMD", ";p_{z} [GeV/c]", 100, 0.0, 320.0);
TH2D* h_omd_occupancy_map = new TH2D("OMD_occupancy_map", "hit y [mm];hit x [mm]", 100, -150, 150, 100, -70, 70);
//B0 tracker hits
TH2D* h_B0_occupancy_map_layer_0 = new TH2D("B0_occupancy_map_0", "B0_occupancy_map_0", 100, -400, 0, 100, -170, 170);
TH2D* h_B0_occupancy_map_layer_1 = new TH2D("B0_occupancy_map_1", "B0_occupancy_map_1", 100, -400, 0, 100, -170, 170);
TH2D* h_B0_occupancy_map_layer_2 = new TH2D("B0_occupancy_map_2", "B0_occupancy_map_2", 100, -400, 0, 100, -170, 170);
TH2D* h_B0_occupancy_map_layer_3 = new TH2D("B0_occupancy_map_3", "B0_occupancy_map_3", 100, -400, 0, 100, -170, 170);
TH1D* h_B0_hit_energy_deposit = new TH1D("B0_tracker_hit_energy_deposit", ";Deposited Energy [keV]", 100, 0.0, 500.0);
//B0 EMCAL clusters
TH2D* h_B0_emcal_occupancy_map = new TH2D("B0_emcal_occupancy_map", "B0_emcal_occupancy_map", 100, -400, 0, 100, -170, 170);
TH1D* h_B0_emcal_cluster_energy = new TH1D("B0_emcal_cluster_energy", ";Cluster Energy [GeV]", 100, 0.0, 100.0);
//Reconstructed tracks (for usage with B0 too!!)
TH1D* h_px_reco_track = new TH1D("px_reco_track", ";p_{x} [GeV/c]", 100, -10.0, 10.0);
TH1D* h_py_reco_track = new TH1D("py_reco_track", ";p_{y} [GeV/c]", 100, -10.0, 10.0);
TH1D* h_pt_reco_track = new TH1D("pt_reco_track", ";p_{t} [GeV/c]", 100, 0.0, 2.0);
TH1D* h_pz_reco_track = new TH1D("pz_reco_track", ";p_{z} [GeV/c]", 100, 0.0, 320.0);
//ZDC EMCAL clusters
TH2D* h_ZDC_emcal_occupancy_map = new TH2D("ZDC_emcal_occupancy_map", "ZDC_emcal_occupancy_map", 100, -1150, -1050, 100, -60, 60);
TH1D* h_ZDC_emcal_cluster_energy = new TH1D("ZDC_emcal_cluster_energy", ";Cluster Energy [GeV]", 100, 0.0, 100.0);
int fileCounter = 0;
int iEvent = 0;
while(getline(fileListStream, fileName)){
TString tmp = fileName;
cout << "Input file " << fileCounter << ": " << fileName << endl;
inputRootFile = new TFile(tmp);
if(!inputRootFile){ cout << "MISSING_ROOT_FILE"<< fileName << endl; continue;}
fileCounter++;
TTree * evtTree = (TTree*)inputRootFile->Get("events");
int numEvents = evtTree->GetEntries();
TTreeReader tree_reader(evtTree); // !the tree reader
//MC particles
TTreeReaderArray<float> mc_px_array = {tree_reader, "MCParticles.momentum.x"};
TTreeReaderArray<float> mc_py_array = {tree_reader, "MCParticles.momentum.y"};
TTreeReaderArray<float> mc_pz_array = {tree_reader, "MCParticles.momentum.z"};
TTreeReaderArray<double> mc_mass_array = {tree_reader, "MCParticles.mass"};
TTreeReaderArray<int> mc_pdg_array = {tree_reader, "MCParticles.PDG"};
//Roman pots -- momentum vector
TTreeReaderArray<float> reco_RP_px = {tree_reader, "ForwardRomanPotRecParticles.momentum.x"};
TTreeReaderArray<float> reco_RP_py = {tree_reader, "ForwardRomanPotRecParticles.momentum.y"};
TTreeReaderArray<float> reco_RP_pz = {tree_reader, "ForwardRomanPotRecParticles.momentum.z"};
//Off-Momentum -- momentum vector
TTreeReaderArray<float> reco_OMD_px = {tree_reader, "ForwardOffMRecParticles.momentum.x"};
TTreeReaderArray<float> reco_OMD_py = {tree_reader, "ForwardOffMRecParticles.momentum.y"};
TTreeReaderArray<float> reco_OMD_pz = {tree_reader, "ForwardOffMRecParticles.momentum.z"};
//hit locations (for debugging)
TTreeReaderArray<float> global_hit_RP_x = {tree_reader, "ForwardRomanPotRecParticles.referencePoint.x"};
TTreeReaderArray<float> global_hit_RP_y = {tree_reader, "ForwardRomanPotRecParticles.referencePoint.y"};
TTreeReaderArray<float> global_hit_RP_z = {tree_reader, "ForwardRomanPotRecParticles.referencePoint.z"};
//hit locations (for debugging)
TTreeReaderArray<float> global_hit_OMD_x = {tree_reader, "ForwardOffMRecParticles.referencePoint.x"};
TTreeReaderArray<float> global_hit_OMD_y = {tree_reader, "ForwardOffMRecParticles.referencePoint.y"};
TTreeReaderArray<float> global_hit_OMD_z = {tree_reader, "ForwardOffMRecParticles.referencePoint.z"};
//b0 tracker hits
TTreeReaderArray<float> b0_hits_x = {tree_reader, "B0TrackerRecHits.position.x"};
TTreeReaderArray<float> b0_hits_y = {tree_reader, "B0TrackerRecHits.position.y"};
TTreeReaderArray<float> b0_hits_z = {tree_reader, "B0TrackerRecHits.position.z"};
TTreeReaderArray<float> b0_hits_eDep = {tree_reader, "B0TrackerRecHits.edep"}; //deposited energy per hit
//b0 EMCAL
TTreeReaderArray<float> b0_cluster_x = {tree_reader, "B0ECalClusters.position.x"};
TTreeReaderArray<float> b0_cluster_y = {tree_reader, "B0ECalClusters.position.y"};
TTreeReaderArray<float> b0_cluster_z = {tree_reader, "B0ECalClusters.position.z"};
TTreeReaderArray<float> b0_cluster_energy = {tree_reader, "B0ECalClusters.energy"}; //deposited energy in cluster
//reco tracks (where b0 tracks live!!!)
TTreeReaderArray<float> reco_track_x = {tree_reader, "ReconstructedChargedParticles.momentum.x"};
TTreeReaderArray<float> reco_track_y = {tree_reader, "ReconstructedChargedParticles.momentum.y"};
TTreeReaderArray<float> reco_track_z = {tree_reader, "ReconstructedChargedParticles.momentum.z"};
//ZDC EMCAL
TTreeReaderArray<float> zdc_ecal_cluster_x = {tree_reader, "ZDCEcalClusters.position.x"};
TTreeReaderArray<float> zdc_ecal_cluster_y = {tree_reader, "ZDCEcalClusters.position.y"};
TTreeReaderArray<float> zdc_ecal_cluster_z = {tree_reader, "ZDCEcalClusters.position.z"};
TTreeReaderArray<float> zdc_ecal_cluster_energy = {tree_reader, "ZDCEcalClusters.energy"}; //deposited energy in cluster
cout << "file has " << evtTree->GetEntries() << " events..." << endl;
tree_reader.SetEntriesRange(0, evtTree->GetEntries());
while (tree_reader.Next()) {
cout << "Reading event: " << iEvent << endl;
//MCParticles
//finding the far-forward proton;
//TLorentzVector scatMC(0,0,0,0);
TVector3 mctrk;
TVector3 rptrk;
double maxPt=-99.;
for(int imc=0;imc<mc_px_array.GetSize();imc++){
mctrk.SetXYZ(mc_px_array[imc], mc_py_array[imc], mc_pz_array[imc]);
if(mc_pdg_array[imc] == 2212){ //only checking for protons here -- change as desired
mctrk.RotateY(0.025);
h_eta_MC->Fill(mctrk.Eta());
h_px_MC->Fill(mctrk.Px());
h_py_MC->Fill(mctrk.Py());
h_pt_MC->Fill(mctrk.Perp());
h_pz_MC->Fill(mctrk.Pz());
}
}
//roman pots reco tracks
for(int iRPPart = 0; iRPPart < reco_RP_px.GetSize(); iRPPart++){
TVector3 prec_romanpots(reco_RP_px[iRPPart], reco_RP_py[iRPPart], reco_RP_pz[iRPPart]);
h_px_RomanPots->Fill(prec_romanpots.Px());
h_py_RomanPots->Fill(prec_romanpots.Py());
h_pt_RomanPots->Fill(prec_romanpots.Perp());
h_pz_RomanPots->Fill(prec_romanpots.Pz());
h_rp_occupancy_map->Fill(global_hit_RP_x[iRPPart], global_hit_RP_y[iRPPart]);
}
//OMD reco tracks
for(int iOMDPart = 0; iOMDPart < reco_OMD_px.GetSize(); iOMDPart++){
TVector3 prec_omd(reco_OMD_px[iOMDPart], reco_OMD_py[iOMDPart], reco_OMD_pz[iOMDPart]);
h_px_OMD->Fill(prec_omd.Px());
h_py_OMD->Fill(prec_omd.Py());
h_pt_OMD->Fill(prec_omd.Perp());
h_pz_OMD->Fill(prec_omd.Pz());
h_omd_occupancy_map->Fill(global_hit_OMD_x[iOMDPart], global_hit_OMD_y[iOMDPart]);
}
double hit_x = -9999.;
double hit_y = -9999.;
double hit_z = -9999.;
double hit_deposited_energy = -9999.;
for(int b0cluster = 0; b0cluster < b0_cluster_x.GetSize(); b0cluster++){
hit_x = b0_cluster_x[b0cluster];
hit_y = b0_cluster_y[b0cluster];
hit_z = b0_cluster_z[b0cluster];
hit_deposited_energy = b0_cluster_energy[b0cluster]*1.246; //poor man's calibration constant, for now
h_B0_emcal_occupancy_map->Fill(hit_x, hit_y);
h_B0_emcal_cluster_energy->Fill(hit_deposited_energy);
}
//b0 tracker hits -- for debugging or external tracking
for(int b0hit = 0; b0hit < b0_hits_x.GetSize(); b0hit++){
hit_x = b0_hits_x[b0hit];
hit_y = b0_hits_y[b0hit];
hit_z = b0_hits_z[b0hit];
hit_deposited_energy = b0_hits_eDep[b0hit]*1e6; //convert GeV --> keV
h_B0_hit_energy_deposit->Fill(hit_deposited_energy);
if(hit_deposited_energy < 10.0){ continue; } //threshold value -- 10 keV, arbitrary, for now
//ACLGAD layout
if(hit_z > 5700 && hit_z < 5990){ h_B0_occupancy_map_layer_0->Fill(hit_x, hit_y); }
if(hit_z > 6100 && hit_z < 6200){ h_B0_occupancy_map_layer_1->Fill(hit_x, hit_y); }
if(hit_z > 6400 && hit_z < 6500){ h_B0_occupancy_map_layer_2->Fill(hit_x, hit_y); }
if(hit_z > 6700 && hit_z < 6750){ h_B0_occupancy_map_layer_3->Fill(hit_x, hit_y); }
}
//reconstructed tracks with ACTS -- used for B0
for(int iRecoTrk = 0; iRecoTrk < reco_track_x.GetSize(); iRecoTrk++){
TVector3 prec_reco_tracks(reco_track_x[iRecoTrk], reco_track_y[iRecoTrk], reco_track_z[iRecoTrk]);
prec_reco_tracks.RotateY(0.025); //remove crossing angle for B0!!!
h_px_reco_track->Fill(prec_reco_tracks.Px());
h_py_reco_track->Fill(prec_reco_tracks.Py());
h_pt_reco_track->Fill(prec_reco_tracks.Perp());
h_pz_reco_track->Fill(prec_reco_tracks.Pz());
}
for(int iZdcEMCALcluster = 0; iZdcEMCALcluster < zdc_ecal_cluster_x.GetSize(); iZdcEMCALcluster++){
hit_x = zdc_ecal_cluster_x[iZdcEMCALcluster];
hit_y = zdc_ecal_cluster_y[iZdcEMCALcluster];
hit_z = zdc_ecal_cluster_z[iZdcEMCALcluster];
hit_deposited_energy = zdc_ecal_cluster_energy[iZdcEMCALcluster]*1.246; //poor man's calibration constant, for now
h_ZDC_emcal_occupancy_map->Fill(hit_x, hit_y);
h_ZDC_emcal_cluster_energy->Fill(hit_deposited_energy);
}
iEvent++;
}// event loop
inputRootFile->Close();
}// input file loop
cout << "Check integrals: " << endl;
cout << "pt_mc integral = " << h_pt_MC->Integral() << endl;
cout << "pt_RP_reco integral = " << h_pt_RomanPots->Integral() << endl;
TFile * outputFile = new TFile(outputFileName, "RECREATE");
h_px_MC->Write();
h_py_MC->Write();
h_pt_MC->Write();
h_pz_MC->Write();
h_px_RomanPots->Write();
h_py_RomanPots->Write();
h_pt_RomanPots->Write();
h_pz_RomanPots->Write();
h_rp_occupancy_map->Write();
h_px_OMD->Write();
h_py_OMD->Write();
h_pt_OMD->Write();
h_pz_OMD->Write();
h_omd_occupancy_map->Write();
h_B0_occupancy_map_layer_0->Write();
h_B0_occupancy_map_layer_1->Write();
h_B0_occupancy_map_layer_2->Write();
h_B0_occupancy_map_layer_3->Write();
h_B0_hit_energy_deposit->Write();
h_B0_emcal_occupancy_map->Write();
h_B0_emcal_cluster_energy->Write();
h_px_reco_track->Write();
h_py_reco_track->Write();
h_pt_reco_track->Write();
h_pz_reco_track->Write();
h_ZDC_emcal_occupancy_map->Write();
h_ZDC_emcal_cluster_energy->Write();
outputFile->Close();
return;
}