| 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 |
lpisotope= 0.; |
| 47 |
lpetot = 0.; |
| 48 |
lppskip = 0.; |
| 49 |
|
| 50 |
memset(calorimetro,0,44*2*sizeof(Float_t)); |
| 51 |
memset(spessore,0,4*sizeof(Float_t)); |
| 52 |
memset(estremi,0,2*2*sizeof(Float_t)); |
| 53 |
Integrale=0.; |
| 54 |
|
| 55 |
for(Int_t l=0;l<44;l++){ |
| 56 |
calorimetro[l][0]=-1.; |
| 57 |
} |
| 58 |
|
| 59 |
}; |
| 60 |
|
| 61 |
void CaloBragg::Print(){ |
| 62 |
// |
| 63 |
|
| 64 |
if(!debug) Process(); |
| 65 |
// |
| 66 |
printf("========================================================================\n"); |
| 67 |
printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack); |
| 68 |
printf(" first plane: %f \n", estremi[0][0]); |
| 69 |
printf(" last plane: %f \n", estremi[1][0]); |
| 70 |
// printf(" chi 2 from truncated mean: %f \n", qtchi2); |
| 71 |
// printf(" Z from truncated mean %f: \n", qtz); |
| 72 |
// printf(" energy from truncated mean %f: \n", qtetot); |
| 73 |
// printf(" plane not used for truncated mean %f: \n", qtpskip); |
| 74 |
printf(" chi 2 from loop %f: \n", lpchi2); |
| 75 |
printf(" Z from loop %f: \n", lpz); |
| 76 |
printf(" isotope from loop %f: \n", lpisotope); |
| 77 |
printf(" energy from loop %f: \n", lpetot); |
| 78 |
printf(" plane not used for loop %f: \n", lppskip); |
| 79 |
printf(" ndf: %i \n",ndf); |
| 80 |
printf("========================================================================\n"); |
| 81 |
// |
| 82 |
}; |
| 83 |
|
| 84 |
void CaloBragg::Delete(){ |
| 85 |
Clear(); |
| 86 |
//delete this; |
| 87 |
}; |
| 88 |
|
| 89 |
|
| 90 |
void CaloBragg::Process(){ |
| 91 |
Process(-1); |
| 92 |
}; |
| 93 |
|
| 94 |
|
| 95 |
void CaloBragg::CleanPlanes(Float_t epiano[22][2], Bool_t zpiano[22][2]){ |
| 96 |
// return; |
| 97 |
Int_t hitplanes = 0; |
| 98 |
Float_t f5 = 0.; |
| 99 |
for (Int_t i = 0; i<22; i++){ |
| 100 |
for (Int_t j = 1; j>=0; j--){ |
| 101 |
zpiano[i][j] = false; |
| 102 |
if ( epiano[i][j] > 0.7 ){ |
| 103 |
if ( hitplanes < 100 ) f5 += epiano[i][j]; |
| 104 |
hitplanes++; |
| 105 |
}; |
| 106 |
}; |
| 107 |
}; |
| 108 |
Int_t atl5 = TMath::Min(hitplanes,100); |
| 109 |
atl5 = TMath::Max(atl5,1); |
| 110 |
Float_t lowlim = 0.85; |
| 111 |
//Float_t lowlim = 1.; |
| 112 |
Float_t dedxone = 0.; |
| 113 |
// Float_t step1 = 0.8*L2->GetCaloLevel2()->qtot/(Float_t)hitplanes; |
| 114 |
Float_t step1 = 0.8*f5/atl5; |
| 115 |
// while ( dedxone < step1 ){ |
| 116 |
for (Int_t i = 0; i<22; i++){ |
| 117 |
for (Int_t j = 1; j>=0; j--){ |
| 118 |
if (debug) printf("Acleanplanes: i %i j %i step1 %f dedxone %f epiano[i][j] %f \n",i,j,step1,dedxone,epiano[i][j]); |
| 119 |
if ( epiano[i][j] >= step1 && dedxone < 0.7 ) dedxone = epiano[i][j]; |
| 120 |
if ( dedxone >= step1 ) break; // new |
| 121 |
}; |
| 122 |
if ( dedxone >= step1 ) break; // new |
| 123 |
}; |
| 124 |
// } |
| 125 |
if ( dedxone < 0.7 ){ // here we could have instead while dedxone == 0. ... perhaps better... |
| 126 |
for (Int_t i = 0; i<22; i++){ |
| 127 |
for (Int_t j = 1; j>=0; j--){ |
| 128 |
if (debug) printf("Bcleanplanes dedxone < 0.7: i %i j %i step1 %f dedxone %f epiano[i][j] %f \n",i,j,step1,dedxone,epiano[i][j]); |
| 129 |
if ( epiano[i][j] > 0. && dedxone < 0.7 ) dedxone = epiano[i][j]; |
| 130 |
if ( dedxone >= 0.7 ) break; // new |
| 131 |
}; |
| 132 |
if ( dedxone >= 0.7 ) break; // new |
| 133 |
}; |
| 134 |
} |
| 135 |
// |
| 136 |
// printf(" dedxone = %f step1 %f \n",dedxone,step1); |
| 137 |
Bool_t revulsera = false; |
| 138 |
Bool_t nullius = false; |
| 139 |
Int_t nulliferus = 0; |
| 140 |
for (Int_t i = 0; i<22; i++){ |
| 141 |
for (Int_t j = 1; j>=0; j--){ |
| 142 |
if ( epiano[i][j] < dedxone*lowlim ){ |
| 143 |
if ( debug ) printf("Ccleanplanes: %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius); |
| 144 |
// epiano[i][j] = 0.; |
| 145 |
zpiano[i][j] = true; |
| 146 |
if ( epiano[i][j] < dedxone*0.05 ) epiano[i][j] = 0.; |
| 147 |
} else { |
| 148 |
if ( debug ) printf("Dcleanplanes else: %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius); |
| 149 |
nulliferus = 0; |
| 150 |
revulsera = true; |
| 151 |
}; |
| 152 |
// if ( epiano[i][j] < 0.7 && revulsera ) nulliferus++; |
| 153 |
if ( (zpiano[i][j] || epiano[i][j] < 0.7 ) && revulsera ) nulliferus++; |
| 154 |
if ( nulliferus > 10 ) nullius = true; |
| 155 |
// if ( nullius ) epiano[i][j] = 0.; |
| 156 |
if ( nullius ) zpiano[i][j] = true; |
| 157 |
}; |
| 158 |
}; |
| 159 |
|
| 160 |
} |
| 161 |
|
| 162 |
|
| 163 |
void CaloBragg::Process(Int_t ntr){ |
| 164 |
// |
| 165 |
if ( !L2 ){ |
| 166 |
printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n"); |
| 167 |
printf(" ERROR: CaloBragg variables not filled \n"); |
| 168 |
return; |
| 169 |
}; |
| 170 |
// |
| 171 |
Bool_t newentry = false; |
| 172 |
// |
| 173 |
if ( L2->IsORB() ){ |
| 174 |
if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || ntr != sntr ){ |
| 175 |
newentry = true; |
| 176 |
OBT = L2->GetOrbitalInfo()->OBT; |
| 177 |
PKT = L2->GetOrbitalInfo()->pkt_num; |
| 178 |
atime = L2->GetOrbitalInfo()->absTime; |
| 179 |
sntr = ntr; |
| 180 |
}; |
| 181 |
} else { |
| 182 |
newentry = true; |
| 183 |
}; |
| 184 |
// |
| 185 |
if ( !newentry ) return; |
| 186 |
// |
| 187 |
tr = ntr; |
| 188 |
// |
| 189 |
if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime); |
| 190 |
// |
| 191 |
Clear(); |
| 192 |
|
| 193 |
// |
| 194 |
// |
| 195 |
// |
| 196 |
Int_t view = 0; |
| 197 |
Int_t plane = 0; |
| 198 |
Int_t strip = 0; |
| 199 |
Float_t mip = 0.; |
| 200 |
Float_t epiano[22][2]; |
| 201 |
memset(epiano,0,22*2*sizeof(Float_t)); |
| 202 |
for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){ |
| 203 |
// |
| 204 |
mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); |
| 205 |
// |
| 206 |
if ( !usepl18x && view==0 && plane==18 ) mip = 0.; |
| 207 |
// |
| 208 |
epiano[plane][view]+=mip; |
| 209 |
// |
| 210 |
// |
| 211 |
}; |
| 212 |
// |
| 213 |
Bool_t zpiano[22][2]; |
| 214 |
this->CleanPlanes(*&epiano, *&zpiano); |
| 215 |
// |
| 216 |
PamTrack *ptrack = 0; |
| 217 |
CaloTrkVar *track = 0; |
| 218 |
// |
| 219 |
if ( usetrack ){ |
| 220 |
if ( ntr >= 0 ){ |
| 221 |
ptrack = L2->GetTrack(ntr); |
| 222 |
if ( ptrack ) track = ptrack->GetCaloTrack(); |
| 223 |
} else { |
| 224 |
track = L2->GetCaloStoredTrack(ntr); //al momento e' vera solo questa riga |
| 225 |
}; |
| 226 |
// |
| 227 |
if ( !track && ntr >= 0 ){ |
| 228 |
printf(" ERROR: cannot find any track!\n"); |
| 229 |
printf(" ERROR: CaloBragg variables not completely filled \n"); |
| 230 |
return; |
| 231 |
}; |
| 232 |
} else { |
| 233 |
if ( ntr >= 0 ){ |
| 234 |
if ( debug ) printf(" ERROR: you asked not to use a track but you are looking for track number %i !\n",ntr); |
| 235 |
if ( debug ) printf(" ERROR: CaloBragg variables not completely filled \n"); |
| 236 |
return; |
| 237 |
}; |
| 238 |
}; |
| 239 |
// |
| 240 |
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 |
| 241 |
|
| 242 |
// |
| 243 |
Bool_t zcalo[44]; |
| 244 |
for(Int_t p=0; p<22; p++){ |
| 245 |
for(Int_t v=0; v<2; v++){ |
| 246 |
/*per usare traccia non del calo camboare cibar*/ |
| 247 |
calorimetro[(2*p)+1-v][0] = L2->GetCaloLevel2()->cibar[p][v];//strip attraversata |
| 248 |
calorimetro[(2*p)+1-v][1] = epiano[p][v]; //energia del piano //(epiano[p][v])/0.89 |
| 249 |
zcalo[(2*p)+1-v] = zpiano[p][v]; |
| 250 |
if ( debug ) printf(" idx %i %f %i \n",(2*p)+1-v,epiano[p][v], zpiano[p][v]); |
| 251 |
}; |
| 252 |
}; |
| 253 |
|
| 254 |
/*per ogni evento calcolo la conversione mip e w attraversato in equivalente Si*/ |
| 255 |
conversione(); // out: 1) g/cm2 Si , 2) spessoreW equivalente in Si, 3)Mip corretta per inclinazione |
| 256 |
|
| 257 |
/*settaggio della soglia per il loop sulla determinazione del piano di partenza */ |
| 258 |
Float_t ordplane[44];//mi serve per la media troncata |
| 259 |
memset(ordplane,0,44*sizeof(Float_t)); |
| 260 |
|
| 261 |
for(Int_t ipla=0; ipla< 2*NPLA; ipla++) ordplane[ipla]=calorimetro[ipla][1]; //energia del piano |
| 262 |
|
| 263 |
|
| 264 |
//ordino tutte le energie dei piani in ordine crescente |
| 265 |
|
| 266 |
Long64_t work[200]; |
| 267 |
Int_t ind = 0; |
| 268 |
//Int_t l = 0; |
| 269 |
Int_t RN = 0; |
| 270 |
Float_t sum4 = 0.; |
| 271 |
Float_t qm = 0.; |
| 272 |
while ( RN < 4 && ind < 44 ){ |
| 273 |
qm = TMath::KOrdStat((Long64_t)44,ordplane,(Long64_t)ind,work); |
| 274 |
if (qm >= 0.7 ){ |
| 275 |
if ( RN < 4 ){ |
| 276 |
sum4 += qm; |
| 277 |
RN++; |
| 278 |
}; |
| 279 |
}; |
| 280 |
ind++; |
| 281 |
}; |
| 282 |
// |
| 283 |
//sum4 /= (Float_t)RN; |
| 284 |
Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2]))); |
| 285 |
if(Zmean ==0.) Zmean=1.; |
| 286 |
if ( Zmean < 1. ) Zmean = 1.; |
| 287 |
|
| 288 |
|
| 289 |
/*trova primo e ultimo piano attraversati*/ |
| 290 |
Int_t p = 0;//contatore piani |
| 291 |
//per il primo parte da 0 e va in giu' |
| 292 |
while( estremi[0][1] <= 0. && p<(2*NPLA) ){ // era ==0 ma ricorda i problemi con Float == !!!!! |
| 293 |
// if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){ |
| 294 |
// if( (calorimetro[p][0] >0) && (calorimetro[p][1]*MIP >0.3)){ //0.7 mip = 70MeV soglia minima |
| 295 |
if( (calorimetro[p][0] >0) && (calorimetro[p][1]*MIP >Zmean*0.7)){ // 70% della MIP |
| 296 |
estremi[0][0]=p; |
| 297 |
estremi[0][1]=calorimetro[p][1] *MIP; //energia in MeV |
| 298 |
}; |
| 299 |
p++; |
| 300 |
}; |
| 301 |
|
| 302 |
//ultimo parte da 44 e sale |
| 303 |
p=43; |
| 304 |
while( (estremi[1][1] <= 0.) && (p>(int)estremi[0][0]) ){ |
| 305 |
if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){ |
| 306 |
estremi[1][0]=p;// |
| 307 |
estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV |
| 308 |
}; |
| 309 |
p = p-1; |
| 310 |
}; |
| 311 |
// |
| 312 |
|
| 313 |
Float_t lastok = 0.; |
| 314 |
// if ( false ){ |
| 315 |
// Bool_t goback = false; |
| 316 |
for ( int o = 0; o < estremi[1][0]; o++ ){ |
| 317 |
// |
| 318 |
if (debug) printf(" goforth1: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok); |
| 319 |
if ( calorimetro[o][1] > 0.7 && !zcalo[o] ) lastok = calorimetro[o][1]; |
| 320 |
if ( (zcalo[o] || calorimetro[o][1] < 0.7) && lastok > 0. ){ |
| 321 |
if ( fabs(calorimetro[o][1]-lastok)/calorimetro[o][1] > 0.5 ) { |
| 322 |
if (debug) printf(" goforthchange %f %f \n",calorimetro[o][1],lastok); |
| 323 |
calorimetro[o][1] = lastok; |
| 324 |
if (debug) printf(" goforthchang+ %f %f \n",calorimetro[o][1],lastok); |
| 325 |
} |
| 326 |
} |
| 327 |
if (debug) printf(" goforth2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok); |
| 328 |
// if ( calorimetro[o][1] < 0.7 ) goback = true; |
| 329 |
// |
| 330 |
}; |
| 331 |
lastok = 0.; |
| 332 |
// if ( goback ){ |
| 333 |
for ( int o = estremi[1][0]; o >= 0; o-- ){ |
| 334 |
// |
| 335 |
if (debug) printf(" goback1: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok); |
| 336 |
if ( o < estremi[1][0] && calorimetro[o][1] > calorimetro[o+1][1]*1.2 && lastok > 0. ) calorimetro[o][1] = lastok; |
| 337 |
if ( calorimetro[o][1] > 0.7 && !zcalo[o] ) lastok = calorimetro[o][1]; |
| 338 |
if ( (zcalo[o] || calorimetro[o][1] < 0.7) && lastok > 0. ){ |
| 339 |
if ( fabs(calorimetro[o][1]-lastok)/calorimetro[o][1] > 0.5 ) { |
| 340 |
if ( debug ) printf(" gobackchange \n"); |
| 341 |
calorimetro[o][1] = lastok; |
| 342 |
} |
| 343 |
} |
| 344 |
if (debug) printf(" goback2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok); |
| 345 |
// |
| 346 |
}; |
| 347 |
// }; |
| 348 |
//} |
| 349 |
|
| 350 |
if ( startZero ) { |
| 351 |
estremi[0][0] = 0.; |
| 352 |
// estremi[0][1] = 0.; |
| 353 |
} |
| 354 |
|
| 355 |
/*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/ |
| 356 |
for(Int_t pl=0; pl<(2*NPLA); pl++){ |
| 357 |
// printf(" integrale: calorimetro %f \n",calorimetro[pl][1]); |
| 358 |
//calcolo intergale in unita di spessori di silicio |
| 359 |
Integrale += calorimetro[pl][1] * MIP;//piano di silicio |
| 360 |
// se non e'il 1o dopo l'Y (tutti i pari) c'e' il W |
| 361 |
if(pl%2!=0){ //equival W in Si |
| 362 |
Integrale+= 0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1]); |
| 363 |
}; |
| 364 |
}; |
| 365 |
//Integrale=24000;//Integrale*1000; |
| 366 |
Integrale *= 1000.; |
| 367 |
|
| 368 |
/*z ed energia con media troncata*/ |
| 369 |
// mediatroncata(); // out: 1)chi2, 2)z, 3)Etot, 4)Pskip |
| 370 |
|
| 371 |
/*z ed energia con loop*/ |
| 372 |
if ( debug ) printf(" call Zdaloop with integrale %f \n",Integrale); |
| 373 |
Zdaloop(); // out: 1)chi2, 2)z, 3)Etot, 4)Pskip |
| 374 |
|
| 375 |
|
| 376 |
if ( debug ) this->Print(); |
| 377 |
if ( debug ) printf(" fine evento \n"); |
| 378 |
// |
| 379 |
}; |
| 380 |
|
| 381 |
|
| 382 |
Float_t CaloBragg::Integral(){ |
| 383 |
Process(); |
| 384 |
|
| 385 |
Float_t dEpianiloop[44]; |
| 386 |
Int_t tz1=(Int_t)lpz; |
| 387 |
Int_t ti1=(Int_t)lpisotope; |
| 388 |
|
| 389 |
Enetrack(&tz1, &ti1 , &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop |
| 390 |
|
| 391 |
|
| 392 |
Float_t integ = 0.; |
| 393 |
for(Int_t i=0;i<=estremi[1][0];i++){ |
| 394 |
// integ += dEplan[i]; |
| 395 |
//printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]); |
| 396 |
integ += dEpianiloop[i]; |
| 397 |
// printf(" step %i integ %f deplan %f \n",i,integ,dEpianiloop[i]); |
| 398 |
} |
| 399 |
return integ; |
| 400 |
} |
| 401 |
|
| 402 |
Float_t CaloBragg::LastIntegral(){ |
| 403 |
Process(); |
| 404 |
|
| 405 |
Float_t integ = 0.; |
| 406 |
for(Int_t i=0;i<=estremi[1][0];i++){ |
| 407 |
integ += dEplan[i]; |
| 408 |
//printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]); |
| 409 |
} |
| 410 |
return integ; |
| 411 |
} |
| 412 |
|
| 413 |
|
| 414 |
void CaloBragg::Draw(){ |
| 415 |
|
| 416 |
Process(); |
| 417 |
|
| 418 |
this->Draw(0.,0.,0.); |
| 419 |
|
| 420 |
} |
| 421 |
|
| 422 |
void CaloBragg::Draw(Int_t Z, Int_t isotope, Float_t enetot){ |
| 423 |
|
| 424 |
// Float_t dEpianimean[44]; |
| 425 |
Float_t dEpianiloop[44]; |
| 426 |
Float_t Depth[44]; |
| 427 |
// Int_t tz=(Int_t)qtz; |
| 428 |
Int_t tz1= Z; |
| 429 |
Int_t ti1= isotope; |
| 430 |
Float_t enet = enetot; |
| 431 |
// Float_t enet = lpetot; |
| 432 |
|
| 433 |
if ( Z > 0. && enetot > 0. ){ |
| 434 |
estremi[0][0] = 0; |
| 435 |
estremi[1][0] = 43; |
| 436 |
|
| 437 |
|
| 438 |
Float_t ytgx = 0.; |
| 439 |
Float_t ytgy = 0.; |
| 440 |
|
| 441 |
//lunghezza effettiva di silicio attraversata (mm) |
| 442 |
Float_t SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy)); |
| 443 |
|
| 444 |
spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2 |
| 445 |
|
| 446 |
/*tungsteno*/ |
| 447 |
|
| 448 |
//rapporto tra rilasci energetici nei due materiali |
| 449 |
Float_t WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad |
| 450 |
//gcm2W = WCross/10. * rhoW; |
| 451 |
|
| 452 |
// (g/cm2W)/(g/cm2Si) |
| 453 |
spessore[3] = (WCross/10.) * rhoW; |
| 454 |
Float_t a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664); //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664); |
| 455 |
spessore[1] = a; |
| 456 |
//riscala mip allo spessore attraversato |
| 457 |
spessore[2] = MIP*(SiCross/ySi); |
| 458 |
|
| 459 |
} else { |
| 460 |
tz1=(Int_t)lpz; |
| 461 |
ti1=(Int_t)lpisotope; |
| 462 |
enet = lpetot; |
| 463 |
// Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata |
| 464 |
|
| 465 |
} |
| 466 |
Enetrack(&tz1, &ti1, &enet, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop |
| 467 |
|
| 468 |
Float_t sp= spessore[0]*spessore[1]; |
| 469 |
for(Int_t i=0;i<44;i++)Depth[i]=i*sp; |
| 470 |
// |
| 471 |
gStyle->SetLabelSize(0.04); |
| 472 |
gStyle->SetNdivisions(510,"XY"); |
| 473 |
// |
| 474 |
TString hid = Form("cCaloBragg"); |
| 475 |
TCanvas *tc = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid)); |
| 476 |
if ( tc ){ |
| 477 |
// tc->Clear(); |
| 478 |
} else { |
| 479 |
tc = new TCanvas(hid,hid); |
| 480 |
// tc->Divide(1,2); |
| 481 |
}; |
| 482 |
// |
| 483 |
// TString thid = Form("hCaloBragg"); |
| 484 |
// TH2F *th = dynamic_cast<TH2F*>(gDirectory->FindObject(thid)); |
| 485 |
// if ( th ) th->Delete(); |
| 486 |
// th->Clear(); |
| 487 |
// th->Reset(); |
| 488 |
// } else { |
| 489 |
// th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.); |
| 490 |
// th->SetMarkerStyle(20); |
| 491 |
// }; |
| 492 |
// |
| 493 |
tc->cd(); |
| 494 |
TString thid2 = Form("hCaloBragg2"); |
| 495 |
TH2F *th2 = dynamic_cast<TH2F*>(gDirectory->FindObject(thid2)); |
| 496 |
if ( th2 ) th2->Delete(); |
| 497 |
th2 = new TH2F(thid2,thid2,300,-0.5,300.,1000,0.,25.); //150 |
| 498 |
th2->SetMarkerStyle(20); |
| 499 |
th2->SetMarkerColor(kRed); |
| 500 |
// |
| 501 |
TString thid3 = Form("hCaloBragg3"); |
| 502 |
TH2F *th3 = dynamic_cast<TH2F*>(gDirectory->FindObject(thid3)); |
| 503 |
if ( th3 ) th3->Delete(); |
| 504 |
th3 = new TH2F(thid3,thid3,300,-0.5,300.,1000,0.,25.);//150. |
| 505 |
th3->SetMarkerStyle(20); |
| 506 |
th3->SetMarkerColor(kBlue); |
| 507 |
|
| 508 |
|
| 509 |
// tc->cd(1); |
| 510 |
// |
| 511 |
// for(Int_t i=0;i<=estremi[1][0];i++)th->Fill(Depth[i],dEpianimean[i]); |
| 512 |
for(Int_t i=0;i<=estremi[1][0];i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP); |
| 513 |
// th->Draw(); |
| 514 |
th2->Draw("same"); |
| 515 |
|
| 516 |
// tc->cd(2); |
| 517 |
tc->cd(); |
| 518 |
// |
| 519 |
for(Int_t i=0;i<=estremi[1][0];i++){ |
| 520 |
th3->Fill(Depth[i],dEpianiloop[i]); |
| 521 |
// printf(" i %i Depth %f depianiloop %f \n",i,Depth[i],dEpianiloop[i]); |
| 522 |
} |
| 523 |
th3->Draw(); |
| 524 |
th2->Draw("same"); |
| 525 |
|
| 526 |
tc->Modified(); |
| 527 |
tc->Update(); |
| 528 |
|
| 529 |
// |
| 530 |
gStyle->SetLabelSize(0); |
| 531 |
gStyle->SetNdivisions(1,"XY"); |
| 532 |
// |
| 533 |
}; |
| 534 |
|
| 535 |
|
| 536 |
|
| 537 |
void CaloBragg::LoadParam(){ |
| 538 |
|
| 539 |
// elem[Z-1][isotop] 0 is the most common one |
| 540 |
// |
| 541 |
|
| 542 |
elem[0][0] = 1.00782; //H 1 |
| 543 |
elem[0][1] = 2.01410; // 2H (Isotope) |
| 544 |
elem[0][2] = -1.; |
| 545 |
elem[0][3] = -1.; |
| 546 |
elem[0][4] = -1.; |
| 547 |
elem[0][5] = -1.; |
| 548 |
elem[0][6] = -1.; |
| 549 |
|
| 550 |
elem[1][0] = 4.002602; //He 2 |
| 551 |
elem[1][1] = 3.016029; // 3He (Isotope) |
| 552 |
elem[1][2] = -1.; |
| 553 |
elem[1][3] = -1.; |
| 554 |
elem[1][4] = -1.; |
| 555 |
elem[1][5] = -1.; |
| 556 |
elem[1][6] = -1.; |
| 557 |
|
| 558 |
elem[2][0] = 7.016004; //Li 3 |
| 559 |
elem[2][1] = 6.015123; //6Li (Isotope) |
| 560 |
elem[2][2] = -1.; |
| 561 |
elem[2][3] = -1.; |
| 562 |
elem[2][4] = -1.; |
| 563 |
elem[2][5] = -1.; |
| 564 |
elem[2][6] = -1.; |
| 565 |
|
| 566 |
elem[3][0] = 9.012182; //Be 4 |
| 567 |
elem[3][1] = 10.01353; //10Be (Isotope) (most stable) |
| 568 |
elem[3][2] = 7.01693; //9Be no EC in space? |
| 569 |
elem[3][3] = -1.; |
| 570 |
elem[3][4] = -1.; |
| 571 |
elem[3][5] = -1.; |
| 572 |
elem[3][6] = -1.; |
| 573 |
|
| 574 |
elem[4][0] = 11.00930; //B 5 |
| 575 |
elem[4][1] = 10.01294; //10B (Isotope) |
| 576 |
elem[4][2] = -1.; |
| 577 |
elem[4][3] = -1.; |
| 578 |
elem[4][4] = -1.; |
| 579 |
elem[4][5] = -1.; |
| 580 |
elem[4][5] = -1.; |
| 581 |
|
| 582 |
elem[5][0] = 12.0107; //C 6 |
| 583 |
elem[5][1] = 13.00335; //13C (Isotope) |
| 584 |
elem[5][2] = -1.; |
| 585 |
elem[5][3] = -1.; |
| 586 |
elem[5][4] = -1.; |
| 587 |
elem[5][5] = -1.; |
| 588 |
elem[5][5] = -1.; |
| 589 |
|
| 590 |
elem[6][0] = 14.00674; //N 7 |
| 591 |
elem[6][1] = 15.00011; //15N (Isotope) |
| 592 |
elem[6][2] = -1.; |
| 593 |
elem[6][3] = -1.; |
| 594 |
elem[6][4] = -1.; |
| 595 |
elem[6][5] = -1.; |
| 596 |
elem[6][5] = -1.; |
| 597 |
|
| 598 |
elem[7][0] = 15.99491; //O 8 |
| 599 |
elem[7][1] = 17.99916; //18O (Isotope) |
| 600 |
elem[7][2] = 16.99916; //17O (Isotope) |
| 601 |
elem[7][3] = -1.; |
| 602 |
elem[7][4] = -1.; |
| 603 |
elem[7][5] = -1.; |
| 604 |
elem[7][5] = -1.; |
| 605 |
|
| 606 |
elem[8][0] = 18.99840; //F 9 |
| 607 |
elem[8][1] = -1.; |
| 608 |
elem[8][2] = -1.; |
| 609 |
elem[8][3] = -1.; |
| 610 |
elem[8][4] = -1.; |
| 611 |
elem[8][5] = -1.; |
| 612 |
elem[8][5] = -1.; |
| 613 |
|
| 614 |
elem[9][0] = 19.99244; //Ne 10 |
| 615 |
elem[9][1] = 21.99138; //22Ne (Isotope) |
| 616 |
elem[9][2] = 20.99384; //21Ne 10 |
| 617 |
elem[9][3] = -1.; |
| 618 |
elem[9][4] = -1.; |
| 619 |
elem[9][5] = -1.; |
| 620 |
elem[9][6] = -1.; |
| 621 |
|
| 622 |
elem[10][0] = 22.98977; //Na 11 |
| 623 |
elem[10][1] = 21.99444; //22Na (Isotope) (most stable) |
| 624 |
elem[10][2] = -1.; |
| 625 |
elem[10][3] = -1.; |
| 626 |
elem[10][4] = -1.; |
| 627 |
elem[10][5] = -1.; |
| 628 |
elem[10][6] = -1.; |
| 629 |
|
| 630 |
elem[11][0] = 23.98504; //Mg 12 |
| 631 |
elem[11][1] = 25.98259; //26Mg (Isotope) |
| 632 |
elem[11][2] = 24.98504; //25Mg (Isotope) |
| 633 |
elem[11][3] = -1.; |
| 634 |
elem[11][4] = -1.; |
| 635 |
elem[11][5] = -1.; |
| 636 |
elem[11][6] = -1.; |
| 637 |
|
| 638 |
elem[12][0] = 26.98154; //Al 13 |
| 639 |
elem[12][1] = 25.98489; //26Al (Isotope) (most stable) |
| 640 |
elem[12][2] = -1.; |
| 641 |
elem[12][3] = -1.; |
| 642 |
elem[12][4] = -1.; |
| 643 |
elem[12][5] = -1.; |
| 644 |
elem[12][6] = -1.; |
| 645 |
|
| 646 |
elem[13][0] = 27.97692; //Si 14 |
| 647 |
elem[13][1] = 28.97649; //29Si (Isotope) |
| 648 |
elem[13][2] = 29.97377; //30Si (Isotope) |
| 649 |
elem[13][3] = -1.; |
| 650 |
elem[13][4] = -1.; |
| 651 |
elem[13][5] = -1.; |
| 652 |
elem[13][6] = -1.; |
| 653 |
|
| 654 |
elem[14][0] = 30.97376; //P 15 |
| 655 |
elem[14][1] = -1.; |
| 656 |
elem[14][2] = -1.; |
| 657 |
elem[14][3] = -1.; |
| 658 |
elem[14][4] = -1.; |
| 659 |
elem[14][5] = -1.; |
| 660 |
elem[14][6] = -1.; |
| 661 |
|
| 662 |
elem[15][0] = 31.97207; //S 16 |
| 663 |
elem[15][1] = 33.96787; //34S (Isotope) |
| 664 |
elem[15][2] = 32.97146; //33S (Isotope) |
| 665 |
elem[15][3] = 35.96708; //36S (Isotope) |
| 666 |
elem[15][4] = -1.; |
| 667 |
elem[15][5] = -1.; |
| 668 |
elem[15][6] = -1.; |
| 669 |
|
| 670 |
elem[16][0] = 34.96885; //Cl 17 |
| 671 |
elem[16][1] = 36.96831; //37Cl 17 |
| 672 |
elem[16][2] = 35.96890; //36Cl (Isotope) |
| 673 |
elem[16][3] = -1.; |
| 674 |
elem[16][4] = -1.; |
| 675 |
elem[16][5] = -1.; |
| 676 |
elem[16][6] = -1.; |
| 677 |
|
| 678 |
elem[17][0] = 39.962383; //Ar 18 |
| 679 |
elem[17][1] = 35.967545; //36Ar (Isotope) |
| 680 |
elem[17][2] = 37.962732; //38Ar (Isotope) |
| 681 |
elem[17][3] = 38.964313; //39Ar (Isotope) |
| 682 |
elem[17][4] = -1.; |
| 683 |
elem[17][5] = -1.; |
| 684 |
elem[17][6] = -1.; |
| 685 |
|
| 686 |
elem[18][0] = 38.963707; //K 19 |
| 687 |
elem[18][1] = 40.961825; //41K (Isotope) |
| 688 |
elem[18][2] = 39.963998; //40K (Isotope) |
| 689 |
elem[18][3] = -1.; |
| 690 |
elem[18][4] = -1.; |
| 691 |
elem[18][5] = -1.; |
| 692 |
elem[18][6] = -1.; |
| 693 |
|
| 694 |
elem[19][0] = 39.962590; //Ca 20 |
| 695 |
elem[19][1] = 43.955482; //44Ca (Isotope) |
| 696 |
elem[19][2] = 41.958618; //42Ca (Isotope) |
| 697 |
elem[19][3] = 42.958767; //43Ca (Isotope) |
| 698 |
elem[19][4] = 45.953693; //46Ca (Isotope) |
| 699 |
elem[19][5] = 40.962278; //41Ca (Isotope) |
| 700 |
elem[19][6] = -1.; |
| 701 |
|
| 702 |
elem[20][0] = 44.955912;//Sc 21 |
| 703 |
elem[20][1] = -1.; |
| 704 |
elem[20][2] = -1.; |
| 705 |
elem[20][3] = -1.; |
| 706 |
elem[20][4] = -1.; |
| 707 |
elem[20][5] = -1.; |
| 708 |
elem[20][6] = -1.; |
| 709 |
|
| 710 |
elem[21][0] = 47.947946; //Ti 22 |
| 711 |
elem[21][1] = 45.952632; //46Ti (Isotope) |
| 712 |
elem[21][2] = 46.951763; //47Ti (Isotope) |
| 713 |
elem[21][3] = 48.947870; //49Ti (Isotope) |
| 714 |
elem[21][4] = 49.944791; //50Ti (Isotope) |
| 715 |
elem[21][5] = 43.959690; //44Ti (Isotope) (half life 60y) |
| 716 |
elem[21][6] = -1.; |
| 717 |
|
| 718 |
elem[22][0] = 50.943960; //V 23 |
| 719 |
elem[22][1] = 49.947158; //50V (Isotope) |
| 720 |
elem[22][2] = -1.; |
| 721 |
elem[22][3] = -1.; |
| 722 |
elem[22][4] = -1.; |
| 723 |
elem[22][5] = -1.; |
| 724 |
elem[22][6] = -1.; |
| 725 |
|
| 726 |
elem[23][0] = 51.940507; //Cr 24 |
| 727 |
elem[23][1] = 52.940649; //53Cr (Isotope) |
| 728 |
elem[23][2] = 49.946044; //50Cr (Isotope) |
| 729 |
elem[23][3] = 53.938880; //54Cr (Isotope) |
| 730 |
elem[23][4] = -1.; |
| 731 |
elem[23][5] = -1.; |
| 732 |
elem[23][6] = -1.; |
| 733 |
|
| 734 |
elem[24][0] = 54.938049;//Mn 25 |
| 735 |
elem[24][1] = 52.941290;//53Mn (Isotope) |
| 736 |
elem[24][2] = -1.; |
| 737 |
elem[24][3] = -1.; |
| 738 |
elem[24][4] = -1.; |
| 739 |
elem[24][5] = -1.; |
| 740 |
elem[24][6] = -1.; |
| 741 |
|
| 742 |
elem[25][0] = 55.934937; //Fe 26 |
| 743 |
elem[25][1] = 53.939610; //54Fe (Isotope) |
| 744 |
elem[25][2] = 56.935394; //57Fe (Isotope) |
| 745 |
elem[25][3] = 57.933276; //58Fe (Isotope) |
| 746 |
elem[25][4] = 59.934072; //58Fe (Isotope) |
| 747 |
|
| 748 |
elem[26][0] = 58.933195; //Co 27 |
| 749 |
elem[26][1] = 59.933817; //60Co (Isotope) |
| 750 |
elem[26][2] = -1.; |
| 751 |
elem[26][3] = -1.; |
| 752 |
elem[26][4] = -1.; |
| 753 |
elem[26][5] = -1.; |
| 754 |
elem[26][6] = -1.; |
| 755 |
|
| 756 |
|
| 757 |
elem[27][0] = 57.935343; //Ni 28 |
| 758 |
elem[27][1] = 61.928345; //62Ni (Isotope) |
| 759 |
elem[27][2] = 59.930786; //60Ni (Isotope) |
| 760 |
elem[27][3] = 60.931056; //61Ni (Isotope) |
| 761 |
elem[27][4] = 63.927966; //64Ni (Isotope) |
| 762 |
elem[27][5] = 58.934346; //59Ni (Isotope) |
| 763 |
elem[27][6] = -1.; |
| 764 |
|
| 765 |
elem[28][0] = 62.929597; //Cu 29 |
| 766 |
elem[28][1] = 64.927789; //65Cu (Isotope) |
| 767 |
elem[28][2] = -1.; |
| 768 |
elem[28][3] = -1.; |
| 769 |
elem[28][4] = -1.; |
| 770 |
elem[28][5] = -1.; |
| 771 |
elem[28][6] = -1.; |
| 772 |
|
| 773 |
elem[29][0] = 63.929142; //Zn 30 |
| 774 |
elem[29][1] = 65.926033; //66Zn (Isotope) |
| 775 |
elem[29][2] = 67.924844; //68Zn (Isotope) |
| 776 |
elem[29][3] = 66.927127; //67Zn (Isotope) |
| 777 |
elem[29][4] = 69.925319; //70Zn (Isotope) |
| 778 |
elem[29][5] = -1.; |
| 779 |
elem[29][6] = -1.; |
| 780 |
|
| 781 |
elem[30][0] = 68.925573; //Ga 31 |
| 782 |
elem[30][1] = 70.924701; //71Ga (Isotope) |
| 783 |
elem[30][2] = -1.; |
| 784 |
elem[30][3] = -1.; |
| 785 |
elem[30][4] = -1.; |
| 786 |
elem[30][5] = -1.; |
| 787 |
elem[30][6] = -1.; |
| 788 |
|
| 789 |
elem[31][0] = 73.921177; //Ge 32 |
| 790 |
elem[31][1] = 71.922075; //72Ge (Isotope) |
| 791 |
elem[31][2] = 69.924247; //70Ge (Isotope) |
| 792 |
elem[31][3] = 75.921403; //76Ge (Isotope) |
| 793 |
elem[31][4] = 73.923459; //73Ge (Isotope) |
| 794 |
elem[31][5] = -1.; |
| 795 |
elem[31][6] = -1.; |
| 796 |
|
| 797 |
|
| 798 |
//parametri calorimetro |
| 799 |
NPLA = 22; |
| 800 |
NCHA = 96; |
| 801 |
nView = 2; |
| 802 |
|
| 803 |
AA = 0.96;//mm larghezza strip |
| 804 |
ADIST = 80.5;//mm distanza tra pad |
| 805 |
PIANO = 8.59;//mm distanza |
| 806 |
|
| 807 |
ySi = 0.38;//mm spessore silicio |
| 808 |
yW = 2.66;//mm spessore tungsteno |
| 809 |
rhoSi = 2.33;//g/cm3 densita' silicio |
| 810 |
rhoW = 19.3;//g/cm3 densita' tugsteno |
| 811 |
MIP = 0.106;//Mev g/cm2 energia al minimo nel silicio per 0.38 mm |
| 812 |
|
| 813 |
emin = 0.; |
| 814 |
|
| 815 |
//parametri bethe-bloch |
| 816 |
pigr = 3.1415; |
| 817 |
Na = 6.02e-23; |
| 818 |
ZA = 0.49; /*Z/A per Si*/ |
| 819 |
// ISi =182e-06; /*MeV*/ |
| 820 |
ISi = 171e-06; /*MeV*/ |
| 821 |
IW = 735e-06; /*MeV*/ |
| 822 |
// ISi =0.0001059994; /*GeV!!*/ no era giusto!! |
| 823 |
Me = 0.511; /* MeV*/ |
| 824 |
MassP = 931.27;/*MeV*/ |
| 825 |
r2 = 7.95e-26; /*ro*ro in cm */ |
| 826 |
|
| 827 |
}; |
| 828 |
|
| 829 |
|
| 830 |
|
| 831 |
// |
| 832 |
void CaloBragg::conversione(){ |
| 833 |
|
| 834 |
// calcolo spessore Si attraverato in funzione dell'inclinazione |
| 835 |
// e conversione dello spessore di W in Si e correzione del valore |
| 836 |
// della Mip pe lo spessore effettivo |
| 837 |
// |
| 838 |
// in : evento |
| 839 |
// |
| 840 |
// out: out[0] = gcm2Si = spessore silicio attraversato nel piano |
| 841 |
// out[1] = WinSi = spessore equivalente in Si del W attraversato |
| 842 |
// out[2] = Mip = fattore conversione energia riscalato allo spessore attrversatonel piano |
| 843 |
|
| 844 |
Float_t SiCross=0.; |
| 845 |
Float_t WCross = 0.; |
| 846 |
Float_t ytgx = 0; |
| 847 |
Float_t ytgy = 0; |
| 848 |
Float_t a = 0.; |
| 849 |
|
| 850 |
/*silicio*/ |
| 851 |
ytgx = ySi * L2->GetCaloLevel2()->tanx[0]; |
| 852 |
ytgy = ySi * L2->GetCaloLevel2()->tany[0]; |
| 853 |
|
| 854 |
//lunghezza effettiva di silicio attraversata (mm) |
| 855 |
SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy)); |
| 856 |
|
| 857 |
spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2 |
| 858 |
|
| 859 |
/*tungsteno*/ |
| 860 |
ytgx = yW * L2->GetCaloLevel2()->tanx[0]; |
| 861 |
ytgy = yW * L2->GetCaloLevel2()->tany[0]; |
| 862 |
|
| 863 |
//rapporto tra rilasci energetici nei due materiali |
| 864 |
WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad |
| 865 |
//gcm2W = WCross/10. * rhoW; |
| 866 |
|
| 867 |
// (g/cm2W)/(g/cm2Si) |
| 868 |
spessore[3] = (WCross/10.) * rhoW; |
| 869 |
a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664); //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664); |
| 870 |
spessore[1] = a; |
| 871 |
//riscala mip allo spessore attraversato |
| 872 |
spessore[2] = MIP*(SiCross/ySi); |
| 873 |
};//end conversione |
| 874 |
|
| 875 |
|
| 876 |
|
| 877 |
|
| 878 |
|
| 879 |
void CaloBragg::BetheBloch(Float_t *x, Float_t *z, Float_t *Mass, Float_t *gam, Float_t *Bet, Float_t *out, Float_t II){ |
| 880 |
|
| 881 |
//rilascio energetico con bethe bloch con correzioni |
| 882 |
//in: x: g/cm2 |
| 883 |
// z: carica |
| 884 |
// Mass: Massa uma |
| 885 |
// Ene: energia particella MeV//tolta |
| 886 |
// gam: (etot/massa) |
| 887 |
// Bet: rad((g2-1)/g2) |
| 888 |
// |
| 889 |
//out: energia rilasciata MeV |
| 890 |
|
| 891 |
|
| 892 |
Float_t eta =0.; |
| 893 |
Float_t Wmax =0.; |
| 894 |
Float_t lg =0.; |
| 895 |
Float_t Energia=0.; |
| 896 |
Float_t C=0.; |
| 897 |
Float_t INo = ISi; |
| 898 |
|
| 899 |
if ( usenewBB ) INo = II; |
| 900 |
|
| 901 |
eta = (*gam)*(*Bet); |
| 902 |
|
| 903 |
//Bet=3/gam; SQ(*gam) * SQ(*Bet) |
| 904 |
Wmax = 2.* Me * SQ(eta) / (1. + 2.*(*gam)*Me/(*Mass) + SQ(Me)/SQ(*Mass)); |
| 905 |
|
| 906 |
lg = 2.* Me * SQ(eta) * Wmax / SQ(INo); |
| 907 |
// Energia = x* 2 * pigr * Na * r2 * Me * rhoSi *ZA* SQ(z)/SQ(Bet) * lg; |
| 908 |
C=(0.42237*pow(eta,-2.) + 0.0304*pow(eta,-4.) - 0.00038*pow(eta,-6.))*pow(10.,-6.)* pow(INo,2.) + |
| 909 |
(3.858*pow(eta,-2.) - 0.1668*pow(eta,-4.) + 0.00158*pow(eta,-6.))*pow(10.,-9.)*pow(INo,3.); |
| 910 |
|
| 911 |
if(eta <= 0.13) C= C * log(eta/0.0653)/log(0.13/0.0653); |
| 912 |
|
| 913 |
Energia = (*x) * 0.307/28.09 * 14. *SQ(*z)/SQ(*Bet)*(0.5*log(lg) - SQ(*Bet) - C/14.); |
| 914 |
|
| 915 |
*out =Energia;//out |
| 916 |
|
| 917 |
};//end Bethebloch |
| 918 |
|
| 919 |
|
| 920 |
|
| 921 |
|
| 922 |
void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Int_t *isotope, Float_t *Etot, Float_t *out, Float_t II){ |
| 923 |
|
| 924 |
/*perdita di energia per ioni pesanti (come da routine geant)*/ |
| 925 |
// in : dx => spessore g/cm2 |
| 926 |
// Z => carica |
| 927 |
// Etot => energia perticella |
| 928 |
// |
| 929 |
// out: energia persa |
| 930 |
|
| 931 |
|
| 932 |
Float_t Q=0.; |
| 933 |
Float_t v=0.; |
| 934 |
Float_t gam=0.; |
| 935 |
Float_t Bet=0.; |
| 936 |
Float_t dEP=0.; |
| 937 |
|
| 938 |
// gamma // Mass = A * MassP; /*in Mev/c2*/ |
| 939 |
gam = (*Etot)/(elem[*Z-1][*isotope]*MassP); // E = gamma M c2 |
| 940 |
|
| 941 |
|
| 942 |
Bet = sqrt((SQ(gam) -1.)/SQ(gam)); |
| 943 |
|
| 944 |
// v= 121.4139*(Bet/pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet/pow((*Z),(2./3.)))); |
| 945 |
v= 121.4139*(Bet*pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet*pow((*Z),(2./3.)))); // EMI AAAAGGH!! |
| 946 |
|
| 947 |
//carica effettiva |
| 948 |
Q= (*Z)*(1- (1.034 - 0.1777*exp(-0.08114*(*Z)))*exp(-v)); |
| 949 |
|
| 950 |
//perdita energia per un protone |
| 951 |
Float_t protone =1.; |
| 952 |
// Float_t Mass=(elem[*Z-1]*MassP); //EMI |
| 953 |
// BetheBloch(dx, &protone, &Mass, &gam, &Bet, &dEP);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP); |
| 954 |
|
| 955 |
BetheBloch(dx, &protone, &MassP, &gam, &Bet, &dEP, II);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP); //EMI |
| 956 |
|
| 957 |
*out= (SQ(Q)*(dEP));//*dx; |
| 958 |
|
| 959 |
|
| 960 |
};//end ELOSS |
| 961 |
|
| 962 |
|
| 963 |
|
| 964 |
|
| 965 |
void CaloBragg::Enetrack(Int_t* Z, Int_t* isotope, Float_t* E0, Float_t* primo,Float_t* ultimo, Float_t out[]){ |
| 966 |
|
| 967 |
//calcola energia rilasciata sulla traccia (usa ELOSS) |
| 968 |
// in : Z =>carica |
| 969 |
// E0 =>energia |
| 970 |
// spess2[3] => conversione spessore Si, Si in W, mip |
| 971 |
// primo => posizione primo piano attraversato |
| 972 |
// |
| 973 |
// out: array[44] =>rilasci energetici calcolati per ogni piano[44] dopo il primo(estremi[0][0]) |
| 974 |
|
| 975 |
|
| 976 |
|
| 977 |
Float_t dE=0.; //energia rilasciata |
| 978 |
Float_t Ezero= *E0;//energia iniziale |
| 979 |
|
| 980 |
//azzero energia rilasciata sui piani |
| 981 |
memset(out, 0, 2*NPLA*sizeof(Float_t)); |
| 982 |
|
| 983 |
Float_t Massa = (elem[(*Z)-1][*isotope] * MassP); |
| 984 |
|
| 985 |
for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){ |
| 986 |
dE=0.; |
| 987 |
//spessore silicio corretto x inclinazione, z, energia, out:rilascio |
| 988 |
ELOSS(&spessore[0], Z , isotope , &Ezero, &dE, ISi);//spessore in g/cm2!! |
| 989 |
|
| 990 |
if(dE!=dE) return; //controlla che non sia un NaN |
| 991 |
|
| 992 |
if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop |
| 993 |
out[ipla] = Ezero - Massa; //MeV |
| 994 |
return; |
| 995 |
|
| 996 |
}else{ |
| 997 |
out[ipla] = dE; //MeV |
| 998 |
Ezero = Ezero - dE;//energia residua |
| 999 |
// if ( debug ) printf(" zompa %i out %f dE %f ezero %f \n",ipla,out[ipla],dE,Ezero); |
| 1000 |
}; |
| 1001 |
//se sono su un piano Y (tutti i pari) dopo c'e' il tungsteno |
| 1002 |
if(ipla%2 == 0){ |
| 1003 |
/*tungsteno*/ |
| 1004 |
dE=0.; |
| 1005 |
Float_t sp = 0.; |
| 1006 |
Float_t II = ISi; |
| 1007 |
if ( usenewBB ){ |
| 1008 |
sp = spessore[3]; |
| 1009 |
II = IW; |
| 1010 |
} else { |
| 1011 |
sp = spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato in g/cm2 |
| 1012 |
} |
| 1013 |
// printf(" sp %f II %f \n",sp,II); |
| 1014 |
ELOSS(&sp, Z, isotope , &Ezero, &dE,II); |
| 1015 |
if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop |
| 1016 |
return; |
| 1017 |
}else{ |
| 1018 |
Ezero = Ezero -dE;//energia residua |
| 1019 |
}; |
| 1020 |
}; |
| 1021 |
|
| 1022 |
};//fine loop piani |
| 1023 |
|
| 1024 |
|
| 1025 |
};//end Enetrack |
| 1026 |
|
| 1027 |
|
| 1028 |
|
| 1029 |
void CaloBragg::chiquadro(Float_t dE[], Float_t out[]){ |
| 1030 |
|
| 1031 |
// calcola chi2 tra energia calcolata e misurata |
| 1032 |
// in : dE[44] =>energia calcolata |
| 1033 |
// calo3[44][2]=> [0]strip attraversata [1]energia misurata per ogni piano |
| 1034 |
// estr2 => array con primo[0][0] e ultimo[1][0] piano attraversati ed energie[][1] |
| 1035 |
// |
| 1036 |
// out: array[3]=> (chi2; piani scartati consecutivi(79= >3 quindi frammentato); piani scartati totale) |
| 1037 |
|
| 1038 |
|
| 1039 |
Float_t sum = 0.; |
| 1040 |
Float_t PianoPrecedente=0.; |
| 1041 |
Float_t badplane=0.; |
| 1042 |
Float_t badplanetot=0.; |
| 1043 |
Float_t w,wi; |
| 1044 |
// |
| 1045 |
if ( newchi2 ){ |
| 1046 |
ndf = 0; |
| 1047 |
sum = 0.; |
| 1048 |
for( Int_t ipla=((int)(estremi[0][0])); ipla<= ((int)(estremi[1][0])); ipla++){ |
| 1049 |
sum += pow((dE[ipla] - (calorimetro[ipla][1] * spessore[2]))/(0.05*dE[ipla]),2.); |
| 1050 |
// printf(" quiqui: dE %f calor %f spessore[2] %f \n",dE[ipla],spessore[2]*calorimetro[ipla][1],spessore[2]); |
| 1051 |
ndf++; |
| 1052 |
} |
| 1053 |
ndf -= 2; |
| 1054 |
if ( ndf > 0 ) sum /= (float)ndf; |
| 1055 |
out[0] = sum; |
| 1056 |
out[1] = 0.; |
| 1057 |
out[2] = (int)(estremi[1][0])-ndf; |
| 1058 |
// printf(" sum %f ndf %i \n ",sum,ndf); |
| 1059 |
} else { |
| 1060 |
for(Int_t ipla=0; ipla<2*NPLA; ipla++){ |
| 1061 |
//tutti i piani attraversati dalla traiettoria |
| 1062 |
if(calorimetro[ipla][0] != -1.){ // |
| 1063 |
w=0.; //normalizzazione; |
| 1064 |
wi=1.;//peso |
| 1065 |
|
| 1066 |
//tolgo piani attraversati dalla traccia ma precedenti il piano individuato come ingresso |
| 1067 |
if (ipla<estremi[0][0]) wi=0.; |
| 1068 |
|
| 1069 |
//tolgo piani attraversati da traccia ma successivi all'ultimo se sono diversi da 0 |
| 1070 |
//if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.; |
| 1071 |
if((ipla>estremi[1][0])) wi=0.; |
| 1072 |
|
| 1073 |
//normalizzazione |
| 1074 |
if (calorimetro[ipla][1] != 0.) w=1./(calorimetro[ipla][1]* MIP); // |
| 1075 |
|
| 1076 |
//tolgo piani con rilasci inferiori al 30% del precedente |
| 1077 |
if(calorimetro[ipla][1] < (0.7*PianoPrecedente)){ // cosi' i piani senza rilascio non vengono considerati nel calcolo del chi2 |
| 1078 |
wi=0.; |
| 1079 |
//se sono piani intermedi (non si e' fermta) li considero non buoni |
| 1080 |
if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=0.)){// |
| 1081 |
badplane+=1.; |
| 1082 |
badplanetot+=1.; |
| 1083 |
}; |
| 1084 |
}; |
| 1085 |
|
| 1086 |
//meno peso ai piani con rilasci maggiori di 1000 MIP |
| 1087 |
// if(calorimetro[ipla][1] > 1000) wi=0.5; |
| 1088 |
if(calorimetro[ipla][1] > 1200.) wi=0.5; |
| 1089 |
if(debug) printf("chiquadro start \n "); |
| 1090 |
Float_t arg = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP)); |
| 1091 |
|
| 1092 |
sum += SQ(arg); // w*wi*(dEpiani[p][v]-(eplane[p][v]*MIP))));//( dEpiani[p][v] - (eplane[p][v]*MIP)); |
| 1093 |
if(debug){ |
| 1094 |
printf("dedx calcolata %f e reale %f \n",dE[ipla],(calorimetro[ipla][1] * MIP)); |
| 1095 |
} |
| 1096 |
//se trovo piano non buono (tolto quindi wi=0) non modifico il piano precedente |
| 1097 |
if(wi != 0.){// |
| 1098 |
PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente |
| 1099 |
badplane = 0.;//azzero contatore piani scartati consecutivi |
| 1100 |
}; |
| 1101 |
}; |
| 1102 |
|
| 1103 |
//da Emi |
| 1104 |
if(badplane > 2){ |
| 1105 |
// printf(" AAAAAAAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG\n"); |
| 1106 |
out[1] =79.; |
| 1107 |
break; |
| 1108 |
}; |
| 1109 |
|
| 1110 |
};//fine loop piani |
| 1111 |
//chi2,frammentato,pskip |
| 1112 |
out[0]=sum; |
| 1113 |
out[2]=badplanetot; |
| 1114 |
} |
| 1115 |
};//end chiquadro |
| 1116 |
|
| 1117 |
|
| 1118 |
|
| 1119 |
void CaloBragg::loopze( Float_t step, Float_t E0,Float_t Zstart, Float_t Zlimite, Int_t nostep = 1000){ |
| 1120 |
// |
| 1121 |
//loop su z ed energie per trovare miglior z (ed energia) |
| 1122 |
//in: nloop => energia massima da provare (nloop x E0) |
| 1123 |
// E0 => energia iniziale (intergale) |
| 1124 |
// Zstart => minimo z da cui patire |
| 1125 |
// Zlimite => z a cui fermarsi (z al minimo di ionizz sul 1o piano) |
| 1126 |
// |
| 1127 |
//out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2 |
| 1128 |
// |
| 1129 |
|
| 1130 |
//printf("entrato"); |
| 1131 |
memset(dEplan,0,2*NPLA*sizeof(Float_t)); |
| 1132 |
|
| 1133 |
Int_t Z = 0;// z iniziale |
| 1134 |
|
| 1135 |
Int_t isotope=0; |
| 1136 |
|
| 1137 |
Float_t Massa = 0.; |
| 1138 |
|
| 1139 |
Float_t Stepint =(step)/(Float_t)nostep;//passo per il calcolo di energia |
| 1140 |
|
| 1141 |
Float_t energia =0.;//energia del loop |
| 1142 |
|
| 1143 |
Float_t chi2[3] = {0,0,0};//out dal calcolo chi2: chi2, piani consecutivi saltati, piani totali saltati |
| 1144 |
|
| 1145 |
Int_t zmin = (int)Zstart; |
| 1146 |
Int_t max=32;//max z di cui so la massa :P |
| 1147 |
if((Zlimite)<=31) max=(int)(Zlimite) + 1; |
| 1148 |
|
| 1149 |
if(debug) printf("loopze inizio max %d \n",max); |
| 1150 |
if ( fzeta > 0. ){ |
| 1151 |
zmin = fzeta; |
| 1152 |
max = fzeta+1; |
| 1153 |
} |
| 1154 |
|
| 1155 |
Int_t colmax=32; |
| 1156 |
Int_t rowmax=3000; |
| 1157 |
Int_t isomax=7; |
| 1158 |
|
| 1159 |
Float_t matrixchi2[colmax][isomax][rowmax][3]; |
| 1160 |
memset(matrixchi2, 0, colmax*isomax*rowmax*3*sizeof(Float_t)); |
| 1161 |
|
| 1162 |
Int_t imin = 1-nostep/2; |
| 1163 |
Int_t imax = nostep/2; |
| 1164 |
|
| 1165 |
//loop elementi |
| 1166 |
for(Int_t inucl=zmin; inucl<max; inucl++){ |
| 1167 |
|
| 1168 |
Z= inucl; |
| 1169 |
|
| 1170 |
//loop isotopi |
| 1171 |
while ( elem[inucl-1][isotope] > 0. ){ |
| 1172 |
|
| 1173 |
if( fiso != -1 ){ |
| 1174 |
isotope=fiso; |
| 1175 |
if(debug) printf("In Loopze - Isotope N %d",isotope); |
| 1176 |
} |
| 1177 |
Massa = elem[inucl-1][isotope]*MassP; |
| 1178 |
|
| 1179 |
//loop energia |
| 1180 |
Int_t iene2 = -1; |
| 1181 |
|
| 1182 |
// 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 <=?? |
| 1183 |
|
| 1184 |
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 <=?? |
| 1185 |
iene2++; |
| 1186 |
energia= Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica.. |
| 1187 |
|
| 1188 |
|
| 1189 |
if( fene > 0. ) energia=fene; //forza l'energia |
| 1190 |
if (debug) printf("loopze energia %f, z %d, isotopo %d ,iene %d\n",energia,Z,isotope,iene); |
| 1191 |
// printf(" energia %f , forzata %f \n",energia,fene); |
| 1192 |
Enetrack(&Z, &isotope, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani |
| 1193 |
|
| 1194 |
chiquadro(dEplan,chi2); //calcolo chi2 |
| 1195 |
if (debug) printf("loopze chi %f \n",chi2[0]); |
| 1196 |
if(debug && TMath::Finite(chi2[0])==1 && (TMath::IsNaN(chi2[0])!=1) ) printf("loopze fin mat %f \n",chi2[0]); |
| 1197 |
// printf(" last deplan from: Z = %i iene %i energia %f chi2 %f \n",inucl,iene,energia,chi2[0]); |
| 1198 |
if( (chi2[1] != 79.) ){//salto quelli che frammentano |
| 1199 |
matrixchi2[inucl][isotope][iene2][0]=chi2[0];//valore chi2 per questo z a questa energia |
| 1200 |
matrixchi2[inucl][isotope][iene2][1]=energia;//energia per questo chi2 |
| 1201 |
matrixchi2[inucl][isotope][iene2][2]=chi2[2];//piani saltati nel chi2 |
| 1202 |
if( fene > 0. ) break; |
| 1203 |
} else { |
| 1204 |
matrixchi2[inucl][isotope][iene2][0]=numeric_limits<Float_t>::max();//valore chi2 per questo z a questa energia |
| 1205 |
matrixchi2[inucl][isotope][iene2][1]=numeric_limits<Float_t>::max();//energia per questo chi2 |
| 1206 |
matrixchi2[inucl][isotope][iene2][2]=numeric_limits<Float_t>::max();//piani saltati nel chi2 |
| 1207 |
break; |
| 1208 |
} |
| 1209 |
|
| 1210 |
}//fine loop energia |
| 1211 |
|
| 1212 |
if( fiso != -1 ){ |
| 1213 |
if(debug) printf("exited form isotopes loop"); |
| 1214 |
break; |
| 1215 |
} |
| 1216 |
|
| 1217 |
isotope++; //incremento il contatore isotopi |
| 1218 |
}//fine loop isotopi |
| 1219 |
isotope=0; //riazzero il contatore isotopi |
| 1220 |
|
| 1221 |
}//fine loop z |
| 1222 |
|
| 1223 |
isotope=0;//non dovrebbe servire |
| 1224 |
|
| 1225 |
//Emi |
| 1226 |
for (Int_t nu=zmin; nu<max; nu++){ |
| 1227 |
if( fiso != -1 ){ |
| 1228 |
isotope=fiso; |
| 1229 |
if(debug) printf("In Loopze EMI - Isotope N %d",isotope); |
| 1230 |
} |
| 1231 |
while(elem[nu-1][isotope]> 0.){ |
| 1232 |
for (Int_t en=0; en<nostep; en++){ |
| 1233 |
if((matrixchi2[nu][isotope][en][0]<bestchi2[0]) && (matrixchi2[nu][isotope][en][0] >0.)){ |
| 1234 |
bestchi2[0]= matrixchi2[nu][isotope][en][0];// chi2 |
| 1235 |
bestchi2[1]= (Float_t)nu; // z |
| 1236 |
bestchi2[2]= matrixchi2[nu][isotope][en][1];//energia; |
| 1237 |
bestchi2[3]= matrixchi2[nu][isotope][en][2];// totale piani saltati |
| 1238 |
bestchi2[4]= (Float_t)isotope; //isotopo |
| 1239 |
} |
| 1240 |
} |
| 1241 |
|
| 1242 |
if( fiso != -1 ){ |
| 1243 |
if(debug) printf("exited form isotopes loop"); |
| 1244 |
break; |
| 1245 |
} |
| 1246 |
|
| 1247 |
isotope++; |
| 1248 |
} |
| 1249 |
isotope=0; |
| 1250 |
} |
| 1251 |
|
| 1252 |
};//endloopze |
| 1253 |
|
| 1254 |
|
| 1255 |
|
| 1256 |
|
| 1257 |
|
| 1258 |
// void CaloBragg::mediatroncata(){ |
| 1259 |
// //calcolo Z con media troncata e utilizzo questo Z per trovare l'energia migliore |
| 1260 |
// //in: ordplane[44] => array con energia dei piani |
| 1261 |
// // spess[3] => conversioni spessore di silicio, w, mip |
| 1262 |
// // estr[2][2] => primo[0][0] e ultimo[1][0] piano attraversati ed energie[][1] |
| 1263 |
// // calo[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano |
| 1264 |
// // integrale => energia totale nel calorimetro considerando il W |
| 1265 |
// // |
| 1266 |
// // out[4] chi2,z,Etot,Pskip |
| 1267 |
|
| 1268 |
// Float_t ordplane[44];//mi serve per la media troncata |
| 1269 |
// memset(ordplane,0,44*sizeof(Float_t)); |
| 1270 |
|
| 1271 |
// for(Int_t ipla=0; ipla< 2*NPLA; ipla++) ordplane[ipla]=calorimetro[ipla][1]; //energia del piano |
| 1272 |
|
| 1273 |
|
| 1274 |
// //ordino tutte le energie dei piani in ordine crescente |
| 1275 |
|
| 1276 |
// Long64_t work[200]; |
| 1277 |
// Int_t ind = 0; |
| 1278 |
// //Int_t l = 0; |
| 1279 |
// Int_t RN = 0; |
| 1280 |
// Float_t sum4 = 0.; |
| 1281 |
// Float_t qm = 0.; |
| 1282 |
// // |
| 1283 |
// //Float_t qmt = ethr*0.8; // *0.9 |
| 1284 |
// // |
| 1285 |
// //Int_t uplim = TMath::Max(3,N); |
| 1286 |
// // |
| 1287 |
// while ( RN < 4 && ind < 44 ){ |
| 1288 |
// qm = TMath::KOrdStat(44,ordplane,ind,work); |
| 1289 |
// if (qm >= 0.7 ){ |
| 1290 |
// if ( RN < 4 ){ |
| 1291 |
// sum4 += qm; |
| 1292 |
// RN++; |
| 1293 |
// }; |
| 1294 |
// // l++; |
| 1295 |
// // if ( debug ) printf(" value no %i qm %f sum4 %f \n",l,qm,sum4); |
| 1296 |
// }; |
| 1297 |
// ind++; |
| 1298 |
// }; |
| 1299 |
// // |
| 1300 |
// sum4 /= (Float_t)RN; |
| 1301 |
// Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));//ma non e'/1?? |
| 1302 |
// if(Zmean ==0.) Zmean=1.; |
| 1303 |
// if ( Zmean < 1. ) Zmean = 1.; |
| 1304 |
|
| 1305 |
|
| 1306 |
// // Zmean =round(Zmean); |
| 1307 |
// // if(Zmean <1.) Zmean=1.; |
| 1308 |
|
| 1309 |
// // if(Zmean >0.)Zmean =round(Zmean); |
| 1310 |
|
| 1311 |
// //======== per i nuclei======= |
| 1312 |
// if (Zmean >=2.){ |
| 1313 |
// ind = 0; |
| 1314 |
// RN = 0; |
| 1315 |
// sum4 = 0.; |
| 1316 |
// qm = 0.; |
| 1317 |
// while ( RN < 4 && ind < 44 ){ |
| 1318 |
// qm = TMath::KOrdStat(44,ordplane,ind,work); |
| 1319 |
// if (qm >= (Zmean*Zmean)-Zmean*Zmean*0.2 ){ |
| 1320 |
// if ( RN < 4 ){ |
| 1321 |
// sum4 += qm; |
| 1322 |
// RN++; |
| 1323 |
// }; |
| 1324 |
// }; |
| 1325 |
// ind++; |
| 1326 |
// }; |
| 1327 |
// // |
| 1328 |
// sum4 /= (Float_t)RN; |
| 1329 |
// Zmean = (sqrt((sum4*MIP)/(4.*spessore[2])));//ma non e' /1?? |
| 1330 |
// } |
| 1331 |
|
| 1332 |
|
| 1333 |
// //calcolo energia migliore per Z trovato con media troncata |
| 1334 |
// // Float_t zmin=Zmean; |
| 1335 |
// Float_t zmin=round(Zmean); |
| 1336 |
|
| 1337 |
// bestchi2[0]=10000.; |
| 1338 |
// bestchi2[1]=0.; |
| 1339 |
// bestchi2[2]=0.; |
| 1340 |
// bestchi2[3]=0.; |
| 1341 |
// Float_t zero=0.; |
| 1342 |
|
| 1343 |
// // step energia zstart zstop |
| 1344 |
// loopze(Integrale,zero,zmin,zmin); |
| 1345 |
|
| 1346 |
|
| 1347 |
// qtchi2=bestchi2[0]; |
| 1348 |
// qtz=bestchi2[1]; |
| 1349 |
// qtetot=bestchi2[2]; |
| 1350 |
// qtpskip=bestchi2[3]; |
| 1351 |
// };//end mediatroncata |
| 1352 |
|
| 1353 |
|
| 1354 |
|
| 1355 |
void CaloBragg::Zdaloop(){ |
| 1356 |
//calcolo Z con un loop su tutti i possibli Z ed energie |
| 1357 |
//in: ordplane[44]=> array con energia dei piani |
| 1358 |
// spess1[3]=> conversioni spessore di silicio, w e mip |
| 1359 |
// estr3[2][2]=> primo[0][0] e ultimo[1][0] piano ed energie |
| 1360 |
// calo1[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano |
| 1361 |
// integrale=> energia totale nel calorimetro considerando il W |
| 1362 |
// |
| 1363 |
// out[4] chi2,z,Etot,Pskip |
| 1364 |
|
| 1365 |
|
| 1366 |
/*z se particella fosse al minimo*/ //energia1piano/mip corretta |
| 1367 |
// Float_t zmax = round(sqrt(estremi[0][1]/spessore[2])); |
| 1368 |
// if(zmax<31)zmax=zmax+1; |
| 1369 |
|
| 1370 |
/*calcolo Z ed E con loop sui vari elementi ed energie*/ |
| 1371 |
|
| 1372 |
Float_t zmin=1.; |
| 1373 |
Float_t zmax=32.; |
| 1374 |
Float_t bestchitemp[5] = {0,0,0,0,0}; |
| 1375 |
|
| 1376 |
bestchi2[0]=numeric_limits<Float_t>::max(); |
| 1377 |
bestchi2[1]=0.; |
| 1378 |
bestchi2[2]=0.; |
| 1379 |
bestchi2[3]=0.; |
| 1380 |
bestchi2[4]=0.; |
| 1381 |
// Float_t zero=0.; |
| 1382 |
//------------primo loop ---------------------- |
| 1383 |
// energia ezero, zstart zstop |
| 1384 |
// loopze(Integrale,zero,zmin,zmax); |
| 1385 |
|
| 1386 |
//-> loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,50); |
| 1387 |
loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,200); |
| 1388 |
|
| 1389 |
// loopze(Integrale*2.,Integrale/100.,zmin,zmax); |
| 1390 |
if ( debug) printf("Zdaloop start Integrale %f , outene %f \n",Integrale,bestchi2[2]); |
| 1391 |
|
| 1392 |
//------------secondo loop ---------------------- |
| 1393 |
for(Int_t i=0;i<5;i++) bestchitemp[i]=bestchi2[i]; |
| 1394 |
bestchi2[0]=numeric_limits<Float_t>::max(); |
| 1395 |
bestchi2[1] = 0.; |
| 1396 |
bestchi2[2] = 0.; |
| 1397 |
bestchi2[3] = 0.; |
| 1398 |
bestchi2[4] = 0.;//riazzero |
| 1399 |
|
| 1400 |
Float_t step = bestchitemp[2];// |
| 1401 |
// zero=0.; // qualsiasi altro valore peggiora le cose |
| 1402 |
// zmin=zmax=bestchitemp[1]; |
| 1403 |
zmin=bestchitemp[1]-1; |
| 1404 |
zmax=bestchitemp[1]+1; |
| 1405 |
//loopze(step,zero,zmin,zmax); // |
| 1406 |
|
| 1407 |
//-> loopze(step,step/2.,zmin,zmax,200); // |
| 1408 |
loopze(step,step/2.,zmin,zmax,500); // |
| 1409 |
|
| 1410 |
//step = bestchitemp[2];// |
| 1411 |
|
| 1412 |
//loopze(step/2,3*step/4.,zmin,zmax,500); // |
| 1413 |
|
| 1414 |
if ( debug ) printf("Zdaloop Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step); |
| 1415 |
|
| 1416 |
//chi2,z,Etot,Pskip |
| 1417 |
lpchi2=bestchi2[0]; |
| 1418 |
lpz=bestchi2[1]; |
| 1419 |
lpetot=bestchi2[2]; |
| 1420 |
lppskip=bestchi2[3]; |
| 1421 |
lpisotope=bestchi2[4]; |
| 1422 |
};//endZdaloop |
| 1423 |
|
| 1424 |
|
| 1425 |
|
| 1426 |
|
| 1427 |
|
| 1428 |
|
| 1429 |
|
| 1430 |
|
| 1431 |
|
| 1432 |
|
| 1433 |
|
| 1434 |
|
Parent Directory
| 
Revision Log