/[PAMELA software]/calo/flight/CaloBragg/src/CaloBragg.cpp
ViewVC logotype

Contents of /calo/flight/CaloBragg/src/CaloBragg.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.13 - (show annotations) (download)
Thu Oct 6 09:27:39 2011 UTC (13 years, 3 months ago) by mocchiut
Branch: MAIN
Changes since 1.12: +348 -88 lines
Andrea's isotope implementation + bugs fixed

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

  ViewVC Help
Powered by ViewVC 1.1.23