1 |
mocchiut |
1.1 |
#include <CaloProfile.h> |
2 |
|
|
|
3 |
|
|
//////////////////////////////////////////////////////////////////////// |
4 |
|
|
/** |
5 |
|
|
* 1-dimension function describing lateral distribution of the |
6 |
|
|
* shower as viewed by calorimeter |
7 |
|
|
* (projection of 3d function in one direction) |
8 |
|
|
* |
9 |
|
|
* xi[0] = x or y coordinate relative to shower axis |
10 |
|
|
* parmin[0] = rt |
11 |
|
|
* parmin[1] = p |
12 |
|
|
* parmin[2] = rc |
13 |
|
|
* |
14 |
|
|
*/ |
15 |
|
|
//////////////////////////////////////////////////////////////////////// |
16 |
|
|
Double_t cfradx(Double_t *xi, Double_t *parmin) { |
17 |
|
|
|
18 |
|
|
double fradxmin2,p,rt,rc,es,x,pig,norm,c; |
19 |
|
|
x=*xi; |
20 |
|
|
pig = acos(-1.); |
21 |
|
|
rt=parmin[0]; |
22 |
|
|
p=parmin[1]; |
23 |
|
|
rc=parmin[2]; |
24 |
|
|
norm=parmin[3]; |
25 |
|
|
c=parmin[4]; |
26 |
|
|
x=x-c; |
27 |
|
|
es=1.5; |
28 |
|
|
fradxmin2=p*pig*pow(rc,2)/pow((pow(x,2)+pow(rc,2)),es); |
29 |
|
|
fradxmin2=fradxmin2+(1-p)*pig*pow(rt,2)/pow((pow(x,2)+pow(rt,2)),es); |
30 |
|
|
fradxmin2=norm*fradxmin2/(2*pig); |
31 |
|
|
//cout<<"x,fradxmin2 "<< x<<" "<<fradxmin2 <<endl; |
32 |
|
|
return fradxmin2; |
33 |
|
|
} |
34 |
|
|
|
35 |
|
|
|
36 |
|
|
Double_t cfunc(Double_t *tin,Double_t *par) |
37 |
|
|
{ |
38 |
|
|
|
39 |
|
|
Double_t gaa,a,tm,et,value,t,t0; |
40 |
|
|
t=*tin; |
41 |
|
|
et=par[0]; |
42 |
|
|
a=par[1]; |
43 |
|
|
tm=par[2]; |
44 |
|
|
t0=par[3]; |
45 |
|
|
gaa=TMath::Gamma(a); |
46 |
|
|
|
47 |
|
|
value=et*((a-1)/(tm*gaa))*pow(((a-1)*(t-t0)/tm),(a-1))*exp(-(a-1)*(t-t0)/tm); |
48 |
|
|
return value; |
49 |
|
|
} |
50 |
|
|
|
51 |
|
|
|
52 |
|
|
//-------------------------------------- |
53 |
|
|
/** |
54 |
|
|
* Default constructor |
55 |
|
|
*/ |
56 |
|
|
CaloLat::CaloLat(){ |
57 |
|
|
Clear(); |
58 |
|
|
}; |
59 |
|
|
|
60 |
|
|
/** |
61 |
|
|
* Default constructor |
62 |
|
|
*/ |
63 |
|
|
CaloLong::CaloLong(){ |
64 |
|
|
Clear(); |
65 |
|
|
}; |
66 |
|
|
|
67 |
|
|
CaloLat::CaloLat(PamLevel2 *l2p){ |
68 |
|
|
// |
69 |
|
|
Clear(); |
70 |
|
|
// |
71 |
|
|
L2 = l2p; |
72 |
|
|
// |
73 |
|
|
if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n"); |
74 |
|
|
// |
75 |
|
|
OBT = 0; |
76 |
|
|
PKT = 0; |
77 |
|
|
atime = 0; |
78 |
|
|
// |
79 |
|
|
debug = false; |
80 |
|
|
// |
81 |
|
|
}; |
82 |
|
|
|
83 |
|
|
CaloLong::CaloLong(PamLevel2 *l2p){ |
84 |
|
|
// |
85 |
|
|
Clear(); |
86 |
|
|
// |
87 |
|
|
L2 = l2p; |
88 |
|
|
// |
89 |
|
|
if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n"); |
90 |
|
|
// |
91 |
|
|
OBT = 0; |
92 |
|
|
PKT = 0; |
93 |
|
|
atime = 0; |
94 |
|
|
// |
95 |
|
|
debug = false; |
96 |
|
|
// |
97 |
|
|
}; |
98 |
|
|
|
99 |
|
|
void CaloLat::Clear(){ |
100 |
|
|
// |
101 |
|
|
}; |
102 |
|
|
|
103 |
|
|
void CaloLong::Clear(){ |
104 |
|
|
// |
105 |
|
|
}; |
106 |
|
|
|
107 |
|
|
void CaloLat::Print(){ |
108 |
|
|
// |
109 |
|
|
Process(); |
110 |
|
|
// |
111 |
|
|
printf("==================== Calorimeter Lateral Profile =======================\n"); |
112 |
|
|
printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime); |
113 |
|
|
// printf(" nx [number of X combination]:.. %i\n",nx); |
114 |
|
|
// printf(" ny [number of Y combination]:.. %i\n",ny); |
115 |
|
|
printf("========================================================================\n"); |
116 |
|
|
// |
117 |
|
|
}; |
118 |
|
|
|
119 |
|
|
void CaloLong::Print(){ |
120 |
|
|
// |
121 |
|
|
Process(); |
122 |
|
|
// |
123 |
|
|
printf("==================== Calorimeter Lateral Profile =======================\n"); |
124 |
|
|
printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime); |
125 |
|
|
// printf(" nx [number of X combination]:.. %i\n",nx); |
126 |
|
|
// printf(" ny [number of Y combination]:.. %i\n",ny); |
127 |
|
|
printf("========================================================================\n"); |
128 |
|
|
// |
129 |
|
|
}; |
130 |
|
|
|
131 |
|
|
void CaloLat::Draw(){ |
132 |
|
|
// |
133 |
|
|
Process(); |
134 |
|
|
Draw(-1,-1); |
135 |
|
|
}; |
136 |
|
|
|
137 |
|
|
void CaloLong::Draw(){ |
138 |
|
|
// |
139 |
|
|
Process(); |
140 |
|
|
Draw(-1); |
141 |
|
|
}; |
142 |
|
|
|
143 |
|
|
|
144 |
|
|
void CaloLat::Draw(Int_t view,Int_t plane){ |
145 |
|
|
// |
146 |
|
|
Int_t minv = 0; |
147 |
|
|
Int_t maxv = 0; |
148 |
|
|
Int_t minp = 0; |
149 |
|
|
Int_t maxp = 0; |
150 |
|
|
// |
151 |
|
|
if ( view == -1 ){ |
152 |
|
|
minv = 0; |
153 |
|
|
maxv = 2; |
154 |
|
|
} else { |
155 |
|
|
minv = view; |
156 |
|
|
maxv = view+1; |
157 |
|
|
}; |
158 |
|
|
|
159 |
|
|
if ( plane == -1 ){ |
160 |
|
|
minp = 0; |
161 |
|
|
maxp = 22; |
162 |
|
|
} else { |
163 |
|
|
minp = plane; |
164 |
|
|
maxp = plane+1; |
165 |
|
|
}; |
166 |
|
|
// |
167 |
|
|
Process(); |
168 |
|
|
// |
169 |
|
|
gStyle->SetLabelSize(0.04); |
170 |
|
|
gStyle->SetNdivisions(510,"XY"); |
171 |
|
|
// |
172 |
|
|
for (Int_t v=minv; v<maxv;v++){ |
173 |
|
|
for (Int_t p=minp; p<maxp;p++){ |
174 |
|
|
TString hid = Form("clatv%ip%i",v,p); |
175 |
|
|
TCanvas *tc = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid)); |
176 |
|
|
if ( tc ){ |
177 |
|
|
// tc->Clear(); |
178 |
|
|
} else { |
179 |
|
|
tc = new TCanvas(hid,hid); |
180 |
|
|
}; |
181 |
|
|
// |
182 |
|
|
TString thid = Form("hlatv%ip%i",v,p); |
183 |
|
|
TH1F *th = dynamic_cast<TH1F*>(gDirectory->FindObject(thid)); |
184 |
|
|
if ( th ) th->Delete(); |
185 |
|
|
// th->Clear(); |
186 |
|
|
// th->Reset(); |
187 |
|
|
// } else { |
188 |
|
|
th = new TH1F(thid,thid,96,-0.5,95.5); |
189 |
|
|
// }; |
190 |
|
|
tc->cd(); |
191 |
|
|
// |
192 |
|
|
for (Int_t st=0;st<96;st++){ |
193 |
|
|
th->Fill(st,estrip[v][p][st]); |
194 |
|
|
}; |
195 |
|
|
th->Draw(); |
196 |
|
|
tc->Modified(); |
197 |
|
|
tc->Update(); |
198 |
|
|
}; |
199 |
|
|
}; |
200 |
|
|
// |
201 |
|
|
gStyle->SetLabelSize(0); |
202 |
|
|
gStyle->SetNdivisions(1,"XY"); |
203 |
|
|
// |
204 |
|
|
}; |
205 |
|
|
|
206 |
|
|
void CaloLong::Draw(Int_t view){ |
207 |
|
|
// |
208 |
|
|
Int_t minv = 0; |
209 |
|
|
Int_t maxv = 0; |
210 |
|
|
// |
211 |
|
|
if ( view == -1 ){ |
212 |
|
|
maxv = -1; |
213 |
|
|
} else { |
214 |
|
|
minv = view; |
215 |
|
|
maxv = view+1; |
216 |
|
|
}; |
217 |
|
|
// |
218 |
|
|
Process(); |
219 |
|
|
// |
220 |
|
|
gStyle->SetLabelSize(0.04); |
221 |
|
|
gStyle->SetNdivisions(510,"XY"); |
222 |
|
|
// |
223 |
|
|
if ( maxv != -1 ){ |
224 |
|
|
for (Int_t v=minv; v<maxv;v++){ |
225 |
|
|
TString hid = Form("clongv%i",v); |
226 |
|
|
TCanvas *tc = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid)); |
227 |
|
|
if ( tc ){ |
228 |
|
|
// tc->Clear(); |
229 |
|
|
} else { |
230 |
|
|
tc = new TCanvas(hid,hid); |
231 |
|
|
}; |
232 |
|
|
// |
233 |
|
|
TString thid = Form("hlongv%i",v); |
234 |
|
|
TH1F *th = dynamic_cast<TH1F*>(gDirectory->FindObject(thid)); |
235 |
|
|
if ( th ) th->Delete(); |
236 |
|
|
// th->Clear(); |
237 |
|
|
// th->Reset(); |
238 |
|
|
// } else { |
239 |
|
|
th = new TH1F(thid,thid,22,-0.5,21.5); |
240 |
|
|
// }; |
241 |
|
|
tc->cd(); |
242 |
|
|
// |
243 |
|
|
for (Int_t st=0;st<22;st++){ |
244 |
|
|
th->Fill(st,eplane[v][st]); |
245 |
|
|
}; |
246 |
|
|
th->Draw(); |
247 |
|
|
tc->Modified(); |
248 |
|
|
tc->Update(); |
249 |
|
|
}; |
250 |
|
|
} else { |
251 |
|
|
// |
252 |
|
|
TString hid = Form("clongvyvx"); |
253 |
|
|
TCanvas *tc = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid)); |
254 |
|
|
if ( tc ){ |
255 |
|
|
} else { |
256 |
|
|
tc = new TCanvas(hid,hid); |
257 |
|
|
}; |
258 |
|
|
// |
259 |
|
|
TString thid = Form("hlongvyvx"); |
260 |
|
|
TH1F *th = dynamic_cast<TH1F*>(gDirectory->FindObject(thid)); |
261 |
|
|
if ( th ) th->Delete(); |
262 |
|
|
th = new TH1F(thid,thid,44,-0.5,43.5); |
263 |
|
|
tc->cd(); |
264 |
|
|
Int_t pp=0; |
265 |
|
|
for (Int_t st=0;st<22;st++){ |
266 |
|
|
for (Int_t v=1; v>=0;v--){ |
267 |
|
|
// |
268 |
|
|
th->Fill(pp,eplane[v][st]); |
269 |
|
|
// |
270 |
|
|
pp++; |
271 |
|
|
}; |
272 |
|
|
}; |
273 |
|
|
th->Draw(); |
274 |
|
|
tc->Modified(); |
275 |
|
|
tc->Update(); |
276 |
|
|
}; |
277 |
|
|
// |
278 |
|
|
gStyle->SetLabelSize(0); |
279 |
|
|
gStyle->SetNdivisions(1,"XY"); |
280 |
|
|
// |
281 |
|
|
}; |
282 |
|
|
|
283 |
|
|
void CaloLat::Delete(){ |
284 |
|
|
Clear(); |
285 |
|
|
//delete this; |
286 |
|
|
}; |
287 |
|
|
|
288 |
|
|
void CaloLong::Delete(){ |
289 |
|
|
Clear(); |
290 |
|
|
//delete this; |
291 |
|
|
}; |
292 |
|
|
|
293 |
|
|
void CaloLat::Process(){ |
294 |
|
|
// |
295 |
|
|
if ( !L2 ){ |
296 |
|
|
printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n"); |
297 |
|
|
printf(" ERROR: CaloHough variables not filled \n"); |
298 |
|
|
return; |
299 |
|
|
}; |
300 |
|
|
// |
301 |
|
|
Bool_t newentry = false; |
302 |
|
|
// |
303 |
|
|
if ( L2->IsORB() ){ |
304 |
|
|
if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){ |
305 |
|
|
newentry = true; |
306 |
|
|
OBT = L2->GetOrbitalInfo()->OBT; |
307 |
|
|
PKT = L2->GetOrbitalInfo()->pkt_num; |
308 |
|
|
atime = L2->GetOrbitalInfo()->absTime; |
309 |
|
|
}; |
310 |
|
|
} else { |
311 |
|
|
newentry = true; |
312 |
|
|
}; |
313 |
|
|
// |
314 |
|
|
if ( !newentry ) return; |
315 |
|
|
// |
316 |
|
|
if ( debug ) printf(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime); |
317 |
|
|
// |
318 |
|
|
Clear(); |
319 |
|
|
// |
320 |
|
|
// let's start |
321 |
|
|
// |
322 |
|
|
memset(estrip,0, 4224*sizeof(Float_t)); |
323 |
|
|
Float_t mip1 = 0.; |
324 |
|
|
Int_t view1 = 0; |
325 |
|
|
Int_t plane1 = 0; |
326 |
|
|
Int_t strip1 = 0; |
327 |
|
|
// |
328 |
|
|
for (Int_t i=0; i<L2->GetCaloLevel1()->istrip ; i++){ |
329 |
|
|
mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1); |
330 |
|
|
estrip[view1][plane1][strip1] = mip1; |
331 |
|
|
}; |
332 |
|
|
// |
333 |
|
|
if ( debug ) this->Print(); |
334 |
|
|
if ( debug ) printf(" exit \n"); |
335 |
|
|
// |
336 |
|
|
}; |
337 |
|
|
|
338 |
|
|
void CaloLong::Process(){ |
339 |
|
|
// |
340 |
|
|
if ( !L2 ){ |
341 |
|
|
printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n"); |
342 |
|
|
printf(" ERROR: CaloHough variables not filled \n"); |
343 |
|
|
return; |
344 |
|
|
}; |
345 |
|
|
// |
346 |
|
|
Bool_t newentry = false; |
347 |
|
|
// |
348 |
|
|
if ( L2->IsORB() ){ |
349 |
|
|
if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){ |
350 |
|
|
newentry = true; |
351 |
|
|
OBT = L2->GetOrbitalInfo()->OBT; |
352 |
|
|
PKT = L2->GetOrbitalInfo()->pkt_num; |
353 |
|
|
atime = L2->GetOrbitalInfo()->absTime; |
354 |
|
|
}; |
355 |
|
|
} else { |
356 |
|
|
newentry = true; |
357 |
|
|
}; |
358 |
|
|
// |
359 |
|
|
if ( !newentry ) return; |
360 |
|
|
// |
361 |
|
|
if ( debug ) printf(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime); |
362 |
|
|
// |
363 |
|
|
Clear(); |
364 |
|
|
// |
365 |
|
|
// let's start |
366 |
|
|
// |
367 |
|
|
memset(eplane,0, 2*22*sizeof(Float_t)); |
368 |
|
|
// |
369 |
|
|
for (Int_t v=0; v<2; v++){ |
370 |
|
|
for (Int_t i=0; i<22; i++){ |
371 |
|
|
eplane[v][i] = L2->GetCaloLevel1()->qtotpl(v,i); |
372 |
|
|
}; |
373 |
|
|
}; |
374 |
|
|
// |
375 |
|
|
if ( debug ) this->Print(); |
376 |
|
|
if ( debug ) printf(" exit \n"); |
377 |
|
|
// |
378 |
|
|
}; |
379 |
|
|
|