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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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