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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.11 - (show annotations) (download)
Mon Dec 14 14:35:53 2009 UTC (14 years, 11 months ago) by mocchiut
Branch: MAIN
CVS Tags: HEAD
Changes since 1.10: +6 -0 lines
Error occurred while calculating annotation data.
Do not use plane 18x by default

1 /**
2 * \file CaloFranzini.cpp
3 * \author Emiliano Mocchiutti (2007/07/18)
4 */
5 //
6 // headers
7 //
8 #include <CaloFranzini.h>
9 //--------------------------------------
10 /**
11 * Default constructor
12 */
13 CaloFranzini::CaloFranzini(){
14 Clear();
15 }
16
17 CaloFranzini::CaloFranzini(PamLevel2 *l2p){
18 //
19 lfile = NULL;
20 ffile = NULL;
21 brig = NULL;
22 brigm = NULL;
23 nbin = 0;
24 //
25 L2 = l2p;
26 //
27 if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
28 //
29 // Default variables
30 //
31 debug = false;
32 dolong = true;
33 dofull = false;
34 sntr = 0;
35 OBT = 0;
36 PKT = 0;
37 atime = 0;
38 //
39 usepl18x = false;
40 //
41 crig = false;
42 sel = true;
43 cont = false;
44 N = 0;
45 NC = 43;
46 //
47 mask18b = -1;
48 //
49 Clear();
50 //
51 }
52
53 void CaloFranzini::Clear(){
54 //
55 longtzeta = 0.;
56 fulltzeta = 0.;
57 degfre = 0;
58 fdegfre = 0;
59 memset(estrip, 0, 2*22*96*sizeof(Float_t));
60 memset(qplane, 0, 43*sizeof(Float_t));
61 //
62 numneg = 0;
63 numpos = 0;
64 negfulltzeta = 0.;
65 posfulltzeta = 0.;
66 aveposvar = 0.;
67 avenegvar = 0.;
68 minsvalue = numeric_limits<UInt_t>::max();
69 maxsvalue = numeric_limits<UInt_t>::min();
70 //
71 }
72
73 void CaloFranzini::Print(){
74 //
75 Process();
76 //
77 printf("========================================================================\n");
78 printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
79 printf(" debug [debug flag]:.. %i\n",debug);
80 printf(" long degree of freedom :.. %i\n",this->GetLongDegreeOfFreedom());
81 printf(" longtzeta :.. %f\n",longtzeta);
82 printf(" longtzeta normalized :.. %f\n",this->GetNormLongTZeta());
83 printf(" full degree of freedom :.. %i\n",this->GetFullDegreeOfFreedom());
84 printf(" fulltzeta :.. %f\n",fulltzeta);
85 printf(" fulltzeta normalized :.. %f\n",this->GetNormFullTZeta());
86 printf(" fulltzeta negative contribute :.. %f\n",negfulltzeta);
87 printf(" fulltzeta positive contribute :.. %f\n",posfulltzeta);
88 printf(" fulltzeta number of negatives :.. %i\n",numneg);
89 printf(" fulltzeta number of positives :.. %i\n",numpos);
90 printf(" fulltzeta minimum variance :.. %f\n",minsvalue);
91 printf(" fulltzeta maximum variance :.. %f\n",maxsvalue);
92 printf(" fulltzeta average positive var :.. %f\n",aveposvar);
93 printf(" fulltzeta average negative var :.. %f\n",avenegvar);
94 printf("========================================================================\n");
95 //
96 }
97
98 void CaloFranzini::Delete(){
99 //
100 if ( ffile ) ffile->Close();
101 if ( lfile ) lfile->Close();
102 //
103 Clear();
104 //
105 }
106
107 void CaloFranzini::SetNoWpreSampler(Int_t n){
108 Int_t nc2 = NC/2;
109 if ( NC >= 37 ) nc2 = (NC+1)/2;
110 if ( nc2+n < 23 ){
111 N = n;
112 } else {
113 printf(" ERROR! Calorimeter is made of 22 W planes\n");
114 printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
115 printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
116 NC = 43;
117 N = 0;
118 };
119 }
120
121 void CaloFranzini::SetNoWcalo(Int_t n){
122 if ( N+n < 23 ){
123 NC = n*2;
124 if ( NC >37 ) NC--;
125 } else {
126 printf(" ERROR! Calorimeter is made of 22 W planes\n");
127 printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
128 printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
129 NC = 43;
130 N = 0;
131 };
132 }
133
134 void CaloFranzini::SplitInto(Int_t NoWpreSampler, Int_t NoWcalo){
135 this->SetNoWpreSampler(0);
136 this->SetNoWcalo(0);
137 if ( NoWpreSampler < NoWcalo ){
138 this->SetNoWpreSampler(NoWpreSampler);
139 this->SetNoWcalo(NoWcalo);
140 } else {
141 this->SetNoWcalo(NoWcalo);
142 this->SetNoWpreSampler(NoWpreSampler);
143 };
144 }
145
146
147 void CaloFranzini::Process(){
148 this->Process(0);
149 }
150
151 void CaloFranzini::Process(Int_t itr){
152 //
153 if ( !L2 ){
154 printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
155 printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
156 return;
157 };
158 //
159 if ( !nbin || !brig || (!lfile && !ffile) ){
160 printf(" ERROR: it seems covariance matrix file has not been opened (CaloFranzini::Open()) \n");
161 printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
162 return;
163 };
164 //
165 Bool_t newentry = false;
166 //
167 if ( L2->IsORB() ){
168 if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || itr != sntr ){
169 newentry = true;
170 OBT = L2->GetOrbitalInfo()->OBT;
171 PKT = L2->GetOrbitalInfo()->pkt_num;
172 atime = L2->GetOrbitalInfo()->absTime;
173 sntr = itr;
174 };
175 } else {
176 newentry = true;
177 };
178 //
179 if ( !newentry ) return;
180 //
181 // Some variables
182 //
183 if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
184 //
185 this->Clear();
186 //
187 Float_t rig = L2->GetTrack(itr)->GetTrkTrack()->GetRigidity();
188 if ( crig ) rig = L2->GetCaloLevel2()->qtot/260.;
189 //
190 // Fill the estrip matrix
191 //
192 Int_t nplane = 0;
193 Int_t view = 0;
194 Int_t plane = 0;
195 Int_t strip = 0;
196 Float_t mip = 0.;
197 //
198 //
199 //
200 if ( dolong ){
201 //
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 estrip[view][plane][strip] = mip;
209 //
210 nplane = 1 - view + 2 * (plane - N);
211 if ( nplane > (37-(2*N)) ) nplane--;
212 //
213 // if ( plane == (18+N) ) mip = 0.;
214 if ( nplane > -1 ) qplane[nplane] += mip;
215 //
216 };
217 //
218 if ( cont ){
219 for (Int_t i=0; i<22; i++){
220 if ( i == (18+N) ){
221 mask18b = 1 + 2 * (i - N);
222 break;
223 };
224 };
225 };
226 //
227 if ( sel ){
228 for (Int_t i=0; i<22; i++){
229 if ( i == (18-N) ){
230 mask18b = 1 + 2 * (i - N);
231 break;
232 };
233 };
234 };
235 //
236 if ( mask18b == 37 ) mask18b = -1;
237 //
238 Int_t dgf = 43;
239 //
240 for (Int_t i=0; i < 22; i++){
241 if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){
242 dgf = 2 * i;
243 break;
244 };
245 if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){
246 dgf = 1 + 2 * i;
247 break;
248 };
249 };
250 //
251 if ( dgf < 43 && dgf > 37 ) dgf--;
252 //
253 degfre = TMath::Min(dgf,NC);
254 //
255 // Float_t longzdiag = 0.;
256 // Float_t longzout = 0.;
257 //
258 if ( degfre > 0 ){
259 for (Int_t i = 0; i < degfre; i++){
260 if ( i != mask18b ){
261 for (Int_t j = 0; j < degfre; j++){
262 if ( j != mask18b ){
263 // if ( i == j ){
264 // longzdiag += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
265 // if ( debug ) printf(" %i %i %f %f %f %f %f\n",i,j,qplane[i],this->GetAverageAt(i,rig),this->GetHmatrixAt(i,j,rig),qplane[j],this->GetAverageAt(j,rig));
266 // } else {
267 // longzout += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
268 // if ( debug && i == (j+1) ) printf(" %i %i %f %f %f %f %f\n",i,j,qplane[i],this->GetAverageAt(i,rig),this->GetHmatrixAt(i,j,rig),qplane[j],this->GetAverageAt(j,rig));
269 // };
270 //
271 longtzeta += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
272 //
273 };
274 };
275 };
276 };
277 // if ( debug ) printf(" diagonal total %f out of diagonal total %f \n",longzdiag,longzout);
278 };
279 };
280 if ( dofull ){
281 // printf(" ERROR: NOT IMPLEMENTED YET\n");
282 // printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
283 //
284 // FULL CALORIMETER
285 //
286 CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
287 //
288 Int_t dgf = 0;
289 Int_t cs = 0;
290 Int_t cd = 0;
291 Int_t mstrip = 0;
292 //
293 Int_t maxnpl = 0;
294 Float_t mipv[43][11];
295 memset(mipv,0,43*11*sizeof(Float_t));
296 //
297 for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
298 //
299 mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
300 //
301 if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
302 //
303 // nplane = 1 - view + 2 * plane;
304 // if ( nplane > 37 ) nplane--;
305 nplane = 1 - view + 2 * (plane - N);
306 if ( nplane > (37-(2*N)) ) nplane--;
307 //
308 cs = ct->tibar[plane][view] - 1;
309 //
310 if ( ct->tibar[plane][view] > -1 ){
311 //
312 cd = 95 - cs;
313 //
314 mstrip = cd + strip;
315 //
316 Int_t lstr = this->ConvertStrip(mstrip);
317 //
318 if ( nplane > maxnpl ) maxnpl = nplane;
319 if ( nplane > -1 ) mipv[nplane][lstr] += mip;
320 //
321 };
322 // mipv[nplane][lstr] += mip;
323 //
324 };
325 //
326 Float_t mip1 = 1.;
327 Int_t cs1;
328 Int_t cd1;
329 Float_t mip2 = 1.;
330 Int_t cs2;
331 Int_t cd2;
332 Int_t mi = -1;
333 Int_t mj = -1;
334 Int_t nn1 = 0;
335 Int_t pl1 = 0;
336 Int_t vi1 = 0;
337 Int_t nn2 = 0;
338 Int_t pl2 = 0;
339 Int_t vi2 = 0;
340 Int_t mstrip1min = 0;
341 Int_t mstrip1max = 0;
342 Int_t mstrip2min = 0;
343 Int_t mstrip2max = 0;
344 //
345 Int_t tonpl = TMath::Min(maxnpl,NC);
346 //
347 Int_t rbi = 0;
348 for (Int_t i = 0; i<nbin-1; i++){
349 if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
350 rbi = i;
351 break;
352 };
353 };
354 //
355 //
356 Int_t therigb = rbi;
357 //
358 if ( rig < brigm->At(2) ){
359 // therigb = 0;
360 therigb = 2;
361 };
362 // if ( rig >= brigm->At(nbin-5) ){
363 if ( rig >= brigm->At(nbin-2) ){
364 therigb = nbin - 3;
365 // therigb = nbin - 5;
366 };
367 Int_t mtherig = therigb;
368 if ( mtherig >= 13 ) mtherig = 12;
369 //
370 if ( debug ) printf(" rig %f brigm0 %f brigm n2 %f nbin %i \n",rig,brigm->At(0),brigm->At(nbin-2),nbin);
371 //
372 if ( cont ){
373 for (Int_t i=0; i<22; i++){
374 if ( i == (18+N) ){
375 mask18b = 1 + 2 * (i - N);
376 break;
377 };
378 };
379 };
380 //
381 if ( sel ){
382 for (Int_t i=0; i<22; i++){
383 if ( i == (18-N) ){
384 mask18b = 1 + 2 * (i - N);
385 break;
386 };
387 };
388 };
389 //
390 if ( mask18b == 37 ) mask18b = -1;
391 //
392 if ( debug ) printf("Start main loop \n");
393 //
394 for (Int_t nplane1=0; nplane1<tonpl; nplane1++){
395 //
396 if ( nplane1 != mask18b ){
397 //
398 if ( nplane1 >= 37 ) nn1 = nplane1 + 1;
399 vi1 = 1;
400 if ( nn1%2 ) vi1 = 0;
401 pl1 = (nn1 - 1 + vi1)/2;
402 //
403 cs1 = ct->tibar[pl1][vi1] - 1;
404 //
405 if ( ct->tibar[pl1][vi1] > -1 ){
406 //
407 dgf++;
408 //
409 cd1 = 95 - cs1;
410 //
411 Int_t at1 = TMath::Max(0,(cd1+0));
412 Int_t at2 = TMath::Min(190,(cd1+95));
413 mstrip1min = this->ConvertStrip(at1);
414 mstrip1max = this->ConvertStrip(at2) + 1;
415 //
416 for (Int_t mstrip1=mstrip1min; mstrip1<mstrip1max; mstrip1++){
417 //
418 mj = -1;
419 //
420 // if ( debug ) printf(" get mip1 %i %i %f %i \n",nplane1,mstrip1,rig,therigb);
421 mip1 = mipv[nplane1][mstrip1] - this->GetFullAverageAt(nplane1,mstrip1,rig,therigb);
422 if ( mip1 < 0 ){
423 numneg++;
424 avenegvar += mip1;
425 };
426 if ( mip1 >= 0 ){
427 numpos++;
428 aveposvar += mip1;
429 };
430 if ( mip1 < minsvalue ) minsvalue = mip1;
431 if ( mip1 >= maxsvalue ) maxsvalue = mip1;
432 //
433 mi = (nplane1 * 11) + mstrip1;
434 //
435 if ( mip1 != 0. ){
436 //
437 for (Int_t nplane2=0; nplane2<tonpl; nplane2++){
438 //
439 if ( nplane2 != mask18b ){
440 //
441 if ( nplane2 >= 37 ) nn2 = nplane2 + 1;
442 vi2 = 1;
443 if ( nn2%2 ) vi2 = 0;
444 pl1 = (nn2 - 1 + vi2)/2;
445 //
446 cs2 = ct->tibar[pl2][vi2] - 1;
447 //
448 if ( ct->tibar[pl2][vi2] > -1 ){
449 //
450 cd2 = 95 - cs2;
451 //
452 Int_t t1 = TMath::Max(0,(cd2+0));
453 Int_t t2 = TMath::Min(190,(cd2+95));
454 mstrip2min = this->ConvertStrip(t1);
455 mstrip2max = this->ConvertStrip(t2) + 1;
456 //
457 for (Int_t mstrip2=mstrip2min; mstrip2<mstrip2max; mstrip2++){
458 //
459 mip2 = mipv[nplane2][mstrip2] - this->GetFullAverageAt(nplane2,mstrip2,rig,therigb);
460 //
461 if ( mip2 != 0. ){
462 //
463 // Int_t sh = -5 + nplane1;
464 // if ( sh > 5 ) sh -= 11*nplane1;
465 // //
466 // mj = (nplane2 * 11) + mstrip2 + sh;
467 // //
468 // if ( mj < 0 ) mj += 473;
469 // if ( mj >= 473 ) mj -= 473;
470 // //
471 mj = (nplane2 * 11) + mstrip2;
472 //
473 Float_t svalue = mip1 * this->GetFullHmatrixAt(mi,mj,rig,mtherig,therigb) * mip2;
474 // fulltzeta += mip1 * this->GetFullHmatrixAt(mi,mj,rig,therigb) * mip2;
475 fulltzeta += svalue;
476 if ( svalue < 0. ) negfulltzeta += svalue;
477 if ( svalue > 0. ) posfulltzeta += svalue;
478 // (*fmatrix[rbi])[mi][mj] += (mip1 * mip2);
479 // if ( mstrip1 == mstrip2 ) printf("\n\n=> nplane1 %i nplane2 %i mstrip1 %i mstrip2 %i \n => mip1 %f H %f mip2 %f \n => mipv1 %f ave1 %f mipv2 %f ave2 %f \n => rig %f rigbin %i mtherigb %i\n",nplane1,nplane2,mstrip1,mstrip2,mip1,this->GetFullHmatrixAt(mi,mj,rig,mtherig),mip2,mipv[nplane1][mstrip1],this->GetFullAverageAt(nplane1,mstrip1,rig,therigb),mipv[nplane2][mstrip2],this->GetFullAverageAt(nplane2,mstrip2,rig,therigb),rig,therigb,mtherig);
480 //
481 };
482 };
483 };
484 };
485 };
486 };
487 };
488 };
489 };
490 };
491 //
492 fdegfre = dgf*11;
493 if ( numpos ) aveposvar /= numpos;
494 if ( numneg ) avenegvar /= numneg;
495 //
496 };
497 //
498 // if ( debug ) this->Print();
499 if ( debug ) printf(" exit \n");
500 }
501
502
503 Float_t CaloFranzini::GetNormLongTZeta(){
504 Process();
505 Float_t normz = 0.;
506 if ( degfre != 0 ) normz = longtzeta/(Float_t)degfre;
507 return normz;
508 }
509
510 Float_t CaloFranzini::GetNormFullTZeta(){
511 Process();
512 Float_t normz = 0.;
513 if ( fdegfre != 0 ) normz = fulltzeta/(Float_t)fdegfre;
514 return normz;
515 }
516
517 Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){
518 //
519 lfile = new TFile(matrixfile.Data(),"READ");
520 //
521 if ( !lfile || lfile->IsZombie() ){
522 if ( dolong ){
523 lfile = new TFile(matrixfile.Data(),"RECREATE");
524 printf(" Create file %s \n",lfile->GetName());
525 };
526 if ( dofull ){
527 ffile = new TFile(matrixfile.Data(),"RECREATE");
528 printf(" Create file %s \n",ffile->GetName());
529 };
530 } else {
531 lfile->Close();
532 printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
533 return(false);
534 };
535 //
536 return(true);
537 //
538 }
539
540 Bool_t CaloFranzini::UpdateMatrixFile(TString matrixfile){
541 //
542 if ( dolong ){
543 lfile = new TFile(matrixfile.Data(),"UPDATE");
544 //
545 if ( !lfile || lfile->IsZombie() ){
546 printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
547 return(false);
548 };
549 };
550 //
551 if ( dofull ){
552 ffile = new TFile(matrixfile.Data(),"UPDATE");
553 //
554 if ( !ffile || ffile->IsZombie() ){
555 printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
556 return(false);
557 };
558 };
559 //
560 return(true);
561 //
562 }
563
564 void CaloFranzini::WriteNumBin(Int_t numbin){
565 if ( dolong ){
566 lfile->cd();
567 TArrayI nbi(1, &numbin);
568 lfile->WriteObject(&nbi, "nbinenergy");
569 };
570 if ( dofull ){
571 ffile->cd();
572 TArrayI nbi(1, &numbin);
573 ffile->WriteObject(&nbi, "nbinenergy");
574 };
575 }
576
577 void CaloFranzini::WriteRigBin(TArrayF *rigbin){
578 if ( dolong ){
579 lfile->cd();
580 // rigbin->Write("binrig");
581 lfile->WriteObject(&(*rigbin), "binrig");
582 };
583 if ( dofull ){
584 ffile->cd();
585 // rigbin->Write("binrig");
586 ffile->WriteObject(&(*rigbin), "binrig");
587 };
588 }
589
590 void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){
591 lfile->cd();
592 TString name = Form("qplmeann%i",bin);
593 lfile->WriteObject(&(*qpl),name.Data());
594 }
595
596 void CaloFranzini::WriteFullMean(TMatrixD *qpl, Int_t bin){
597 ffile->cd();
598 TString name = Form("fqplmeann%i",bin);
599 ffile->WriteObject(&(*qpl),name.Data());
600 ffile->Purge();
601 }
602
603 void CaloFranzini::WriteFullNMean(TMatrixD *qpl, Int_t bin){
604 ffile->cd();
605 TString name = Form("fnqplmeann%i",bin);
606 ffile->WriteObject(&(*qpl),name.Data());
607 ffile->Purge();
608 }
609
610 void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){
611 lfile->cd();
612 TString name = Form("matrixn%i",bin);
613 // mat.Write(name.Data());
614 lfile->WriteObject(&mat,name.Data());
615 }
616
617 void CaloFranzini::WriteInvertedFullMatrix(TMatrixD mat, Int_t bin){
618 ffile->cd();
619 TString name = Form("fmatrixn%i",bin);
620 // mat.Write(name.Data());
621 ffile->WriteObject(&mat,name.Data());
622 }
623
624 void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
625 lfile->cd();
626 TString name = Form("origmatrixn%i",bin);
627 // mat.Write(name.Data());
628 lfile->WriteObject(&(*mat),name.Data());
629 }
630
631 void CaloFranzini::WriteFullMatrix(TMatrixD *mat, Int_t bin){
632 ffile->cd();
633 TString name = Form("origfmatrixn%i",bin);
634 // mat.Write(name.Data());
635 ffile->WriteObject(&(*mat),name.Data());
636 ffile->Purge();
637 }
638
639 void CaloFranzini::WriteFullNMatrix(TMatrixF *mat, Int_t bin){
640 ffile->cd();
641 TString name = Form("origfnmatrixn%i",bin);
642 // mat.Write(name.Data());
643 ffile->WriteObject(&(*mat),name.Data());
644 ffile->Purge();
645 }
646
647 void CaloFranzini::CloseMatrixFile(){
648 if ( dolong ){
649 lfile->cd();
650 lfile->Close();
651 };
652 if ( dofull ){
653 ffile->cd();
654 ffile->Close();
655 };
656 }
657
658 Bool_t CaloFranzini::Open(TString matrixfile){
659 return(this->Open(matrixfile,""));
660 }
661
662 Bool_t CaloFranzini::Open(TString longmatrixfile, TString fullmatrixfile){
663 //
664 // find matrix file
665 //
666 if ( !strcmp(longmatrixfile.Data(),"") ){
667 if (dolong) longmatrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_longel.root";
668 };
669 if ( !strcmp(fullmatrixfile.Data(),"") ){
670 if (dofull) fullmatrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_fullel.root";
671 };
672 //
673 if ( dolong ){
674 //
675 lfile = new TFile(longmatrixfile.Data(),"READ");
676 //
677 if ( !lfile || lfile->IsZombie() ){
678 printf(" ERROR: cannot open file %s \n",longmatrixfile.Data());
679 return(false);
680 };
681 //
682 if ( !this->LoadBin() ){
683 printf(" %s \n",longmatrixfile.Data());
684 return(false);
685 };
686 //
687 if ( !this->LoadLong() ){
688 printf(" %s \n",longmatrixfile.Data());
689 return(false);
690 };
691 //
692 if ( !this->LoadMatrices() ){
693 printf(" %s \n",longmatrixfile.Data());
694 return(false);
695 };
696 };
697 //
698 if ( dofull ){
699 //
700 ffile = new TFile(fullmatrixfile.Data(),"READ");
701 //
702 if ( !ffile || ffile->IsZombie() ){
703 printf(" ERROR: cannot open file %s \n",fullmatrixfile.Data());
704 return(false);
705 };
706 //
707 if ( !dolong ){
708 //
709 if ( !this->LoadBin(true) ){
710 printf(" %s \n",fullmatrixfile.Data());
711 return(false);
712 };
713 //
714 };
715 if ( !this->LoadFull() ){
716 printf(" %s \n",fullmatrixfile.Data());
717 return(false);
718 };
719 //
720 if ( !this->LoadFullMatrices() ){
721 printf(" %s \n",fullmatrixfile.Data());
722 return(false);
723 };
724 };
725 //
726 //
727 return(true);
728 //
729 }
730
731
732 Bool_t CaloFranzini::LoadBin(){
733 return(this->LoadBin(false));
734 }
735
736 Bool_t CaloFranzini::LoadBin(Bool_t full){
737 //
738 if ( !full ) {
739 //
740 TArrayI *numbin = (TArrayI*)lfile->Get("nbinenergy");
741 if ( !numbin ){
742 printf(" ERROR: cannot read number of bins from file ");
743 return(false);
744 };
745 nbin = (Int_t)numbin->At(0);
746 if ( nbin <= 0 ){
747 printf(" ERROR: cannot work with 0 energy bins from file ");
748 return(false);
749 };
750 //
751 brig = (TArrayF*)lfile->Get("binrig");
752 if ( !brig ){
753 printf(" ERROR: cannot read rigidity binning from file ");
754 return(false);
755 };
756 //
757 brigm=(TArrayF*)lfile->Get("binrigmean");
758 if ( !brigm ){
759 printf(" ERROR: cannot read mean rigidity binning from file ");
760 return(false);
761 };
762 //
763 } else {
764 //
765 TArrayI *numbin = (TArrayI*)ffile->Get("nbinenergy");
766 if ( !numbin ){
767 printf(" ERROR: cannot read number of bins from file ");
768 return(false);
769 };
770 nbin = (Int_t)numbin->At(0);
771 if ( nbin <= 0 ){
772 printf(" ERROR: cannot work with 0 energy bins from file ");
773 return(false);
774 };
775 //
776 brig = (TArrayF*)ffile->Get("binrig");
777 if ( !brig ){
778 printf(" ERROR: cannot read rigidity binning from file ");
779 return(false);
780 };
781 //
782 brigm=(TArrayF*)ffile->Get("binrigmean");
783 if ( !brigm ){
784 printf(" ERROR: cannot read mean rigidity binning from file ");
785 return(false);
786 };
787 };
788 //
789 return(true);
790 }
791
792 Bool_t CaloFranzini::LoadLong(){
793 //
794 for (Int_t i=0;i<17;i++){
795 TString name = Form("qplmeann%i",i);
796 qplmean[i] = (TArrayF*)lfile->Get(name.Data());
797 if ( !qplmean[i] ){
798 printf(" ERROR: cannot read average from file ");
799 return(false);
800 };
801 };
802 //
803 return(true);
804 }
805
806 Bool_t CaloFranzini::LoadFull(){
807 //
808 for (Int_t i=0;i<17;i++){
809 TString name = Form("fqplmeann%i",i);
810 fqplmean[i] = (TMatrixD*)ffile->Get(name.Data());
811 if ( !fqplmean[i] ){
812 printf(" ERROR: cannot read average from file ");
813 return(false);
814 };
815 };
816 //
817 return(true);
818 }
819
820 Bool_t CaloFranzini::LoadMatrices(){
821 //
822 for (Int_t i=0;i<17;i++){
823 TString name1 = Form("matrixn%i",i);
824 hmat[i] = (TMatrixD*)lfile->Get(name1.Data());
825 };
826 //
827 return(true);
828 }
829
830 Bool_t CaloFranzini::LoadFullMatrices(){
831 //
832 for (Int_t i=0;i<17;i++){
833 TString name1 = Form("fmatrixn%i",i);
834 hfmat[i] = (TMatrixD*)ffile->Get(name1.Data());
835 };
836 //
837 return(true);
838 }
839
840 TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){
841 //
842 Int_t mv = 0;
843 for (Int_t i = 0; i<nbin-1; i++){
844 if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
845 mv = i;
846 break;
847 };
848 };
849 if ( rig < brig->At(0) ){
850 printf(" WARNING: Event with rigidity lower than the first covariance matrix bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
851 mv = 0;
852 };
853 if ( rig >= brig->At(nbin-1) ){
854 printf(" WARNING: Event with rigidity higher than the last covariance matrix bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
855 mv = nbin-2;
856 };
857 //
858 return(hmat[mv]);
859 //
860 }
861
862 TMatrixD *CaloFranzini::LoadFullAverage(Int_t rigbin){
863 //
864 TString name = Form("fqplmeann%i",rigbin);
865 TMatrixD *fmean=(TMatrixD*)ffile->Get(name.Data());
866 //
867 return(fmean);
868 //
869 }
870
871 TMatrixF *CaloFranzini::LoadFullMatrix(Int_t rigbin){
872 //
873 TString name = Form("origfmatrixn%i",rigbin);
874 TMatrixF *fmatri=(TMatrixF*)ffile->Get(name.Data());
875 //
876 return(fmatri);
877 //
878 }
879
880 void CaloFranzini::LoadFullMatrix(Int_t rigbin, TMatrixF *fmatri){
881 //
882 TString name = Form("origfmatrixn%i",rigbin);
883 fmatri=(TMatrixF*)ffile->Get(name.Data());
884 //
885 }
886
887 void CaloFranzini::UnLoadFullAverage(Int_t rigbin){
888 //
889 TString name = Form("fqplmeann%i",rigbin);
890 ffile->Delete(name.Data());
891 //
892 }
893
894 void CaloFranzini::UnLoadFullMatrix(Int_t rigbin){
895 //
896 TString name = Form("origfmatrixn%i",rigbin);
897 ffile->Delete(name.Data());
898 //
899 }
900
901 TMatrixF *CaloFranzini::LoadFullNMatrix(Int_t rigbin){
902 //
903 TString name = Form("origfnmatrixn%i",rigbin);
904 TMatrixF *fnmatri=(TMatrixF*)ffile->Get(name.Data());
905 //
906 return(fnmatri);
907 //
908 }
909
910 void CaloFranzini::UnLoadFullNMatrix(Int_t rigbin){
911 //
912 TString name = Form("origfnmatrixn%i",rigbin);
913 ffile->Delete(name.Data());
914 //
915 }
916
917 TMatrixD *CaloFranzini::LoadFullNAverage(Int_t rigbin){
918 //
919 TString name = Form("fnqplmeann%i",rigbin);
920 TMatrixD *fnmean=(TMatrixD*)ffile->Get(name.Data());
921 //
922 return(fnmean);
923 //
924 }
925
926 void CaloFranzini::UnLoadFullNAverage(Int_t rigbin){
927 //
928 TString name = Form("fnqplmeann%i",rigbin);
929 ffile->Delete(name.Data());
930 //
931 }
932
933
934 TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){
935 //
936 Int_t mv=0;
937 for (Int_t i = 0; i<nbin-1; i++){
938 if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
939 mv = i;
940 break;
941 };
942 };
943 if ( rig < brig->At(0) ){
944 printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
945 mv = 0;
946 };
947 if ( rig >= brig->At(nbin-1) ){
948 printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
949 mv=nbin-2;
950 };
951 //
952 return(qplmean[mv]);
953 //
954 }
955
956 Float_t CaloFranzini::GetAverageAt(Int_t plane, Float_t rig){
957 //
958 Int_t therigb = 0;
959 for (Int_t i = 0; i<nbin-2; i++){
960 if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
961 therigb = i;
962 break;
963 };
964 };
965 //
966 Float_t minrig;
967 Float_t maxrig;
968 //
969 //
970 if ( rig < brigm->At(0) ){
971 if ( rig < brig->At(0) ){
972 // printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
973 };
974 therigb = 0;
975 };
976 if ( rig >= brigm->At(nbin-2) ){
977 if ( rig >= brig->At(nbin-2) ) {
978 // printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
979 };
980 therigb = nbin - 3;
981 };
982 //
983 minrig = brigm->At(therigb);
984 maxrig = brigm->At(therigb+1);
985 //
986 Float_t minene = (*qplmean[therigb])[plane];
987 Float_t maxene = (*qplmean[therigb+1])[plane];
988 //
989 if ( maxrig == minrig ){
990 printf("Unrecoverable ERROR! Matrix will be screwed up... \n");
991 return(0.);
992 };
993 Float_t b = log(maxene/minene)/(maxrig-minrig);
994 Float_t a = minene/exp(minrig*b);
995 Float_t ave = a*exp(b*rig);
996 if ( a == 0. || minene == 0. || ave != ave ){
997 // if ( a == 0. || minene == 0. ){
998 Float_t m = (maxene-minene)/(maxrig-minrig);
999 Float_t q = minene - m * minrig;
1000 ave = rig * m + q;
1001 };
1002 //
1003 return(ave);
1004 //
1005 }
1006
1007 Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig){
1008 //
1009 Int_t therigb = 0;
1010 for (Int_t i = 0; i<nbin-2; i++){
1011 if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
1012 therigb = i;
1013 break;
1014 };
1015 };
1016 //
1017 //
1018 if ( rig < brigm->At(0) ){
1019 if ( rig < brig->At(0) ){
1020 // printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
1021 };
1022 therigb = 0;
1023 };
1024 if ( rig >= brigm->At(nbin-2) ){
1025 if ( rig >= brig->At(nbin-2) ) {
1026 // printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
1027 };
1028 therigb = nbin - 3;
1029 };
1030 //
1031 return(this->GetFullAverageAt(plane,strip,rig,therigb));
1032 //
1033 }
1034
1035 Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig, Int_t therigb){
1036 //
1037 Float_t minrig;
1038 Float_t maxrig;
1039 //
1040 minrig = brigm->At(therigb);
1041 maxrig = brigm->At(therigb+1);
1042 //
1043 Float_t minene = (*fqplmean[therigb])[plane][strip];
1044 Float_t maxene = (*fqplmean[therigb+1])[plane][strip];
1045 //
1046 if ( maxrig == minrig ){
1047 printf("Unrecoverable ERROR! Matrix will be screwed up... \n");
1048 return(0.);
1049 };
1050 Float_t b = log(maxene/minene)/(maxrig-minrig);
1051 Float_t a = minene/exp(minrig*b);
1052 Float_t ave = a*exp(b*rig);
1053 if ( a == 0. || minene == 0. || ave != ave ){
1054 Float_t m = (maxene-minene)/(maxrig-minrig);
1055 Float_t q = minene - m * minrig;
1056 ave = rig * m + q;
1057 };
1058 //
1059 // ave += (44.-plane)*strip;
1060 //if ( a == 0. ) ave = 0.;
1061 if ( !(ave == ave) ) printf("a %f b %f ave %f maxene %f minene %f maxrig %f minrig %f \n",a,b,ave,maxene,minene,maxrig,minrig);
1062 //
1063 return(ave);
1064 //
1065 }
1066
1067 Float_t CaloFranzini::GetHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
1068 Int_t therigb = 0;
1069 for (Int_t i = 0; i<nbin-2; i++){
1070 if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
1071 therigb = i;
1072 break;
1073 };
1074 };
1075 //
1076 Float_t minrig;
1077 Float_t maxrig;
1078 //
1079 if ( rig < brigm->At(0) ){
1080 if ( rig < brig->At(0) ){
1081 // printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
1082 };
1083 therigb = 0;
1084 };
1085 // if ( rig >= brigm->At(nbin-4) ){ // -2
1086 if ( rig >= brigm->At(nbin-2) ){
1087 if ( rig >= brig->At(nbin-2) ) {
1088 // printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
1089 };
1090 // therigb = nbin - 5;// -3
1091 therigb = nbin - 3;
1092 };
1093 //
1094 if ( therigb < 2 ) therigb = 2;
1095 minrig = brigm->At(therigb);
1096 maxrig = brigm->At(therigb+1);
1097 // printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);
1098 //
1099 Float_t minene = (*hmat[therigb])[iindex][jindex];
1100 Float_t maxene = (*hmat[therigb+1])[iindex][jindex];
1101 // printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,a,b,rig,ave);
1102 //
1103 // Float_t a = 0.;
1104 // Float_t b = 0.;
1105 // Float_t ave = 0.;
1106 // if ( minene == 0. ){
1107 //
1108 // } else {
1109 // b = log(maxene/minene)/(maxrig-minrig);
1110 // a = minene/exp(minrig*b);
1111 // ave = a*exp(b*rig);
1112 // };
1113 //
1114 Float_t m = (maxene-minene)/(maxrig-minrig);
1115 Float_t q = minene - m * minrig;
1116 Float_t ave = rig * m + q;
1117 if ( debug ) printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,m,q,rig,ave);
1118 //
1119 //
1120 return(ave);
1121 //
1122 }
1123
1124 Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
1125 Int_t therigb = 0;
1126 for (Int_t i = 0; i<nbin-2; i++){
1127 if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
1128 therigb = i;
1129 break;
1130 };
1131 };
1132 //
1133 if ( rig < brigm->At(0) ){
1134 if ( rig < brig->At(0) ){
1135 // printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
1136 };
1137 therigb = 0;
1138 };
1139 // if ( rig >= brigm->At(nbin-4) ){ // -2
1140 if ( rig >= brigm->At(nbin-2) ){
1141 //if ( rig >= brigm->At(nbin-5) ){
1142 if ( rig >= brig->At(nbin-2) ) {
1143 // printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
1144 };
1145 // therigb = nbin - 5;// -3
1146 // therigb = nbin - 5;// -3
1147 therigb = nbin - 3;
1148 };
1149 //
1150 if ( therigb < 2 ) therigb = 2;
1151 if ( therigb > 13 ) therigb = 13;
1152 //
1153 return(this->GetFullHmatrixAt(iindex,jindex,rig,therigb));
1154 //
1155 }
1156
1157
1158 Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig, Int_t therigb){
1159 //
1160 return(this->GetFullHmatrixAt(iindex,jindex,rig,therigb,0));
1161 //
1162 }
1163
1164 Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig, Int_t therigb, Int_t mtherig){
1165 //
1166 // Int_t lofit = 12;
1167 Int_t lofit = 3;
1168 //
1169 Float_t lowrig = 0.;
1170 Float_t minrig;
1171 Float_t maxrig;
1172 //
1173 if ( mtherig > lofit ) lowrig = brigm->At(therigb-1);
1174 minrig = brigm->At(therigb);
1175 maxrig = brigm->At(therigb+1);
1176 //if ( therigb > 10 ) printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);
1177 //
1178 Float_t lowene = 0.;
1179 if ( mtherig > lofit ) lowene = (*hfmat[therigb-1])[iindex][jindex];
1180 Float_t minene = (*hfmat[therigb])[iindex][jindex];
1181 Float_t maxene = (*hfmat[therigb+1])[iindex][jindex];
1182 // printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,a,b,rig,ave);
1183 //
1184 Float_t ave = 0.;
1185 //
1186 // Float_t a = 0.;
1187 // Float_t b = 0.;
1188 // Float_t ave = 0.;
1189 // if ( minene == 0. ){
1190 //
1191 // } else {
1192 // b = log(maxene/minene)/(maxrig-minrig);
1193 // a = minene/exp(minrig*b);
1194 // ave = a*exp(b*rig);
1195 // };
1196 //
1197 if ( mtherig > lofit ){
1198 //
1199 // FIT
1200 //
1201 // Float_t x[3]={lowrig,minrig,maxrig};
1202 // Float_t y[3]={lowene,minene,maxene};
1203 // //
1204 // TGraph *tg= new TGraph(3,x,y);
1205 // //
1206 // // gStyle->SetLabelSize(0.04);
1207 // // gStyle->SetNdivisions(510,"XY");
1208 // // TCanvas *c = new TCanvas();
1209 // // c->Draw();
1210 // // c->cd();
1211 // // tg->Draw("AP");
1212 // //
1213 // TF1 *fun = new TF1("fun","pol2");
1214 // tg->Fit("fun","QNC");
1215 // //
1216 // ave = fun->GetParameter(0) + rig * fun->GetParameter(1) + rig * rig * fun->GetParameter(2);
1217 // //
1218 // // printf(" therigb %i ave %f rig %f lowrig %f minrig %f maxrig %f lowene %f minene %f maxene %f \n",therigb,ave,rig,lowrig,minrig,maxrig,lowene,minene,maxene);
1219 // //
1220 // //
1221 // // c->Modified();
1222 // // c->Update();
1223 // // gSystem->ProcessEvents();
1224 // // gSystem->Sleep(6000);
1225 // // c->Close();
1226 // // gStyle->SetLabelSize(0);
1227 // // gStyle->SetNdivisions(1,"XY");
1228 // //
1229 // tg->Delete();
1230 // fun->Delete();
1231 Float_t mrigl = (lowrig+minrig)/2.;
1232 Float_t mrigu = (minrig+maxrig)/2.;
1233 Float_t menel = (lowene+minene)/2.;
1234 Float_t meneu = (minene+maxene)/2.;
1235 Float_t m = (meneu-menel)/(mrigu-mrigl);
1236 Float_t q = menel - m * mrigl;
1237 ave = rig * m + q;
1238 //
1239 } else {
1240 Float_t m = (maxene-minene)/(maxrig-minrig);
1241 Float_t q = minene - m * minrig;
1242 ave = rig * m + q;
1243 if ( debug ) printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,m,q,rig,ave);
1244 };
1245 //
1246 //
1247 return(ave);
1248 //
1249 }
1250
1251 void CaloFranzini::DrawLongAverage(Float_t rig){
1252 //
1253 TArrayF *ll = this->LoadLongAverage(rig);
1254 //
1255 gStyle->SetLabelSize(0.04);
1256 gStyle->SetNdivisions(510,"XY");
1257 //
1258 TString hid = Form("cfralongvyvx");
1259 TCanvas *tcf = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
1260 if ( tcf ){
1261 tcf->cd();
1262 } else {
1263 tcf = new TCanvas(hid,hid);
1264 };
1265 //
1266 TString thid = Form("hfralongvyvx");
1267 TH1F *thf = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
1268 if ( thf ) thf->Delete();
1269 thf = new TH1F(thid,thid,44,-0.5,43.5);
1270 tcf->cd();
1271 Float_t qpl[43];
1272 for (Int_t st=0;st<43;st++){
1273 qpl[st] = ll->At(st);
1274 printf("st %i qpl %f\n",st,qpl[st]);
1275 };
1276 for (Int_t st=0;st<44;st++){
1277 if ( st == 37 ){
1278 thf->Fill(st,0.);
1279 } else {
1280 Int_t ss = st;
1281 if ( st > 37 ) ss--;
1282 thf->Fill(st,qpl[ss]);
1283 };
1284 };
1285 thf->Draw();
1286 tcf->Modified();
1287 tcf->Update();
1288 //
1289 gStyle->SetLabelSize(0);
1290 gStyle->SetNdivisions(1,"XY");
1291 //
1292 };
1293
1294 void CaloFranzini::DrawLongAverage(Int_t bin){
1295 //
1296 TArrayF *ll = this->LoadLongAverage(brigm->At(bin));
1297 //
1298 gStyle->SetLabelSize(0.04);
1299 gStyle->SetNdivisions(510,"XY");
1300 //
1301 TString hid = Form("cfralongvyvx");
1302 TCanvas *tcf = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
1303 if ( tcf ){
1304 tcf->cd();
1305 } else {
1306 tcf = new TCanvas(hid,hid);
1307 };
1308 //
1309 TString thid = Form("hfralongvyvx");
1310 TH1F *thf = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
1311 if ( thf ) thf->Delete();
1312 thf = new TH1F(thid,thid,44,-0.5,43.5);
1313 tcf->cd();
1314 Float_t qpl[43];
1315 for (Int_t st=0;st<43;st++){
1316 qpl[st] = ll->At(st);
1317 printf("st %i qpl %f\n",st,qpl[st]);
1318 };
1319 for (Int_t st=0;st<44;st++){
1320 if ( st == 37 ){
1321 thf->Fill(st,0.);
1322 } else {
1323 Int_t ss = st;
1324 if ( st > 37 ) ss--;
1325 thf->Fill(st,qpl[ss]);
1326 };
1327 };
1328 thf->Draw();
1329 tcf->Modified();
1330 tcf->Update();
1331 //
1332 gStyle->SetLabelSize(0);
1333 gStyle->SetNdivisions(1,"XY");
1334 //
1335 };
1336
1337 Int_t CaloFranzini::ConvertStrip(Int_t mstrip){
1338 //
1339 Int_t lastrip = 0;
1340 // //
1341 // if ( mstrip < 50 ) lastrip = 0;
1342 // //
1343 // if ( mstrip >= 50 && mstrip < 64 ) lastrip = 1;
1344 // //
1345 // if ( mstrip >= 64 && mstrip < 70 ) lastrip = 2;
1346 // //
1347 // if ( mstrip >= 70 && mstrip < 75 ) lastrip = 3;
1348 // //
1349 // if ( mstrip >= 75 && mstrip < 84 ){
1350 // lastrip = (int)trunc((mstrip - 75)/3.) + 4;
1351 // };
1352 // //
1353 // if ( mstrip >= 84 && mstrip < 90 ){
1354 // lastrip = (int)trunc((mstrip - 84)/2.) + 7;
1355 // };
1356 // //
1357 // if ( mstrip >= 90 && mstrip < 101 ){
1358 // lastrip = mstrip - 90 + 10;
1359 // };
1360 // //
1361 // if ( mstrip >= 101 && mstrip < 107 ){
1362 // lastrip = (int)trunc((mstrip - 101)/2.) + 21;
1363 // };
1364 // //
1365 // //
1366 // if ( mstrip >= 107 && mstrip < 116 ){
1367 // lastrip = (int)trunc((mstrip - 107)/3.) + 24;
1368 // };
1369 // //
1370 // if ( mstrip >= 116 && mstrip < 121 ) lastrip = 27;
1371 // //
1372 // if ( mstrip >= 121 && mstrip < 127 ) lastrip = 28;
1373 // //
1374 // if ( mstrip >= 127 && mstrip < 141 ) lastrip = 29;
1375 // //
1376 // if ( mstrip >= 141 ) lastrip = 30;
1377 // //
1378 //
1379 if ( mstrip < 83 ) lastrip = 0;
1380 //
1381 if ( mstrip >= 83 && mstrip < 90 ) lastrip = 1;
1382 //
1383 if ( mstrip >= 90 && mstrip < 93 ) lastrip = 2;
1384 //
1385 if ( mstrip >= 93 && mstrip < 98 ){
1386 lastrip = mstrip - 93 + 3;
1387 };
1388 //
1389 if ( mstrip >= 98 && mstrip < 101 ) lastrip = 8;
1390 //
1391 if ( mstrip >= 101 && mstrip < 107 ) lastrip = 9;
1392 //
1393 if ( mstrip >= 107 ) lastrip = 10;
1394 //
1395 return(lastrip);
1396 }

  ViewVC Help
Powered by ViewVC 1.1.23