1 |
#include <CaloBragg.h> |
2 |
|
3 |
|
4 |
ClassImp(CaloBragg); |
5 |
//-------------------------------------- |
6 |
/* |
7 |
* Default constructor |
8 |
*/ |
9 |
CaloBragg::CaloBragg(){ |
10 |
Clear(); |
11 |
}; |
12 |
|
13 |
CaloBragg::CaloBragg(PamLevel2 *l2p){ |
14 |
// |
15 |
Clear(); |
16 |
LoadParam(); |
17 |
// |
18 |
L2 = l2p; |
19 |
// |
20 |
if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n"); |
21 |
// |
22 |
OBT = 0; |
23 |
PKT = 0; |
24 |
atime = 0; |
25 |
// |
26 |
debug = false; |
27 |
usetrack = false; |
28 |
usepl18x = false; |
29 |
newchi2 = false; |
30 |
usenewBB = false; |
31 |
fzeta = -1.; |
32 |
// |
33 |
}; |
34 |
|
35 |
void CaloBragg::Clear(){ |
36 |
// |
37 |
ndf = 0; |
38 |
tr = 0; |
39 |
sntr = 0; |
40 |
// qtchi2 = 0.; |
41 |
// qtz = 0.; |
42 |
// qtetot = 0.; |
43 |
// qtpskip = 0.; |
44 |
lpchi2 = 0.; |
45 |
lpz = 0.; |
46 |
lpetot = 0.; |
47 |
lppskip = 0.; |
48 |
memset(calorimetro,0,44*2*sizeof(Float_t)); |
49 |
memset(spessore,0,4*sizeof(Float_t)); |
50 |
memset(estremi,0,2*2*sizeof(Float_t)); |
51 |
Integrale=0.; |
52 |
|
53 |
for(Int_t l=0;l<44;l++){ |
54 |
calorimetro[l][0]=-1.; |
55 |
} |
56 |
|
57 |
}; |
58 |
|
59 |
void CaloBragg::Print(){ |
60 |
// |
61 |
|
62 |
if(!debug) Process(); |
63 |
// |
64 |
printf("========================================================================\n"); |
65 |
printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack); |
66 |
printf(" first plane: %f \n", estremi[0][0]); |
67 |
printf(" last plane: %f \n", estremi[1][0]); |
68 |
// printf(" chi 2 from truncated mean: %f \n", qtchi2); |
69 |
// printf(" Z from truncated mean %f: \n", qtz); |
70 |
// printf(" energy from truncated mean %f: \n", qtetot); |
71 |
// printf(" plane not used for truncated mean %f: \n", qtpskip); |
72 |
printf(" chi 2 from loop %f: \n", lpchi2); |
73 |
printf(" Z from loop %f: \n", lpz); |
74 |
printf(" energy from loop %f: \n", lpetot); |
75 |
printf(" plane not used for loop %f: \n", lppskip); |
76 |
printf(" ndf: %i \n",ndf); |
77 |
printf("========================================================================\n"); |
78 |
// |
79 |
}; |
80 |
|
81 |
void CaloBragg::Delete(){ |
82 |
Clear(); |
83 |
//delete this; |
84 |
}; |
85 |
|
86 |
|
87 |
void CaloBragg::Process(){ |
88 |
Process(-1); |
89 |
}; |
90 |
|
91 |
|
92 |
void CaloBragg::CleanPlanes(Float_t epiano[22][2]){ |
93 |
// return; |
94 |
Int_t hitplanes = 0; |
95 |
for (Int_t i = 0; i<22; i++){ |
96 |
for (Int_t j = 1; j>=0; j--){ |
97 |
if ( epiano[i][j] > 0.7 ) hitplanes++; |
98 |
}; |
99 |
}; |
100 |
Float_t lowlim = 0.85; |
101 |
Float_t dedxone = 0.; |
102 |
Float_t step1 = 0.8*L2->GetCaloLevel2()->qtot/(Float_t)hitplanes; |
103 |
while ( dedxone < step1 ){ |
104 |
for (Int_t i = 0; i<22; i++){ |
105 |
for (Int_t j = 1; j>=0; j--){ |
106 |
if ( epiano[i][j] >= step1 && dedxone < 0.7 ) dedxone = epiano[i][j]; |
107 |
}; |
108 |
}; |
109 |
} |
110 |
if ( dedxone < 0.7 ){ |
111 |
for (Int_t i = 0; i<22; i++){ |
112 |
for (Int_t j = 1; j>=0; j--){ |
113 |
if ( epiano[i][j] > 0. && dedxone < 0.7 ) dedxone = epiano[i][j]; |
114 |
}; |
115 |
}; |
116 |
} |
117 |
// |
118 |
// printf(" dedxone = %f step1 %f \n",dedxone,step1); |
119 |
Bool_t revulsera = false; |
120 |
Bool_t nullius = false; |
121 |
Int_t nulliferus = 0; |
122 |
for (Int_t i = 0; i<22; i++){ |
123 |
for (Int_t j = 1; j>=0; j--){ |
124 |
if ( epiano[i][j] < dedxone*lowlim ){ |
125 |
// printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius); |
126 |
epiano[i][j] = 0.; |
127 |
} else { |
128 |
//x printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius); |
129 |
nulliferus = 0; |
130 |
revulsera = true; |
131 |
}; |
132 |
if ( epiano[i][j] < 0.7 && revulsera ) nulliferus++; |
133 |
if ( nulliferus > 10 ) nullius = true; |
134 |
if ( nullius ) epiano[i][j] = 0.; |
135 |
}; |
136 |
}; |
137 |
|
138 |
} |
139 |
|
140 |
|
141 |
void CaloBragg::Process(Int_t ntr){ |
142 |
// |
143 |
if ( !L2 ){ |
144 |
printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n"); |
145 |
printf(" ERROR: CaloBragg variables not filled \n"); |
146 |
return; |
147 |
}; |
148 |
// |
149 |
Bool_t newentry = false; |
150 |
// |
151 |
if ( L2->IsORB() ){ |
152 |
if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || ntr != sntr ){ |
153 |
newentry = true; |
154 |
OBT = L2->GetOrbitalInfo()->OBT; |
155 |
PKT = L2->GetOrbitalInfo()->pkt_num; |
156 |
atime = L2->GetOrbitalInfo()->absTime; |
157 |
sntr = ntr; |
158 |
}; |
159 |
} else { |
160 |
newentry = true; |
161 |
}; |
162 |
// |
163 |
if ( !newentry ) return; |
164 |
// |
165 |
tr = ntr; |
166 |
// |
167 |
if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime); |
168 |
// |
169 |
Clear(); |
170 |
|
171 |
// |
172 |
// |
173 |
// |
174 |
Int_t view = 0; |
175 |
Int_t plane = 0; |
176 |
Int_t strip = 0; |
177 |
Float_t mip = 0.; |
178 |
Float_t epiano[22][2]; |
179 |
memset(epiano,0,22*2*sizeof(Float_t)); |
180 |
for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){ |
181 |
// |
182 |
mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); |
183 |
// |
184 |
if ( !usepl18x && view==0 && plane==18 ) mip = 0.; |
185 |
// |
186 |
epiano[plane][view]+=mip; |
187 |
// |
188 |
// |
189 |
}; |
190 |
// |
191 |
this->CleanPlanes(*&epiano); |
192 |
// |
193 |
PamTrack *ptrack = 0; |
194 |
CaloTrkVar *track = 0; |
195 |
// |
196 |
if ( usetrack ){ |
197 |
if ( ntr >= 0 ){ |
198 |
ptrack = L2->GetTrack(ntr); |
199 |
if ( ptrack ) track = ptrack->GetCaloTrack(); |
200 |
} else { |
201 |
track = L2->GetCaloStoredTrack(ntr); //al momento e' vera solo questa riga |
202 |
}; |
203 |
// |
204 |
if ( !track && ntr >= 0 ){ |
205 |
printf(" ERROR: cannot find any track!\n"); |
206 |
printf(" ERROR: CaloBragg variables not completely filled \n"); |
207 |
return; |
208 |
}; |
209 |
} else { |
210 |
if ( ntr >= 0 ){ |
211 |
if ( debug ) printf(" ERROR: you asked not to use a track but you are looking for track number %i !\n",ntr); |
212 |
if ( debug ) printf(" ERROR: CaloBragg variables not completely filled \n"); |
213 |
return; |
214 |
}; |
215 |
}; |
216 |
// |
217 |
if(L2->GetCaloLevel2()->npcfit[0]==0 && L2->GetCaloLevel2()->npcfit[1]==0 && L2->GetCaloLevel2()->npcfit[2]==0 && L2->GetCaloLevel2()->npcfit[3]==0) return;// controllo sulla traccia nel calorimetro |
218 |
|
219 |
// |
220 |
for(Int_t p=0; p<22; p++){ |
221 |
for(Int_t v=0; v<2; v++){ |
222 |
/*per usare traccia non del calo camboare cibar*/ |
223 |
calorimetro[(2*p)+1-v][0] = L2->GetCaloLevel2()->cibar[p][v];//strip attraversata |
224 |
calorimetro[(2*p)+1-v][1] = epiano[p][v]; //energia del piano //(epiano[p][v])/0.89 |
225 |
}; |
226 |
}; |
227 |
|
228 |
/*per ogni evento calcolo la conversione mip e w attraversato in equivalente Si*/ |
229 |
conversione(); // out: 1) g/cm2 Si , 2) spessoreW equivalente in Si, 3)Mip corretta per inclinazione |
230 |
|
231 |
/*settaggio della soglia per il loop sulla determinazione del piano di partenza */ |
232 |
Float_t ordplane[44];//mi serve per la media troncata |
233 |
memset(ordplane,0,44*sizeof(Float_t)); |
234 |
|
235 |
for(Int_t ipla=0; ipla< 2*NPLA; ipla++) ordplane[ipla]=calorimetro[ipla][1]; //energia del piano |
236 |
|
237 |
|
238 |
//ordino tutte le energie dei piani in ordine crescente |
239 |
|
240 |
Long64_t work[200]; |
241 |
Int_t ind = 0; |
242 |
//Int_t l = 0; |
243 |
Int_t RN = 0; |
244 |
Float_t sum4 = 0.; |
245 |
Float_t qm = 0.; |
246 |
while ( RN < 4 && ind < 44 ){ |
247 |
qm = TMath::KOrdStat((Long64_t)44,ordplane,(Long64_t)ind,work); |
248 |
if (qm >= 0.7 ){ |
249 |
if ( RN < 4 ){ |
250 |
sum4 += qm; |
251 |
RN++; |
252 |
}; |
253 |
}; |
254 |
ind++; |
255 |
}; |
256 |
// |
257 |
//sum4 /= (Float_t)RN; |
258 |
Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2]))); |
259 |
if(Zmean ==0.) Zmean=1.; |
260 |
if ( Zmean < 1. ) Zmean = 1.; |
261 |
|
262 |
|
263 |
/*trova primo e ultimo piano attraversati*/ |
264 |
Int_t p = 0;//contatore piani |
265 |
//per il primo parte da 0 e va in giu' |
266 |
while( estremi[0][1] <= 0. && p<(2*NPLA) ){ // era ==0 ma ricorda i problemi con Float == !!!!! |
267 |
// if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){ |
268 |
// if( (calorimetro[p][0] >0) && (calorimetro[p][1]*MIP >0.3)){ //0.7 mip = 70MeV soglia minima |
269 |
if( (calorimetro[p][0] >0) && (calorimetro[p][1]*MIP >Zmean*0.7)){ // 70% della MIP |
270 |
estremi[0][0]=p; |
271 |
estremi[0][1]=calorimetro[p][1] *MIP; //energia in MeV |
272 |
}; |
273 |
p++; |
274 |
}; |
275 |
|
276 |
//ultimo parte da 44 e sale |
277 |
p=43; |
278 |
while( (estremi[1][1] <= 0.) && (p>(int)estremi[0][0]) ){ |
279 |
if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){ |
280 |
estremi[1][0]=p;// |
281 |
estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV |
282 |
}; |
283 |
p = p-1; |
284 |
}; |
285 |
// |
286 |
|
287 |
Float_t lastok = 0.; |
288 |
// Bool_t goback = false; |
289 |
for ( int o = 0; o < estremi[1][0]; o++ ){ |
290 |
// |
291 |
if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1]; |
292 |
if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok; |
293 |
// if ( calorimetro[o][1] < 0.7 ) goback = true; |
294 |
// |
295 |
}; |
296 |
lastok = 0.; |
297 |
// if ( goback ){ |
298 |
for ( int o = estremi[1][0]; o >= 0; o-- ){ |
299 |
// |
300 |
// printf(" goback1: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok); |
301 |
if ( o < estremi[1][0] && calorimetro[o][1] > calorimetro[o+1][1]*1.2 && lastok > 0. ) calorimetro[o][1] = lastok; |
302 |
if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1]; |
303 |
if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok; |
304 |
// printf(" goback2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok); |
305 |
// |
306 |
}; |
307 |
// }; |
308 |
|
309 |
if ( startZero ) { |
310 |
estremi[0][0] = 0.; |
311 |
// estremi[0][1] = 0.; |
312 |
} |
313 |
|
314 |
/*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/ |
315 |
for(Int_t pl=0; pl<(2*NPLA); pl++){ |
316 |
// printf(" integrale: calorimetro %f \n",calorimetro[pl][1]); |
317 |
//calcolo intergale in unita di spessori di silicio |
318 |
Integrale += calorimetro[pl][1] * MIP;//piano di silicio |
319 |
// se non e'il 1o dopo l'Y (tutti i pari) c'e' il W |
320 |
if(pl%2!=0){ //equival W in Si |
321 |
Integrale+= 0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1]); |
322 |
}; |
323 |
}; |
324 |
//Integrale=24000;//Integrale*1000; |
325 |
Integrale *= 1000.; |
326 |
|
327 |
/*z ed energia con media troncata*/ |
328 |
// mediatroncata(); // out: 1)chi2, 2)z, 3)Etot, 4)Pskip |
329 |
|
330 |
/*z ed energia con loop*/ |
331 |
if ( debug ) printf(" call Zdaloop with integrale %f \n",Integrale); |
332 |
Zdaloop(); // out: 1)chi2, 2)z, 3)Etot, 4)Pskip |
333 |
|
334 |
|
335 |
if ( debug ) this->Print(); |
336 |
if ( debug ) printf(" fine evento \n"); |
337 |
// |
338 |
}; |
339 |
|
340 |
|
341 |
Float_t CaloBragg::Integral(){ |
342 |
Process(); |
343 |
|
344 |
Float_t dEpianiloop[44]; |
345 |
Int_t tz1=(Int_t)lpz; |
346 |
Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop |
347 |
|
348 |
|
349 |
Float_t integ = 0.; |
350 |
for(Int_t i=0;i<=estremi[1][0];i++){ |
351 |
// integ += dEplan[i]; |
352 |
//printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]); |
353 |
integ += dEpianiloop[i]; |
354 |
// printf(" step %i integ %f deplan %f \n",i,integ,dEpianiloop[i]); |
355 |
} |
356 |
return integ; |
357 |
} |
358 |
|
359 |
Float_t CaloBragg::LastIntegral(){ |
360 |
Process(); |
361 |
|
362 |
Float_t integ = 0.; |
363 |
for(Int_t i=0;i<=estremi[1][0];i++){ |
364 |
integ += dEplan[i]; |
365 |
//printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]); |
366 |
} |
367 |
return integ; |
368 |
} |
369 |
|
370 |
void CaloBragg::Draw(){ |
371 |
|
372 |
Process(); |
373 |
|
374 |
// Float_t dEpianimean[44]; |
375 |
Float_t dEpianiloop[44]; |
376 |
Float_t Depth[44]; |
377 |
// Int_t tz=(Int_t)qtz; |
378 |
Int_t tz1=(Int_t)lpz; |
379 |
// Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata |
380 |
Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop |
381 |
|
382 |
Float_t sp= spessore[0]*spessore[1]; |
383 |
for(Int_t i=0;i<44;i++)Depth[i]=i*sp; |
384 |
// |
385 |
gStyle->SetLabelSize(0.04); |
386 |
gStyle->SetNdivisions(510,"XY"); |
387 |
// |
388 |
TString hid = Form("cCaloBragg"); |
389 |
TCanvas *tc = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid)); |
390 |
if ( tc ){ |
391 |
// tc->Clear(); |
392 |
} else { |
393 |
tc = new TCanvas(hid,hid); |
394 |
// tc->Divide(1,2); |
395 |
}; |
396 |
// |
397 |
// TString thid = Form("hCaloBragg"); |
398 |
// TH2F *th = dynamic_cast<TH2F*>(gDirectory->FindObject(thid)); |
399 |
// if ( th ) th->Delete(); |
400 |
// th->Clear(); |
401 |
// th->Reset(); |
402 |
// } else { |
403 |
// th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.); |
404 |
// th->SetMarkerStyle(20); |
405 |
// }; |
406 |
// |
407 |
tc->cd(); |
408 |
TString thid2 = Form("hCaloBragg2"); |
409 |
TH2F *th2 = dynamic_cast<TH2F*>(gDirectory->FindObject(thid2)); |
410 |
if ( th2 ) th2->Delete(); |
411 |
th2 = new TH2F(thid2,thid2,300,-0.5,300.,1000,0.,150.); |
412 |
th2->SetMarkerStyle(20); |
413 |
th2->SetMarkerColor(kRed); |
414 |
// |
415 |
TString thid3 = Form("hCaloBragg3"); |
416 |
TH2F *th3 = dynamic_cast<TH2F*>(gDirectory->FindObject(thid3)); |
417 |
if ( th3 ) th3->Delete(); |
418 |
th3 = new TH2F(thid3,thid3,300,-0.5,300.,1000,0.,150.); |
419 |
th3->SetMarkerStyle(20); |
420 |
th3->SetMarkerColor(kBlue); |
421 |
|
422 |
|
423 |
// tc->cd(1); |
424 |
// |
425 |
// for(Int_t i=0;i<=estremi[1][0];i++)th->Fill(Depth[i],dEpianimean[i]); |
426 |
for(Int_t i=0;i<=estremi[1][0];i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP); |
427 |
// th->Draw(); |
428 |
th2->Draw("same"); |
429 |
|
430 |
// tc->cd(2); |
431 |
tc->cd(); |
432 |
// |
433 |
for(Int_t i=0;i<=estremi[1][0];i++)th3->Fill(Depth[i],dEpianiloop[i]); |
434 |
th3->Draw(); |
435 |
th2->Draw("same"); |
436 |
|
437 |
tc->Modified(); |
438 |
tc->Update(); |
439 |
|
440 |
// |
441 |
gStyle->SetLabelSize(0); |
442 |
gStyle->SetNdivisions(1,"XY"); |
443 |
// |
444 |
}; |
445 |
|
446 |
|
447 |
|
448 |
void CaloBragg::LoadParam(){ |
449 |
|
450 |
// |
451 |
elem[0] = 1.00794; //H 1 |
452 |
elem[1] = 4.002602; //He 2 |
453 |
elem[2] = 6.941; //Li 3 |
454 |
elem[3] = 9.012182;//Be 4 |
455 |
elem[4] = 10.811; //B 5 |
456 |
elem[5] = 12.0107; //C 6 |
457 |
elem[6] = 14.00674;//N 7 |
458 |
elem[7] = 15.9994; //O 8 |
459 |
elem[8] = 18.9984032; //F 9 |
460 |
elem[9] = 20.1797; //Ne 10 |
461 |
elem[10] = 22.98977;//Na 11 |
462 |
elem[11] = 24.3050; //Mg 12 |
463 |
elem[12] = 26.9815; //Al 13 |
464 |
elem[13] = 28.0855; //Si 14 |
465 |
elem[14] = 30.974; //P 15 |
466 |
elem[15] = 32.066; //S 16 |
467 |
elem[16] = 35.4527; //Cl 17 |
468 |
elem[17] = 39.948; //Ar 18 |
469 |
elem[18] = 39.0983; //K 19 |
470 |
elem[19] = 40.078; //Ca 20 |
471 |
elem[20] = 44.95591;//Sc 21 |
472 |
elem[21] = 47.867; //Ti 22 |
473 |
elem[22] = 50.9415; //V 23 |
474 |
elem[23] = 51.9961; //Cr 24 |
475 |
elem[24] = 54.938049;//Mn 25 |
476 |
elem[25] = 55.845; //Fe 26 |
477 |
elem[26] = 58.9332; //Co 27 |
478 |
elem[27] = 58.6934; //Ni 28 |
479 |
elem[28] = 63.546; //Cu 29 |
480 |
elem[29] = 65.39; //Zn 30 |
481 |
elem[30] = 69.723; //Ga 31 |
482 |
elem[31] = 72.61; //Ge 32 |
483 |
|
484 |
|
485 |
//parametri calorimetro |
486 |
NPLA = 22; |
487 |
NCHA = 96; |
488 |
nView = 2; |
489 |
|
490 |
AA = 0.96;//mm larghezza strip |
491 |
ADIST = 80.5;//mm distanza tra pad |
492 |
PIANO = 8.59;//mm distanza |
493 |
|
494 |
ySi = 0.38;//mm spessore silicio |
495 |
yW = 2.66;//mm spessore tungsteno |
496 |
rhoSi = 2.33;//g/cm3 densita' silicio |
497 |
rhoW = 19.3;//g/cm3 densita' tugsteno |
498 |
MIP = 0.106;//Mev g/cm2 energia al minimo nel silicio per 0.38 mm |
499 |
|
500 |
emin = 0.; |
501 |
|
502 |
//parametri bethe-bloch |
503 |
pigr = 3.1415; |
504 |
Na = 6.02e-23; |
505 |
ZA = 0.49; /*Z/A per Si*/ |
506 |
// ISi =182e-06; /*MeV*/ |
507 |
ISi = 171e-06; /*MeV*/ |
508 |
IW = 735e-06; /*MeV*/ |
509 |
// ISi =0.0001059994; /*GeV!!*/ no era giusto!! |
510 |
Me = 0.511; /* MeV*/ |
511 |
MassP = 931.27;/*MeV*/ |
512 |
r2 = 7.95e-26; /*ro*ro in cm */ |
513 |
|
514 |
}; |
515 |
|
516 |
|
517 |
|
518 |
// |
519 |
void CaloBragg::conversione(){ |
520 |
|
521 |
// calcolo spessore Si attraverato in funzione dell'inclinazione |
522 |
// e conversione dello spessore di W in Si e correzione del valore |
523 |
// della Mip pe lo spessore effettivo |
524 |
// |
525 |
// in : evento |
526 |
// |
527 |
// out: out[0] = gcm2Si = spessore silicio attraversato nel piano |
528 |
// out[1] = WinSi = spessore equivalente in Si del W attraversato |
529 |
// out[2] = Mip = fattore conversione energia riscalato allo spessore attrversatonel piano |
530 |
|
531 |
Float_t SiCross=0.; |
532 |
Float_t WCross = 0.; |
533 |
Float_t ytgx = 0; |
534 |
Float_t ytgy = 0; |
535 |
Float_t a = 0.; |
536 |
|
537 |
/*silicio*/ |
538 |
ytgx = ySi * L2->GetCaloLevel2()->tanx[0]; |
539 |
ytgy = ySi * L2->GetCaloLevel2()->tany[0]; |
540 |
|
541 |
//lunghezza effettiva di silicio attraversata (mm) |
542 |
SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy)); |
543 |
|
544 |
spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2 |
545 |
|
546 |
/*tungsteno*/ |
547 |
ytgx = yW * L2->GetCaloLevel2()->tanx[0]; |
548 |
ytgy = yW * L2->GetCaloLevel2()->tany[0]; |
549 |
|
550 |
//rapporto tra rilasci energetici nei due materiali |
551 |
WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad |
552 |
//gcm2W = WCross/10. * rhoW; |
553 |
|
554 |
// (g/cm2W)/(g/cm2Si) |
555 |
spessore[3] = (WCross/10.) * rhoW; |
556 |
a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664); //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664); |
557 |
spessore[1] = a; |
558 |
//riscala mip allo spessore attraversato |
559 |
spessore[2] = MIP*(SiCross/ySi); |
560 |
};//end conversione |
561 |
|
562 |
|
563 |
|
564 |
|
565 |
|
566 |
void CaloBragg::BetheBloch(Float_t *x, Float_t *z, Float_t *Mass, Float_t *gam, Float_t *Bet, Float_t *out, Float_t II){ |
567 |
|
568 |
//rilascio energetico con bethe bloch con correzioni |
569 |
//in: x: g/cm2 |
570 |
// z: carica |
571 |
// Mass: Massa uma |
572 |
// Ene: energia particella MeV//tolta |
573 |
// gam: (etot/massa) |
574 |
// Bet: rad((g2-1)/g2) |
575 |
// |
576 |
//out: energia rilasciata MeV |
577 |
|
578 |
|
579 |
Float_t eta =0.; |
580 |
Float_t Wmax =0.; |
581 |
Float_t lg =0.; |
582 |
Float_t Energia=0.; |
583 |
Float_t C=0.; |
584 |
Float_t INo = ISi; |
585 |
|
586 |
if ( usenewBB ) INo = II; |
587 |
|
588 |
eta = (*gam)*(*Bet); |
589 |
|
590 |
//Bet=3/gam; SQ(*gam) * SQ(*Bet) |
591 |
Wmax = 2.* Me * SQ(eta) / (1. + 2.*(*gam)*Me/(*Mass) + SQ(Me)/SQ(*Mass)); |
592 |
|
593 |
lg = 2.* Me * SQ(eta) * Wmax / SQ(INo); |
594 |
// Energia = x* 2 * pigr * Na * r2 * Me * rhoSi *ZA* SQ(z)/SQ(Bet) * lg; |
595 |
C=(0.42237*pow(eta,-2.) + 0.0304*pow(eta,-4.) - 0.00038*pow(eta,-6.))*pow(10.,-6.)* pow(INo,2.) + |
596 |
(3.858*pow(eta,-2.) - 0.1668*pow(eta,-4.) + 0.00158*pow(eta,-6.))*pow(10.,-9.)*pow(INo,3.); |
597 |
|
598 |
if(eta <= 0.13) C= C * log(eta/0.0653)/log(0.13/0.0653); |
599 |
|
600 |
Energia = (*x) * 0.307/28.09 * 14. *SQ(*z)/SQ(*Bet)*(0.5*log(lg) - SQ(*Bet) - C/14.); |
601 |
|
602 |
*out =Energia;//out |
603 |
|
604 |
};//end Bethebloch |
605 |
|
606 |
|
607 |
|
608 |
|
609 |
void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Float_t *Etot, Float_t *out, Float_t II){ |
610 |
|
611 |
/*perdita di energia per ioni pesanti (come da routine geant)*/ |
612 |
// in : dx => spessore g/cm2 |
613 |
// Z => carica |
614 |
// Etot => energia perticella |
615 |
// |
616 |
// out: energia persa |
617 |
|
618 |
|
619 |
Float_t Q=0.; |
620 |
Float_t v=0.; |
621 |
Float_t gam=0.; |
622 |
Float_t Bet=0.; |
623 |
Float_t dEP=0.; |
624 |
|
625 |
// gamma // Mass = A * MassP; /*in Mev/c2*/ |
626 |
gam = (*Etot)/(elem[*Z-1]*MassP); // E = gamma M c2 |
627 |
|
628 |
|
629 |
Bet = sqrt((SQ(gam) -1.)/SQ(gam)); |
630 |
|
631 |
// v= 121.4139*(Bet/pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet/pow((*Z),(2./3.)))); |
632 |
v= 121.4139*(Bet*pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet*pow((*Z),(2./3.)))); // EMI AAAAGGH!! |
633 |
|
634 |
//carica effettiva |
635 |
Q= (*Z)*(1- (1.034 - 0.1777*exp(-0.08114*(*Z)))*exp(-v)); |
636 |
|
637 |
//perdita energia per un protone |
638 |
Float_t protone =1.; |
639 |
// Float_t Mass=(elem[*Z-1]*MassP); //EMI |
640 |
// BetheBloch(dx, &protone, &Mass, &gam, &Bet, &dEP);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP); |
641 |
|
642 |
BetheBloch(dx, &protone, &MassP, &gam, &Bet, &dEP, II);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP); //EMI |
643 |
|
644 |
*out= (SQ(Q)*(dEP));//*dx; |
645 |
|
646 |
|
647 |
};//end ELOSS |
648 |
|
649 |
|
650 |
|
651 |
|
652 |
void CaloBragg::Enetrack(Int_t* Z, Float_t* E0, Float_t* primo,Float_t* ultimo, Float_t out[]){ |
653 |
|
654 |
//calcola energia rilasciata sulla traccia (usa ELOSS) |
655 |
// in : Z =>carica |
656 |
// E0 =>energia |
657 |
// spess2[3] => conversione spessore Si, Si in W, mip |
658 |
// primo => posizione primo piano attraversato |
659 |
// |
660 |
// out: array[44] =>rilasci energetici calcolati per ogni piano[44] dopo il primo(estremi[0][0]) |
661 |
|
662 |
|
663 |
|
664 |
Float_t dE=0.; //energia rilasciata |
665 |
Float_t Ezero= *E0;//energia iniziale |
666 |
|
667 |
//azzero energia rilasciata sui piani |
668 |
memset(out, 0, 2*NPLA*sizeof(Float_t)); |
669 |
|
670 |
Float_t Massa = (elem[(*Z)-1] * MassP); |
671 |
|
672 |
for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){ |
673 |
dE=0.; |
674 |
//spessore silicio corretto x inclinazione, z, energia, out:rilascio |
675 |
ELOSS(&spessore[0], Z, &Ezero, &dE, ISi);//spessore in g/cm2!! |
676 |
if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop |
677 |
out[ipla] = Ezero - Massa; //MeV |
678 |
return; |
679 |
|
680 |
}else{ |
681 |
out[ipla] = dE; //MeV |
682 |
Ezero = Ezero - dE;//energia residua |
683 |
if ( debug ) printf(" zompa %i out %f dE %f ezero %f \n",ipla,out[ipla],dE,Ezero); |
684 |
}; |
685 |
//se sono su un piano Y (tutti i pari) dopo c'e' il tungsteno |
686 |
if(ipla%2 == 0){ |
687 |
/*tungsteno*/ |
688 |
dE=0.; |
689 |
Float_t sp = 0.; |
690 |
Float_t II = ISi; |
691 |
if ( usenewBB ){ |
692 |
sp = spessore[3]; |
693 |
II = IW; |
694 |
} else { |
695 |
sp = spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato in g/cm2 |
696 |
} |
697 |
// printf(" sp %f II %f \n",sp,II); |
698 |
ELOSS(&sp, Z, &Ezero, &dE,II); |
699 |
if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop |
700 |
return; |
701 |
}else{ |
702 |
Ezero = Ezero -dE;//energia residua |
703 |
}; |
704 |
}; |
705 |
|
706 |
};//fine loop piani |
707 |
|
708 |
|
709 |
};//end Enetrack |
710 |
|
711 |
|
712 |
|
713 |
void CaloBragg::chiquadro(Float_t dE[], Float_t out[]){ |
714 |
|
715 |
// calcola chi2 tra energia calcolata e misurata |
716 |
// in : dE[44] =>energia calcolata |
717 |
// calo3[44][2]=> [0]strip attraversata [1]energia misurata per ogni piano |
718 |
// estr2 => array con primo[0][0] e ultimo[1][0] piano attraversati ed energie[][1] |
719 |
// |
720 |
// out: array[3]=> (chi2; piani scartati consecutivi(79= >3 quindi frammentato); piani scartati totale) |
721 |
|
722 |
|
723 |
Float_t sum = 0.; |
724 |
Float_t PianoPrecedente=0.; |
725 |
Float_t badplane=0.; |
726 |
Float_t badplanetot=0.; |
727 |
Float_t w,wi; |
728 |
// |
729 |
if ( newchi2 ){ |
730 |
ndf = 0; |
731 |
sum = 0.; |
732 |
for( Int_t ipla=((int)(estremi[0][0])); ipla<= ((int)(estremi[1][0])); ipla++){ |
733 |
sum += pow((dE[ipla] - (calorimetro[ipla][1] * spessore[2]))/(0.05*dE[ipla]),2.); |
734 |
// printf(" quiqui: dE %f calor %f spessore[2] %f \n",dE[ipla],spessore[2]*calorimetro[ipla][1],spessore[2]); |
735 |
ndf++; |
736 |
} |
737 |
ndf -= 2; |
738 |
if ( ndf > 0 ) sum /= (float)ndf; |
739 |
out[0] = sum; |
740 |
out[1] = 0.; |
741 |
out[2] = (int)(estremi[1][0])-ndf; |
742 |
// printf(" sum %f ndf %i \n ",sum,ndf); |
743 |
} else { |
744 |
for(Int_t ipla=0; ipla<2*NPLA; ipla++){ |
745 |
//tutti i piani attraversati dalla traiettoria |
746 |
if(calorimetro[ipla][0] != -1.){ // |
747 |
w=0.; //normalizzazione; |
748 |
wi=1.;//peso |
749 |
|
750 |
//tolgo piani attraversati dalla traccia ma precedenti il piano individuato come ingresso |
751 |
if (ipla<estremi[0][0]) wi=0.; |
752 |
|
753 |
//tolgo piani attraversati da traccia ma successivi all'ultimo se sono diversi da 0 |
754 |
//if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.; |
755 |
if((ipla>estremi[1][0])) wi=0.; |
756 |
|
757 |
//normalizzazione |
758 |
if (calorimetro[ipla][1] != 0.) w=1./(calorimetro[ipla][1]* MIP); // |
759 |
|
760 |
//tolgo piani con rilasci inferiori al 30% del precedente |
761 |
if(calorimetro[ipla][1] < (0.7*PianoPrecedente)){ // cosi' i piani senza rilascio non vengono considerati nel calcolo del chi2 |
762 |
wi=0.; |
763 |
//se sono piani intermedi (non si e' fermta) li considero non buoni |
764 |
if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=0.)){// |
765 |
badplane+=1.; |
766 |
badplanetot+=1.; |
767 |
}; |
768 |
}; |
769 |
|
770 |
//meno peso ai piani con rilasci maggiori di 1000 MIP |
771 |
// if(calorimetro[ipla][1] > 1000) wi=0.5; |
772 |
if(calorimetro[ipla][1] > 1200.) wi=0.5; |
773 |
|
774 |
Float_t arg = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP)); |
775 |
|
776 |
sum += SQ(arg); // w*wi*(dEpiani[p][v]-(eplane[p][v]*MIP))));//( dEpiani[p][v] - (eplane[p][v]*MIP)); |
777 |
if(debug){ |
778 |
printf("dedx calcolata %f e reale %f \n",dE[ipla],(calorimetro[ipla][1] * MIP)); |
779 |
} |
780 |
//se trovo piano non buono (tolto quindi wi=0) non modifico il piano precedente |
781 |
if(wi != 0.){// |
782 |
PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente |
783 |
badplane = 0.;//azzero contatore piani scartati consecutivi |
784 |
}; |
785 |
}; |
786 |
|
787 |
//da Emi |
788 |
if(badplane > 2){ |
789 |
// printf(" AAAAAAAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG\n"); |
790 |
out[1] =79.; |
791 |
break; |
792 |
}; |
793 |
|
794 |
};//fine loop piani |
795 |
//chi2,frammentato,pskip |
796 |
out[0]=sum; |
797 |
out[2]=badplanetot; |
798 |
} |
799 |
};//end chiquadro |
800 |
|
801 |
|
802 |
|
803 |
void CaloBragg::loopze( Float_t step, Float_t E0,Float_t Zstart, Float_t Zlimite, Int_t nostep = 1000){ |
804 |
// |
805 |
//loop su z ed energie per trovare miglior z (ed energia) |
806 |
//in: nloop => energia massima da provare (nloop x E0) |
807 |
// E0 => energia iniziale (intergale) |
808 |
// Zstart => minimo z da cui patire |
809 |
// Zlimite => z a cui fermarsi (z al minimo di ionizz sul 1o piano) |
810 |
// |
811 |
//out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2 |
812 |
// |
813 |
|
814 |
|
815 |
memset(dEplan,0,2*NPLA*sizeof(Float_t)); |
816 |
|
817 |
Int_t Z = 0;// z iniziale |
818 |
|
819 |
Float_t Massa = 0.; |
820 |
|
821 |
Float_t Stepint =(step)/(Float_t)nostep;//passo per il calcolo di energia |
822 |
|
823 |
Float_t energia =0.;//energia del loop |
824 |
|
825 |
Float_t chi2[3] = {0,0,0};//out dal calcolo chi2: chi2, piani consecutivi saltati, piani totali saltati |
826 |
|
827 |
Int_t zmin = (int)Zstart; |
828 |
Int_t max=32;//max z di cui so la massa :P |
829 |
if((Zlimite)<=31) max=(int)(Zlimite) + 1; |
830 |
|
831 |
if ( fzeta > 0. ){ |
832 |
zmin = fzeta; |
833 |
max = fzeta+1; |
834 |
} |
835 |
|
836 |
Int_t colmax=32; |
837 |
Int_t rowmax=3000; |
838 |
|
839 |
Float_t matrixchi2[colmax][rowmax][3]; |
840 |
memset(matrixchi2, 0, colmax*rowmax*3*sizeof(Float_t)); |
841 |
|
842 |
Int_t imin = 1-nostep/2; |
843 |
Int_t imax = nostep/2; |
844 |
|
845 |
//loop elementi |
846 |
for(Int_t inucl=zmin; inucl<max; inucl++){ |
847 |
|
848 |
Z= inucl; |
849 |
|
850 |
Massa = elem[inucl-1]*MassP; |
851 |
|
852 |
//loop energia |
853 |
Int_t iene2 = 0; |
854 |
// for(Int_t iene= 0; iene<1000; iene++){// da non cambiare in base a Stepint altrimenti cambia la matrice bestchi2!!!cosi' non raggiungo mai integrale!!!!! mettere <=?? |
855 |
for(Int_t iene= imin; iene<imax; iene++){// da non cambiare in base a Stepint altrimenti cambia la matrice bestchi2!!!cosi' non raggiungo mai integrale!!!!! mettere <=?? |
856 |
|
857 |
iene2++; |
858 |
energia= Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica.. |
859 |
|
860 |
Enetrack(&Z, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani |
861 |
//calcolo chi2 |
862 |
chiquadro(dEplan,chi2); |
863 |
|
864 |
// printf(" last deplan from: Z = %i iene %i energia %f chi2 %f \n",inucl,iene,energia,chi2[0]); |
865 |
|
866 |
if( (chi2[1] != 79.) ){//salto quelli che frammentano |
867 |
matrixchi2[inucl][iene2][0]=chi2[0];//valore chi2 per questo z a questa energia |
868 |
matrixchi2[inucl][iene2][1]=energia;//energia per questo chi2 |
869 |
matrixchi2[inucl][iene2][2]=chi2[2];//piani saltati nel chi2 |
870 |
} else { |
871 |
matrixchi2[inucl][iene2][0]=1000.;//valore chi2 per questo z a questa energia |
872 |
matrixchi2[inucl][iene2][1]=1000.;//energia per questo chi2 |
873 |
matrixchi2[inucl][iene2][2]=1000.;//piani saltati nel chi2 |
874 |
break; |
875 |
} |
876 |
}//fine loop energia |
877 |
|
878 |
|
879 |
};//fine loop z |
880 |
|
881 |
|
882 |
//Emi |
883 |
for (Int_t nu=zmin; nu<max; nu++){ |
884 |
for (Int_t en=0; en<nostep; en++){ |
885 |
if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] >0.)){ |
886 |
bestchi2[0]= matrixchi2[nu][en][0];// chi2 |
887 |
bestchi2[1]= (Float_t)nu; // z |
888 |
bestchi2[2]= matrixchi2[nu][en][1];//energia; |
889 |
bestchi2[3]= matrixchi2[nu][en][2];// totale piani saltati |
890 |
} |
891 |
} |
892 |
} |
893 |
|
894 |
|
895 |
};//endloopze |
896 |
|
897 |
|
898 |
|
899 |
|
900 |
|
901 |
// void CaloBragg::mediatroncata(){ |
902 |
// //calcolo Z con media troncata e utilizzo questo Z per trovare l'energia migliore |
903 |
// //in: ordplane[44] => array con energia dei piani |
904 |
// // spess[3] => conversioni spessore di silicio, w, mip |
905 |
// // estr[2][2] => primo[0][0] e ultimo[1][0] piano attraversati ed energie[][1] |
906 |
// // calo[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano |
907 |
// // integrale => energia totale nel calorimetro considerando il W |
908 |
// // |
909 |
// // out[4] chi2,z,Etot,Pskip |
910 |
|
911 |
// Float_t ordplane[44];//mi serve per la media troncata |
912 |
// memset(ordplane,0,44*sizeof(Float_t)); |
913 |
|
914 |
// for(Int_t ipla=0; ipla< 2*NPLA; ipla++) ordplane[ipla]=calorimetro[ipla][1]; //energia del piano |
915 |
|
916 |
|
917 |
// //ordino tutte le energie dei piani in ordine crescente |
918 |
|
919 |
// Long64_t work[200]; |
920 |
// Int_t ind = 0; |
921 |
// //Int_t l = 0; |
922 |
// Int_t RN = 0; |
923 |
// Float_t sum4 = 0.; |
924 |
// Float_t qm = 0.; |
925 |
// // |
926 |
// //Float_t qmt = ethr*0.8; // *0.9 |
927 |
// // |
928 |
// //Int_t uplim = TMath::Max(3,N); |
929 |
// // |
930 |
// while ( RN < 4 && ind < 44 ){ |
931 |
// qm = TMath::KOrdStat(44,ordplane,ind,work); |
932 |
// if (qm >= 0.7 ){ |
933 |
// if ( RN < 4 ){ |
934 |
// sum4 += qm; |
935 |
// RN++; |
936 |
// }; |
937 |
// // l++; |
938 |
// // if ( debug ) printf(" value no %i qm %f sum4 %f \n",l,qm,sum4); |
939 |
// }; |
940 |
// ind++; |
941 |
// }; |
942 |
// // |
943 |
// sum4 /= (Float_t)RN; |
944 |
// Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));//ma non e'/1?? |
945 |
// if(Zmean ==0.) Zmean=1.; |
946 |
// if ( Zmean < 1. ) Zmean = 1.; |
947 |
|
948 |
|
949 |
// // Zmean =round(Zmean); |
950 |
// // if(Zmean <1.) Zmean=1.; |
951 |
|
952 |
// // if(Zmean >0.)Zmean =round(Zmean); |
953 |
|
954 |
// //======== per i nuclei======= |
955 |
// if (Zmean >=2.){ |
956 |
// ind = 0; |
957 |
// RN = 0; |
958 |
// sum4 = 0.; |
959 |
// qm = 0.; |
960 |
// while ( RN < 4 && ind < 44 ){ |
961 |
// qm = TMath::KOrdStat(44,ordplane,ind,work); |
962 |
// if (qm >= (Zmean*Zmean)-Zmean*Zmean*0.2 ){ |
963 |
// if ( RN < 4 ){ |
964 |
// sum4 += qm; |
965 |
// RN++; |
966 |
// }; |
967 |
// }; |
968 |
// ind++; |
969 |
// }; |
970 |
// // |
971 |
// sum4 /= (Float_t)RN; |
972 |
// Zmean = (sqrt((sum4*MIP)/(4.*spessore[2])));//ma non e' /1?? |
973 |
// } |
974 |
|
975 |
|
976 |
// //calcolo energia migliore per Z trovato con media troncata |
977 |
// // Float_t zmin=Zmean; |
978 |
// Float_t zmin=round(Zmean); |
979 |
|
980 |
// bestchi2[0]=10000.; |
981 |
// bestchi2[1]=0.; |
982 |
// bestchi2[2]=0.; |
983 |
// bestchi2[3]=0.; |
984 |
// Float_t zero=0.; |
985 |
|
986 |
// // step energia zstart zstop |
987 |
// loopze(Integrale,zero,zmin,zmin); |
988 |
|
989 |
|
990 |
// qtchi2=bestchi2[0]; |
991 |
// qtz=bestchi2[1]; |
992 |
// qtetot=bestchi2[2]; |
993 |
// qtpskip=bestchi2[3]; |
994 |
// };//end mediatroncata |
995 |
|
996 |
|
997 |
|
998 |
void CaloBragg::Zdaloop(){ |
999 |
//calcolo Z con un loop su tutti i possibli Z ed energie |
1000 |
//in: ordplane[44]=> array con energia dei piani |
1001 |
// spess1[3]=> conversioni spessore di silicio, w e mip |
1002 |
// estr3[2][2]=> primo[0][0] e ultimo[1][0] piano ed energie |
1003 |
// calo1[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano |
1004 |
// integrale=> energia totale nel calorimetro considerando il W |
1005 |
// |
1006 |
// out[4] chi2,z,Etot,Pskip |
1007 |
|
1008 |
|
1009 |
/*z se particella fosse al minimo*/ //energia1piano/mip corretta |
1010 |
// Float_t zmax = round(sqrt(estremi[0][1]/spessore[2])); |
1011 |
// if(zmax<31)zmax=zmax+1; |
1012 |
|
1013 |
/*calcolo Z ed E con loop sui vari elementi ed energie*/ |
1014 |
|
1015 |
Float_t zmin=1.; |
1016 |
Float_t zmax=32.; |
1017 |
Float_t bestchitemp[4] = {0,0,0,0}; |
1018 |
|
1019 |
bestchi2[0]=10000.; |
1020 |
bestchi2[1]=0.; |
1021 |
bestchi2[2]=0.; |
1022 |
bestchi2[3]=0.; |
1023 |
Float_t zero=0.; |
1024 |
//------------primo loop ---------------------- |
1025 |
// energia ezero, zstart zstop |
1026 |
// loopze(Integrale,zero,zmin,zmax); |
1027 |
|
1028 |
//-> loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,50); |
1029 |
loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,200); |
1030 |
|
1031 |
// loopze(Integrale*2.,Integrale/100.,zmin,zmax); |
1032 |
if ( debug ) printf(" Integrale %f , outene %f \n",Integrale,bestchi2[2]); |
1033 |
|
1034 |
//------------secondo loop ---------------------- |
1035 |
for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i]; |
1036 |
bestchi2[0] = 10000.; |
1037 |
bestchi2[1] = 0.; |
1038 |
bestchi2[2] = 0.; |
1039 |
bestchi2[3] = 0.;//riazzero |
1040 |
|
1041 |
Float_t step = bestchitemp[2];// |
1042 |
zero=0.; // qualsiasi altro valore peggiora le cose |
1043 |
// zmin=zmax=bestchitemp[1]; |
1044 |
zmin=bestchitemp[1]-1; |
1045 |
zmax=bestchitemp[1]+1; |
1046 |
// loopze(step,zero,zmin,zmax); // |
1047 |
|
1048 |
//-> loopze(step,step/2.,zmin,zmax,200); // |
1049 |
loopze(step,step/2.,zmin,zmax,500); // |
1050 |
|
1051 |
if ( debug ) printf(" Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step); |
1052 |
|
1053 |
|
1054 |
//chi2,z,Etot,Pskip |
1055 |
lpchi2=bestchi2[0]; |
1056 |
lpz=bestchi2[1]; |
1057 |
lpetot=bestchi2[2]; |
1058 |
lppskip=bestchi2[3]; |
1059 |
|
1060 |
};//endZdaloop |
1061 |
|
1062 |
|
1063 |
|
1064 |
|
1065 |
|
1066 |
|
1067 |
|
1068 |
|
1069 |
|
1070 |
|
1071 |
|
1072 |
|