/[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.11 by mocchiut, Mon Dec 14 14:35:53 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      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();    Clear();
50    //    //
# Line 37  CaloFranzini::CaloFranzini(PamLevel2 *l2 Line 52  CaloFranzini::CaloFranzini(PamLevel2 *l2
52    
53  void CaloFranzini::Clear(){  void CaloFranzini::Clear(){
54    //    //
55    OBT = 0;    longtzeta = 0.;
56    PKT = 0;    fulltzeta = 0.;
   atime = 0;  
   longtzeta = -100.;  
   fulltzeta = -100.;  
57    degfre = 0;    degfre = 0;
58      fdegfre = 0;
59    memset(estrip, 0, 2*22*96*sizeof(Float_t));    memset(estrip, 0, 2*22*96*sizeof(Float_t));
60    memset(qplane, 0, 2*22*sizeof(Float_t));    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    
# Line 55  void CaloFranzini::Print(){ Line 77  void CaloFranzini::Print(){
77    printf("========================================================================\n");    printf("========================================================================\n");
78    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);    printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
79    printf(" debug                [debug flag]:.. %i\n",debug);    printf(" debug                [debug flag]:.. %i\n",debug);
80      printf(" long  degree of freedom          :.. %i\n",this->GetLongDegreeOfFreedom());  
81    printf(" longtzeta                        :.. %f\n",longtzeta);      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);      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");    printf("========================================================================\n");
95    //    //
96  }  }
97    
98  void CaloFranzini::Delete(){  void CaloFranzini::Delete(){
99    //    //
100    if ( file ) file->Close();    if ( ffile ) ffile->Close();
101      if ( lfile ) lfile->Close();
102    //    //
103    Clear();    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(){  void CaloFranzini::Process(){
148    this->Process(0);    this->Process(0);
# Line 82  void CaloFranzini::Process(Int_t itr){ Line 156  void CaloFranzini::Process(Int_t itr){
156      return;      return;
157    };    };
158    //    //
159    if ( !nbin || !file || !brig ){    if ( !nbin || !brig || (!lfile && !ffile) ){
160      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");
161      printf(" ERROR: CaloFranzini variables _NOT_ filled \n");      printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
162      return;      return;
# Line 111  void CaloFranzini::Process(Int_t itr){ Line 185  void CaloFranzini::Process(Int_t itr){
185    this->Clear();    this->Clear();
186    //    //
187    Float_t rig = L2->GetTrack(itr)->GetTrkTrack()->GetRigidity();    Float_t rig = L2->GetTrack(itr)->GetTrkTrack()->GetRigidity();
188      if ( crig ) rig = L2->GetCaloLevel2()->qtot/260.;
189    //    //
190    // Fill the estrip matrix    // Fill the estrip matrix
191    //    //
# Line 119  void CaloFranzini::Process(Int_t itr){ Line 194  void CaloFranzini::Process(Int_t itr){
194    Int_t plane = 0;    Int_t plane = 0;
195    Int_t strip = 0;    Int_t strip = 0;
196    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;  
     //  
   };  
197    //    //
198    //    //
199    //    //
200    if ( dolong ){    if ( dolong ){
201      //      //
202      degfre = 44;      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++){      for (Int_t i=0; i < 22; i++){
241        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){
242          degfre = 2 * i;          dgf = 2 * i;
243          break;          break;
244        };        };
245        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){        if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){
246        degfre = 1 + 2 * i;        dgf = 1 + 2 * i;
247        break;        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 ){      if ( degfre > 0 ){
       TArrayF *qplmean = this->LoadLongAverage(rig);  
       TMatrixD *lmatrix = this->LoadCovarianceMatrix(rig);  
259        for (Int_t i = 0; i < degfre; i++){        for (Int_t i = 0; i < degfre; i++){
260          for (Int_t j = 0; j < degfre; j++){              if ( i != mask18b ){
261            longtzeta += (qplane[i]-qplmean->At(i)) * (*lmatrix)[i][j] * (qplane[j]-qplmean->At(j));            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 ){    if ( dofull ){
281      printf(" ERROR: NOT IMPLEMENTED YET\n");      //    printf(" ERROR: NOT IMPLEMENTED YET\n");
282      printf(" ERROR: CaloFranzini variables _NOT_ filled \n");      //    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();    //  if ( debug ) this->Print();
499    if ( debug ) printf(" exit \n");    if ( debug ) printf(" exit \n");
500  }  }
501    
# Line 169  void CaloFranzini::Process(Int_t itr){ Line 503  void CaloFranzini::Process(Int_t itr){
503  Float_t CaloFranzini::GetNormLongTZeta(){  Float_t CaloFranzini::GetNormLongTZeta(){
504    Process();    Process();
505    Float_t normz = 0.;    Float_t normz = 0.;
506    if ( degfre > 0 ) normz = longtzeta/(Float_t)degfre;    if ( degfre != 0 ) normz = longtzeta/(Float_t)degfre;
507    return normz;    return normz;
508  }  }
509    
510  Float_t CaloFranzini::GetNormFullTZeta(){  Float_t CaloFranzini::GetNormFullTZeta(){
511    Process();    Process();
512    Float_t normz = 0.;    Float_t normz = 0.;
513    if ( degfre > 0 ) normz = fulltzeta/(Float_t)degfre;    if ( fdegfre != 0 ) normz = fulltzeta/(Float_t)fdegfre;
514    return normz;    return normz;
515  }  }
516    
517  Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){    Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){  
518    //    //
519    file = new TFile(matrixfile.Data(),"READ");    lfile = new TFile(matrixfile.Data(),"READ");
520    //    //
521    if ( !file || file->IsZombie() ){    if ( !lfile || lfile->IsZombie() ){
522      file = new TFile(matrixfile.Data(),"RECREATE");      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 {    } else {
531        lfile->Close();
532      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());
533      return(false);      return(false);
534    };    };
# Line 195  Bool_t CaloFranzini::CreateMatrixFile(TS Line 537  Bool_t CaloFranzini::CreateMatrixFile(TS
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){  void CaloFranzini::WriteNumBin(Int_t numbin){
565    TArrayI nbi(1, &numbin);    if ( dolong ){
566    file->WriteObject(&nbi, "nbinenergy");      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){  void CaloFranzini::WriteRigBin(TArrayF *rigbin){
578    file->WriteObject(&rigbin, "binrig");    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){  void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){
591      lfile->cd();
592    TString name = Form("qplmeann%i",bin);    TString name = Form("qplmeann%i",bin);
593    file->WriteObject(&qpl,name.Data());    lfile->WriteObject(&(*qpl),name.Data());
594  }  }
595    
596  void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){  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);    TString name = Form("matrixn%i",bin);
613    file->WriteObject(&(*mat),name.Data());    //  mat.Write(name.Data());
614      lfile->WriteObject(&mat,name.Data());
615  }  }
616    
617  void CaloFranzini::CloseMatrixFile(){  void CaloFranzini::WriteInvertedFullMatrix(TMatrixD mat, Int_t bin){
618    file->Close();    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){    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    // find matrix file
665    //    //
666    if ( !strcmp(matrixfile.Data(),"") ){    if ( !strcmp(longmatrixfile.Data(),"") ){
667      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";    
668      if (dofull) matrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_fullel.root";        };
669      if ( !strcmp(fullmatrixfile.Data(),"") ){
670        if (dofull) fullmatrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_fullel.root";    
671    };    };
672    //    //
673    file = new TFile(matrixfile.Data(),"READ");    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 ( !file || file->IsZombie() ){    if ( dofull ){
699      printf(" ERROR: cannot open file %s \n",matrixfile.Data());      //
700      return(false);      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    TArrayI *numbin = (TArrayI*)file->Get("nbinenergy");    //
727    if ( !numbin ){    return(true);
728      printf(" ERROR: cannot read number of bins from file %s \n",matrixfile.Data());    //
729      return(false);  }
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    nbin = (Int_t)numbin->At(0);    //
789    if ( nbin <= 0 ){    return(true);
790      printf(" ERROR: cannot work with 0 energy bins (from file %s) \n",matrixfile.Data());  }
791      return(false);  
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    brig = (TArrayF*)file->Get("binrig");    return(true);
804    if ( !brig ){  }
805      printf(" ERROR: cannot read rigidity binning from file %s \n",matrixfile.Data());  
806      return(false);  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);    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){  TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){
841    //    //
842    TString name;    Int_t mv = 0;
843    for (Int_t i = 0; i<nbin-1; i++){    for (Int_t i = 0; i<nbin-1; i++){
844      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
845        name = Form("matrixn%i",i);        mv = i;
846        break;        break;
847      };      };
848    };    };
849    if ( rig < brig->At(0) ){    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));      printf(" WARNING: Event with rigidity lower than the first covariance matrix bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
851      name = "matrixn0";      mv = 0;
     printf(" Using matrix %s \n",name.Data());  
852    };    };
853    if ( rig >= brig->At(nbin) ){    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));      printf(" WARNING: Event with rigidity higher than the last covariance matrix bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
855      name = Form("matrixn%i",nbin-1);      mv = nbin-2;
     printf(" Using matrix %s \n",name.Data());  
856    };    };
857    //    //
858    TMatrixD *matrix = (TMatrixD*)file->Get(name.Data());    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(matrix);    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){  TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){
935    //    //
936    TString name;    Int_t mv=0;
937    for (Int_t i = 0; i<nbin-1; i++){    for (Int_t i = 0; i<nbin-1; i++){
938      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){      if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
939        name = Form("qplmeann%i",i);          mv = i;
940        break;        break;
941      };      };
942    };    };
943    if ( rig < brig->At(0) ){    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));      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
945      name = "qplmeann0";      mv = 0;
946      printf(" Using qplmean %s \n",name.Data());    };
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    if ( rig >= brig->At(nbin) ){    //
1017      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin));    //
1018      name = Form("qplmeann%i",nbin-1);    if ( rig < brigm->At(0) ){
1019      printf(" Using qplmean %s \n",name.Data());      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    TArrayF *qplmean = (TArrayF*)file->Get(name.Data());    if ( mstrip >= 101 && mstrip < 107 ) lastrip = 9;
1392    //    //
1393    return(qplmean);    if ( mstrip >= 107 ) lastrip = 10;
1394    //    //
1395      return(lastrip);
1396  }  }
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.23