1 |
pam-fi |
1.1 |
//////////////////////////////////////////////////////////////// |
2 |
|
|
// |
3 |
|
|
// Histogram set for tracker-efficiency study. |
4 |
|
|
// |
5 |
|
|
// It contains the functions: |
6 |
|
|
// |
7 |
|
|
// - CreateHistos() |
8 |
|
|
// - FillHistos(PamLevel2* event) |
9 |
|
|
// - Save Histos() |
10 |
|
|
// |
11 |
|
|
// Usage (example): |
12 |
|
|
// [] .L My-histo-set-macro.C+ |
13 |
|
|
// [] .L My-selection.C+ |
14 |
|
|
// [] .L Loop.C+ |
15 |
|
|
// [] Loop("data-dir/","file-list.txt",100000,"+ALL fillTree fillHistos","output-file.root") |
16 |
|
|
// |
17 |
|
|
//////////////////////////////////////////////////////////////// |
18 |
|
|
|
19 |
|
|
#if !defined(__CINT__) || defined(__MAKECINT__) |
20 |
|
|
|
21 |
|
|
#include <TString.h> |
22 |
|
|
#include <TH1F.h> |
23 |
|
|
#include <TH2F.h> |
24 |
|
|
#include <TMath.h> |
25 |
|
|
#include <TLine.h> |
26 |
|
|
#include <TPolyMarker.h> |
27 |
|
|
#include <TSelector.h> |
28 |
|
|
#include <TFile.h> |
29 |
|
|
|
30 |
|
|
#include <stdlib.h> |
31 |
|
|
#include <iostream> |
32 |
|
|
using namespace std; |
33 |
|
|
|
34 |
|
|
#include <PamLevel2.h> |
35 |
|
|
#include <CaloAxis.h> |
36 |
|
|
|
37 |
|
|
#endif |
38 |
|
|
//====================== |
39 |
|
|
// HITOGRAMS declaration |
40 |
|
|
//====================== |
41 |
|
|
TH1F *htofresx; |
42 |
|
|
TH1F *htofresy; |
43 |
|
|
TH1F *hcalresx; |
44 |
|
|
TH1F *hcalresy; |
45 |
|
|
TH1F *htrkresx[10]; |
46 |
|
|
TH1F *htrkresy[10]; |
47 |
|
|
TH1F *hxcalotrk[6]; |
48 |
|
|
TH1F *hycalotrk[6]; |
49 |
|
|
TH1F *hresangx; |
50 |
|
|
TH1F *hresangy; |
51 |
|
|
//=============================================================== |
52 |
|
|
// Create histograms |
53 |
|
|
//=============================================================== |
54 |
|
|
void CreateHistos( ){ |
55 |
|
|
|
56 |
|
|
gROOT->cd();//create histos in memory |
57 |
|
|
|
58 |
|
|
gDirectory->Delete("htofresx"); |
59 |
|
|
gDirectory->Delete("htofresy"); |
60 |
|
|
TH1F *htofresx = new TH1F("htofresx","Spatial residuals (S12+S32) - S21 (x)",200,-10.,10.); |
61 |
|
|
TH1F *htofresy = new TH1F("htofresy","Spatial residuals (S11+S31) - S22 (y)",200,-10.,10.); |
62 |
|
|
|
63 |
|
|
|
64 |
|
|
gDirectory->Delete("hcalresx"); |
65 |
|
|
gDirectory->Delete("hcalresy"); |
66 |
|
|
TH1F *hcalresx = new TH1F("hcalresx","Spatial residuals on CAL - S32 (x)",200,-10.,10.); |
67 |
|
|
TH1F *hcalresy = new TH1F("hcalresy","Spatial residuals on CAL - S31 (y)",200,-10.,10.); |
68 |
|
|
|
69 |
|
|
TString hid; |
70 |
|
|
for(Int_t i=0; i<3; i++){ |
71 |
|
|
hid.Form("htrkresx%i",i); |
72 |
|
|
htrkresx[i] = new TH1F(hid.Data(),"Spatial residuals on TRK - TOF (x)",200,-10.,10.); |
73 |
|
|
hid.Form("htrkresy%i",i); |
74 |
|
|
htrkresy[i] = new TH1F(hid.Data(),"Spatial residuals on TRK - TOF (y)",200,-10.,10.); |
75 |
|
|
} |
76 |
|
|
for(Int_t i=3; i<9; i++){ |
77 |
|
|
hid.Form("htrkresx%i",i); |
78 |
|
|
htrkresx[i] = new TH1F(hid.Data(),"Spatial residuals on TRK (x)",5000,-.2,.2); |
79 |
|
|
hid.Form("htrkresy%i",i); |
80 |
|
|
htrkresy[i] = new TH1F(hid.Data(),"Spatial residuals on TRK (y)",5000,-.2,.2); |
81 |
|
|
} |
82 |
|
|
for(Int_t i=9; i<10; i++){ |
83 |
|
|
hid.Form("htrkresx%i",i); |
84 |
|
|
htrkresx[i] = new TH1F(hid.Data(),"Spatial residuals on TRK - CAL (x)",800,-10.,10.); |
85 |
|
|
hid.Form("htrkresy%i",i); |
86 |
|
|
htrkresy[i] = new TH1F(hid.Data(),"Spatial residuals on TRK - CAL (y)",800,-10.,10.); |
87 |
|
|
} |
88 |
|
|
|
89 |
|
|
|
90 |
|
|
for(Int_t ip=0; ip<6; ip++){ |
91 |
|
|
hid.Form("hxcalotrk%i",ip); |
92 |
|
|
hxcalotrk[ip] = new TH1F(hid.Data(),"Calo axis intersaction on Trk plane (x)",500,-25.,25.); |
93 |
|
|
hid.Form("hycalotrk%i",ip); |
94 |
|
|
hycalotrk[ip] = new TH1F(hid.Data(),"Calo axis intersaction on Trk plane (y)",500,-25.,25.); |
95 |
|
|
} |
96 |
|
|
|
97 |
|
|
hresangx = new TH1F("hresangx","Angular residual on the top calo plane (x)",400,-40.,40.); |
98 |
|
|
hresangy = new TH1F("hresangy","Angular residual on the top calo plane (y)",400,-40.,40.); |
99 |
|
|
|
100 |
|
|
} |
101 |
|
|
|
102 |
|
|
//=============================================================== |
103 |
|
|
// Fill histograms (called event-by-event) |
104 |
|
|
//=============================================================== |
105 |
|
|
|
106 |
|
|
void FillHistos( PamLevel2* event ){ |
107 |
|
|
// ---------------------------------- |
108 |
|
|
// retrieve histos |
109 |
|
|
// ---------------------------------- |
110 |
|
|
htofresx = dynamic_cast<TH1F*>(gDirectory->FindObject("htofresx")); |
111 |
|
|
htofresy = dynamic_cast<TH1F*>(gDirectory->FindObject("htofresy")); |
112 |
|
|
// cout << hex << htofresx << dec << endl; |
113 |
|
|
hcalresx = dynamic_cast<TH1F*>(gDirectory->FindObject("hcalresx")); |
114 |
|
|
hcalresy = dynamic_cast<TH1F*>(gDirectory->FindObject("hcalresy")); |
115 |
|
|
TString hid; |
116 |
|
|
for(Int_t i=0; i<10; i++){ |
117 |
|
|
hid.Form("htrkresx%i",i); |
118 |
|
|
htrkresx[i] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data())); |
119 |
|
|
hid.Form("htrkresy%i",i); |
120 |
|
|
htrkresy[i] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data())); |
121 |
|
|
} |
122 |
|
|
for(Int_t ip=0; ip<6; ip++){ |
123 |
|
|
hid.Form("hxcalotrk%i",ip); |
124 |
|
|
hxcalotrk[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data())); |
125 |
|
|
hid.Form("hycalotrk%i",ip); |
126 |
|
|
hycalotrk[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data())); |
127 |
|
|
} |
128 |
|
|
hresangx = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangx")); |
129 |
|
|
hresangy = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangy")); |
130 |
|
|
|
131 |
|
|
// ---------------------------------- |
132 |
|
|
// calo variables |
133 |
|
|
// ---------------------------------- |
134 |
|
|
CaloAxis *x_axis = new CaloAxis(); |
135 |
|
|
CaloAxis *y_axis = new CaloAxis(); |
136 |
|
|
|
137 |
|
|
float rcil = 1.;// tolerance (cm) |
138 |
|
|
x_axis->FitAxis(event->GetCaloLevel1(),0,rcil); |
139 |
|
|
y_axis->FitAxis(event->GetCaloLevel1(),1,rcil); |
140 |
|
|
|
141 |
|
|
float xs31 = x_axis->par[0]+x_axis->par[1]*event->GetZTOF(31); |
142 |
|
|
float xs32 = x_axis->par[0]+x_axis->par[1]*event->GetZTOF(32); |
143 |
|
|
float ys31 = y_axis->par[0]+y_axis->par[1]*event->GetZTOF(31); |
144 |
|
|
float ys32 = y_axis->par[0]+y_axis->par[1]*event->GetZTOF(32); |
145 |
|
|
float qtrack = x_axis->GetQaxis()+y_axis->GetQaxis(); |
146 |
|
|
int ntrack = x_axis->GetN()+y_axis->GetN(); |
147 |
|
|
|
148 |
|
|
CaloStrip top = CaloStrip(); |
149 |
|
|
top.Set(0,0,48); |
150 |
|
|
float zcalotop = top.GetZ(); |
151 |
|
|
|
152 |
|
|
float xcalotop = x_axis->par[0]+x_axis->par[1]*zcalotop; |
153 |
|
|
float ycalotop = y_axis->par[0]+y_axis->par[1]*zcalotop; |
154 |
|
|
|
155 |
|
|
|
156 |
|
|
float xtrk[6]; |
157 |
|
|
float ytrk[6]; |
158 |
|
|
for(Int_t ip=0; ip<6; ip++){ |
159 |
|
|
xtrk[ip] = x_axis->par[0]+x_axis->par[1]*event->GetZTrk(ip+1); |
160 |
|
|
ytrk[ip] = y_axis->par[0]+y_axis->par[1]*event->GetZTrk(ip+1); |
161 |
|
|
} |
162 |
|
|
|
163 |
|
|
float caloangx = (float)atan((double)x_axis->par[1])*180./3.1415026; |
164 |
|
|
float caloangy = (float)atan((double)y_axis->par[1])*180./3.1415026; |
165 |
|
|
|
166 |
|
|
float xtrktop = x_axis->par[0]+x_axis->par[1]*event->GetZTrk(1); |
167 |
|
|
float ytrktop = y_axis->par[0]+y_axis->par[1]*event->GetZTrk(1); |
168 |
|
|
|
169 |
|
|
x_axis->Delete(); |
170 |
|
|
y_axis->Delete(); |
171 |
|
|
|
172 |
|
|
// tolleranza sul contenimento della traccia nel tracker |
173 |
|
|
float dztrk = 44.51; |
174 |
|
|
float dxtrk = 16.14; |
175 |
|
|
float dytrk = 13.14; |
176 |
|
|
float axtoll = 10.; //deg |
177 |
|
|
float aytoll = 5.; //deg |
178 |
|
|
float tollx = tan(atan(dxtrk/dztrk)+axtoll*3.1415026/180.)*dztrk-dxtrk; |
179 |
|
|
float tolly = tan(atan(dytrk/dztrk)+aytoll*3.1415026/180.)*dztrk-dytrk; |
180 |
|
|
// cout << tollx << " "<<tolly << endl; |
181 |
|
|
|
182 |
|
|
// ---------------------------------- |
183 |
|
|
// tof variables |
184 |
|
|
// ---------------------------------- |
185 |
|
|
// get track stored by the tof |
186 |
|
|
ToFTrkVar *tof = event->GetToFStoredTrack(-1); |
187 |
|
|
if(!tof){ |
188 |
|
|
cout << " no ToF-track stored "<<endl; |
189 |
|
|
return; |
190 |
|
|
} |
191 |
|
|
float resxs21 = (event->GetZTOF(21)-event->GetZTOF(32))*(tof->xtofpos[0]-tof->xtofpos[2])/(event->GetZTOF(12)-event->GetZTOF(32))+tof->xtofpos[2]; |
192 |
|
|
float resys22 = (event->GetZTOF(22)-event->GetZTOF(31))*(tof->ytofpos[0]-tof->ytofpos[2])/(event->GetZTOF(11)-event->GetZTOF(31))+tof->ytofpos[2]; |
193 |
|
|
resxs21 -= tof->xtofpos[1]; |
194 |
|
|
resys22 -= tof->ytofpos[1]; |
195 |
|
|
|
196 |
|
|
float resxs32 = xs32 - tof->xtofpos[2]; |
197 |
|
|
float resys31 = ys31 - tof->ytofpos[2]; |
198 |
|
|
|
199 |
|
|
// ---------------------------------- |
200 |
|
|
// tracker variables |
201 |
|
|
// ---------------------------------- |
202 |
|
|
float resxtrk[]={-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.}; |
203 |
|
|
float resytrk[]={-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.,-999.}; |
204 |
|
|
|
205 |
|
|
float resangx = -999; |
206 |
|
|
float resangy = -999; |
207 |
|
|
if( event->GetNTracks()==1 ){ |
208 |
|
|
|
209 |
|
|
PamTrack *track = event->GetTrack(0); |
210 |
|
|
if( track->chi2 >0 ){ |
211 |
|
|
|
212 |
|
|
Float_t ztraj[13]; |
213 |
|
|
|
214 |
|
|
Int_t i=0; |
215 |
|
|
ztraj[i++] = event->GetZTOF(11); |
216 |
|
|
ztraj[i++] = event->GetZTOF(12); |
217 |
|
|
ztraj[i++] = event->GetZTOF(21); |
218 |
|
|
ztraj[i++] = event->GetZTOF(22); |
219 |
|
|
ztraj[i++] = event->GetZTrk(1); |
220 |
|
|
ztraj[i++] = event->GetZTrk(2); |
221 |
|
|
ztraj[i++] = event->GetZTrk(3); |
222 |
|
|
ztraj[i++] = event->GetZTrk(4); |
223 |
|
|
ztraj[i++] = event->GetZTrk(5); |
224 |
|
|
ztraj[i++] = event->GetZTrk(6); |
225 |
|
|
ztraj[i++] = event->GetZTOF(31); |
226 |
|
|
ztraj[i++] = event->GetZTOF(32); |
227 |
|
|
ztraj[i++] = zcalotop; |
228 |
|
|
Trajectory traj = Trajectory(13,ztraj); |
229 |
|
|
traj.DoTrack2( track->al ); |
230 |
|
|
|
231 |
|
|
// --- S1 |
232 |
|
|
resxtrk[0] = track->xtofpos[0] - traj.x[1]; |
233 |
|
|
resytrk[0] = track->ytofpos[0] - traj.y[0]; |
234 |
|
|
// --- S2 |
235 |
|
|
resxtrk[1] = track->xtofpos[1] - traj.x[2]; |
236 |
|
|
resytrk[1] = track->ytofpos[1] - traj.y[3]; |
237 |
|
|
// --- S3 |
238 |
|
|
resxtrk[2] = track->xtofpos[2] - traj.x[11]; |
239 |
|
|
resytrk[2] = track->ytofpos[2] - traj.y[10]; |
240 |
|
|
// --- trk |
241 |
|
|
for(Int_t ip=0; ip<6; ip++){ |
242 |
|
|
if(track->XGood(ip))resxtrk[3+ip] = track->xm[ip] - traj.x[4+ip]; |
243 |
|
|
if(track->YGood(ip))resytrk[3+ip] = track->ym[ip] - traj.y[4+ip]; |
244 |
|
|
} |
245 |
|
|
// --- calo |
246 |
|
|
resxtrk[9] = xcalotop - traj.x[12]; |
247 |
|
|
resytrk[9] = ycalotop - traj.y[12]; |
248 |
|
|
|
249 |
|
|
resangx = caloangx - traj.thx[12]; |
250 |
|
|
resangy = caloangy - traj.thy[12]; |
251 |
|
|
// cout << caloangx << " -x- " << traj.thx[12]<<endl; |
252 |
|
|
// cout << caloangy << " -y- " << traj.thy[12]<<endl; |
253 |
|
|
|
254 |
|
|
} |
255 |
|
|
if(track) delete track; |
256 |
|
|
} |
257 |
|
|
|
258 |
|
|
// ---------------------------------- |
259 |
|
|
// fill spatial residuals |
260 |
|
|
// ---------------------------------- |
261 |
|
|
htofresx->Fill(resxs21); |
262 |
|
|
htofresy->Fill(resys22); |
263 |
|
|
hcalresx->Fill(resxs32); |
264 |
|
|
hcalresy->Fill(resys31); |
265 |
|
|
for(Int_t ih=0; ih<10; ih++){ |
266 |
|
|
htrkresx[ih]->Fill(resxtrk[ih]); |
267 |
|
|
htrkresy[ih]->Fill(resytrk[ih]); |
268 |
|
|
} |
269 |
|
|
for(Int_t ip=0; ip<6; ip++){ |
270 |
|
|
hxcalotrk[ip]->Fill(xtrk[ip]); |
271 |
|
|
hycalotrk[ip]->Fill(ytrk[ip]); |
272 |
|
|
} |
273 |
|
|
hresangx->Fill(resangx); |
274 |
|
|
hresangy->Fill(resangy); |
275 |
|
|
|
276 |
|
|
|
277 |
|
|
} |
278 |
|
|
//=============================================================== |
279 |
|
|
// Save histograms |
280 |
|
|
//=============================================================== |
281 |
|
|
void SaveHistos( TFile * outf ){ |
282 |
|
|
|
283 |
|
|
TString hid; |
284 |
|
|
|
285 |
|
|
gROOT->cd(); |
286 |
|
|
// ---------------------------------- |
287 |
|
|
// retrieve histos |
288 |
|
|
// ---------------------------------- |
289 |
|
|
htofresx = dynamic_cast<TH1F*>(gDirectory->FindObject("htofresx")); |
290 |
|
|
htofresy = dynamic_cast<TH1F*>(gDirectory->FindObject("htofresy")); |
291 |
|
|
// cout << hex << htofresx << dec << endl; |
292 |
|
|
hcalresx = dynamic_cast<TH1F*>(gDirectory->FindObject("hcalresx")); |
293 |
|
|
hcalresy = dynamic_cast<TH1F*>(gDirectory->FindObject("hcalresy")); |
294 |
|
|
for(Int_t i=0; i<10; i++){ |
295 |
|
|
hid.Form("htrkresx%i",i); |
296 |
|
|
htrkresx[i] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data())); |
297 |
|
|
hid.Form("htrkresy%i",i); |
298 |
|
|
htrkresy[i] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data())); |
299 |
|
|
} |
300 |
|
|
for(Int_t ip=0; ip<6; ip++){ |
301 |
|
|
hid.Form("hxcalotrk%i",ip); |
302 |
|
|
hxcalotrk[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data())); |
303 |
|
|
hid.Form("hycalotrk%i",ip); |
304 |
|
|
hycalotrk[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data())); |
305 |
|
|
} |
306 |
|
|
hresangx = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangx")); |
307 |
|
|
hresangy = dynamic_cast<TH1F*>(gDirectory->FindObject("hresangy")); |
308 |
|
|
|
309 |
|
|
// ---------------------------------- |
310 |
|
|
// save on file |
311 |
|
|
// ---------------------------------- |
312 |
|
|
|
313 |
|
|
outf->cd(); |
314 |
|
|
|
315 |
|
|
if(htofresx)htofresx->Write(); |
316 |
|
|
if(htofresy)htofresy->Write(); |
317 |
|
|
if(hcalresx)hcalresx->Write(); |
318 |
|
|
if(hcalresy)hcalresy->Write(); |
319 |
|
|
|
320 |
|
|
for(Int_t i=0; i<10; i++){ |
321 |
|
|
if(htrkresx[i])htrkresx[i]->Write(); |
322 |
|
|
if(htrkresy[i])htrkresy[i]->Write(); |
323 |
|
|
} |
324 |
|
|
|
325 |
|
|
for(Int_t ip=0; ip<6; ip++){ |
326 |
|
|
if(hxcalotrk[ip])hxcalotrk[ip]->Write(); |
327 |
|
|
if(hycalotrk[ip])hycalotrk[ip]->Write(); |
328 |
|
|
} |
329 |
|
|
|
330 |
|
|
if(hresangx) hresangx->Write(); |
331 |
|
|
if(hresangy) hresangy->Write(); |
332 |
|
|
} |
333 |
|
|
|