/[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.7 by mocchiut, Mon Jan 21 10:24:11 2008 UTC revision 1.8 by mocchiut, Fri Jan 25 15:09:07 2008 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;    brigm = NULL;
23    nbin = 0;    nbin = 0;
# Line 52  void CaloFranzini::Clear(){ Line 53  void CaloFranzini::Clear(){
53    longtzeta = 0.;    longtzeta = 0.;
54    fulltzeta = 0.;    fulltzeta = 0.;
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, 43*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    
71  void CaloFranzini::Print(){  void CaloFranzini::Print(){
# Line 64  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(" degree of freedom                :.. %i\n",this->GetDegreeOfFreedom());      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());      printf(" longtzeta normalized             :.. %f\n",this->GetNormLongTZeta());  
81    //  printf(" fulltzeta                        :.. %f\n",fulltzeta);      printf(" full degree of freedom           :.. %i\n",this->GetFullDegreeOfFreedom());  
82    //  printf(" longtzeta normalized             :.. %f\n",this->GetNormFullTZeta());      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    //    //
# Line 120  void CaloFranzini::Process(Int_t itr){ Line 141  void CaloFranzini::Process(Int_t itr){
141      return;      return;
142    };    };
143    //    //
144    if ( !nbin || !file || !brig ){    if ( !nbin || !brig || (!lfile && !ffile) ){
145      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");
146      printf(" ERROR: CaloFranzini variables _NOT_ filled \n");      printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
147      return;      return;
# Line 158  void CaloFranzini::Process(Int_t itr){ Line 179  void CaloFranzini::Process(Int_t itr){
179    Int_t plane = 0;    Int_t plane = 0;
180    Int_t strip = 0;    Int_t strip = 0;
181    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 - N);  
     if ( nplane > (37-(2*N)) ) nplane--;  
     //  
 //    if ( plane == (18+N) ) mip = 0.;  
     if ( nplane > -1 ) qplane[nplane] += mip;  
     //  
   };  
182    //    //
183    //    //
184    //    //
185    if ( dolong ){    if ( dolong ){
186      //      //
187        for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
188          //
189          mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
190          //
191          estrip[view][plane][strip] = mip;
192          //    
193          nplane = 1 - view + 2 * (plane - N);
194          if ( nplane > (37-(2*N)) ) nplane--;
195          //
196          //    if ( plane == (18+N) ) mip = 0.;
197          if ( nplane > -1 ) qplane[nplane] += mip;
198          //
199        };
200        //
201      if ( cont ){      if ( cont ){
202        for (Int_t i=0; i<22; i++){        for (Int_t i=0; i<22; i++){
203          if ( i == (18+N) ){          if ( i == (18+N) ){
# Line 213  void CaloFranzini::Process(Int_t itr){ Line 235  void CaloFranzini::Process(Int_t itr){
235      //      //
236      degfre = TMath::Min(dgf,NC);      degfre = TMath::Min(dgf,NC);
237      //      //
238      Float_t longzdiag = 0.;      //    Float_t longzdiag = 0.;
239      Float_t longzout = 0.;      //    Float_t longzout = 0.;
240      //      //
241      if ( degfre > 0 ){      if ( degfre > 0 ){
242        for (Int_t i = 0; i < degfre; i++){        for (Int_t i = 0; i < degfre; i++){
243          if ( i != mask18b ){          if ( i != mask18b ){
244            for (Int_t j = 0; j < degfre; j++){              for (Int_t j = 0; j < degfre; j++){  
245              if ( j != mask18b ){              if ( j != mask18b ){
246                if ( i == j ){  //            if ( i == j ){
247                  longzdiag += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));  //              longzdiag += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
248                  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));  //              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));
249                } else {  //            } else {
250                  longzout += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));  //              longzout += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
251                  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));  //              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));
252                };  //            };
253                  //
254                longtzeta += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));                longtzeta += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
255                  //
256              };              };
257            };            };
258          };          };
259        };        };
260        if ( debug ) printf(" diagonal total %f out of diagonal total %f \n",longzdiag,longzout);        //      if ( debug ) printf(" diagonal total %f out of diagonal total %f \n",longzdiag,longzout);
261      };      };
262    };    };
263    if ( dofull ){    if ( dofull ){
264      printf(" ERROR: NOT IMPLEMENTED YET\n");      //    printf(" ERROR: NOT IMPLEMENTED YET\n");
265      printf(" ERROR: CaloFranzini variables _NOT_ filled \n");      //    printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
266    };        //
267        // FULL CALORIMETER
268        //
269        CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
270        //
271        Int_t dgf = 0;
272        Int_t cs = 0;
273        Int_t cd = 0;
274        Int_t mstrip = 0;
275        //
276        Int_t maxnpl = 0;
277        Float_t mipv[43][11];
278        memset(mipv,0,43*11*sizeof(Float_t));
279        //
280        for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
281          //
282          mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
283          //
284          //      nplane = 1 - view + 2 * plane;
285          //      if ( nplane > 37 ) nplane--;
286          nplane = 1 - view + 2 * (plane - N);
287          if ( nplane > (37-(2*N)) ) nplane--;
288          //
289          cs = ct->tibar[plane][view] - 1;
290          //
291          if ( ct->tibar[plane][view] > -1 ){
292            //
293            cd = 95 - cs;
294            //
295            mstrip = cd + strip;
296            //
297            Int_t lstr = this->ConvertStrip(mstrip);
298            //
299            if ( nplane > maxnpl ) maxnpl = nplane;
300            if ( nplane > -1 ) mipv[nplane][lstr] += mip;
301            //
302          };
303          //      mipv[nplane][lstr] += mip;
304          //
305        };
306        //
307        Float_t mip1 = 1.;
308        Int_t cs1;
309        Int_t cd1;
310        Float_t mip2 = 1.;
311        Int_t cs2;
312        Int_t cd2;
313        Int_t mi = -1;
314        Int_t mj = -1;
315        Int_t nn1 = 0;
316        Int_t pl1 = 0;
317        Int_t vi1 = 0;
318        Int_t nn2 = 0;
319        Int_t pl2 = 0;
320        Int_t vi2 = 0;
321        Int_t mstrip1min = 0;
322        Int_t mstrip1max = 0;
323        Int_t mstrip2min = 0;
324        Int_t mstrip2max = 0;
325        //
326        Int_t tonpl = TMath::Min(maxnpl,NC);
327        //
328        Int_t rbi = 0;
329        for (Int_t i = 0; i<nbin-1; i++){
330          if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
331            rbi = i;
332            break;
333          };
334        };
335        //
336        //
337        Int_t therigb = rbi;
338        //
339        if ( rig < brigm->At(2) ){
340    //      therigb = 0;
341          therigb = 2;
342        };
343        //    if ( rig >= brigm->At(nbin-5) ){
344        if ( rig >= brigm->At(nbin-2) ){
345          therigb = nbin - 3;
346          //      therigb = nbin - 5;
347        };
348        Int_t mtherig = therigb;
349        if ( mtherig >= 13 ) mtherig = 12;
350        //
351        if ( debug ) printf(" rig %f brigm0 %f brigm n2 %f nbin %i \n",rig,brigm->At(0),brigm->At(nbin-2),nbin);
352        //
353        if ( cont ){
354          for (Int_t i=0; i<22; i++){
355            if ( i == (18+N) ){
356              mask18b =  1 + 2 * (i - N);
357              break;
358            };
359          };
360        };  
361        //  
362        if ( sel ){
363          for (Int_t i=0; i<22; i++){
364            if ( i == (18-N) ){
365              mask18b =  1 + 2 * (i - N);
366              break;
367            };
368          };
369        };
370        //
371        if ( mask18b == 37 ) mask18b = -1;
372        //
373        if ( debug ) printf("Start main loop \n");
374        //
375        for (Int_t nplane1=0; nplane1<tonpl; nplane1++){
376          //
377          if ( nplane1 != mask18b ){
378            //
379            if ( nplane1 >= 37 ) nn1 = nplane1 + 1;
380            vi1 = 1;
381            if ( nn1%2 ) vi1 = 0;
382            pl1 = (nn1 - 1 + vi1)/2;
383            //
384            cs1 = ct->tibar[pl1][vi1] - 1;
385            //
386            if ( ct->tibar[pl1][vi1] > -1 ){
387              //
388              dgf++;
389              //
390              cd1 = 95 - cs1;
391              //
392              Int_t at1 = TMath::Max(0,(cd1+0));
393              Int_t at2 = TMath::Min(190,(cd1+95));
394              mstrip1min = this->ConvertStrip(at1);
395              mstrip1max = this->ConvertStrip(at2) + 1;
396              //
397              for (Int_t mstrip1=mstrip1min; mstrip1<mstrip1max; mstrip1++){
398                //
399                mj = -1;
400                //
401                //      if ( debug ) printf(" get mip1 %i %i %f %i \n",nplane1,mstrip1,rig,therigb);
402                mip1 = mipv[nplane1][mstrip1] - this->GetFullAverageAt(nplane1,mstrip1,rig,therigb);
403                if ( mip1 < 0 ){
404                  numneg++;
405                  avenegvar += mip1;
406                };
407                if ( mip1 >= 0 ){
408                  numpos++;
409                  aveposvar += mip1;
410                };
411                if ( mip1 < minsvalue ) minsvalue = mip1;
412                if ( mip1 >= maxsvalue ) maxsvalue = mip1;
413                //
414                mi = (nplane1 * 11) + mstrip1;
415                //
416                if ( mip1 != 0. ){
417                  //
418                  for (Int_t nplane2=0; nplane2<tonpl; nplane2++){
419                    //
420                    if ( nplane2 != mask18b ){
421                      //
422                      if ( nplane2 >= 37 ) nn2 = nplane2 + 1;
423                      vi2 = 1;
424                      if ( nn2%2 ) vi2 = 0;
425                      pl1 = (nn2 - 1 + vi2)/2;
426                      //
427                      cs2 = ct->tibar[pl2][vi2] - 1;
428                      //
429                      if ( ct->tibar[pl2][vi2] > -1 ){
430                        //
431                        cd2 = 95 - cs2;
432                        //
433                        Int_t t1 = TMath::Max(0,(cd2+0));
434                        Int_t t2 = TMath::Min(190,(cd2+95));
435                        mstrip2min = this->ConvertStrip(t1);
436                        mstrip2max = this->ConvertStrip(t2) + 1;
437                        //
438                        for (Int_t mstrip2=mstrip2min; mstrip2<mstrip2max; mstrip2++){
439                          //              
440                          mip2 = mipv[nplane2][mstrip2] - this->GetFullAverageAt(nplane2,mstrip2,rig,therigb);
441                          //
442                          if ( mip2 != 0. ){
443                            //
444    //                      Int_t sh = -5 + nplane1;
445    //                      if ( sh > 5 ) sh -= 11*nplane1;
446    //                      //
447    //                      mj = (nplane2 * 11) + mstrip2 + sh;
448    //                      //
449    //                      if ( mj < 0 ) mj += 473;
450    //                      if ( mj >= 473 ) mj -= 473;
451    //                      //
452                            mj = (nplane2 * 11) + mstrip2;
453                            //
454                            Float_t svalue =  mip1 * this->GetFullHmatrixAt(mi,mj,rig,mtherig,therigb) * mip2;
455                            //                      fulltzeta += mip1 * this->GetFullHmatrixAt(mi,mj,rig,therigb) * mip2;
456                            fulltzeta += svalue;
457                            if ( svalue < 0. ) negfulltzeta += svalue;
458                            if ( svalue > 0. ) posfulltzeta += svalue;
459                            //              (*fmatrix[rbi])[mi][mj] += (mip1 * mip2);
460                            //                      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);
461                            //
462                          };
463                        };
464                      };
465                    };
466                  };
467                };
468              };
469            };
470          };
471        };
472        //
473        fdegfre = dgf*11;
474        if ( numpos ) aveposvar /= numpos;
475        if ( numneg ) avenegvar /= numneg;
476        //
477      };
478    //    //
479  //  if ( debug ) this->Print();    //  if ( debug ) this->Print();
480    if ( debug ) printf(" exit \n");    if ( debug ) printf(" exit \n");
481  }  }
482    
# Line 249  void CaloFranzini::Process(Int_t itr){ Line 484  void CaloFranzini::Process(Int_t itr){
484  Float_t CaloFranzini::GetNormLongTZeta(){  Float_t CaloFranzini::GetNormLongTZeta(){
485    Process();    Process();
486    Float_t normz = 0.;    Float_t normz = 0.;
487    if ( degfre > 0 ) normz = longtzeta/(Float_t)degfre;    if ( degfre != 0 ) normz = longtzeta/(Float_t)degfre;
488    return normz;    return normz;
489  }  }
490    
491  Float_t CaloFranzini::GetNormFullTZeta(){  Float_t CaloFranzini::GetNormFullTZeta(){
492    Process();    Process();
493    Float_t normz = 0.;    Float_t normz = 0.;
494    if ( degfre > 0 ) normz = fulltzeta/(Float_t)degfre;    if ( fdegfre != 0 ) normz = fulltzeta/(Float_t)fdegfre;
495    return normz;    return normz;
496  }  }
497    
498  Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){    Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){  
499    //    //
500    file = new TFile(matrixfile.Data(),"READ");    lfile = new TFile(matrixfile.Data(),"READ");
501    //    //
502    if ( !file || file->IsZombie() ){    if ( !lfile || lfile->IsZombie() ){
503      file = new TFile(matrixfile.Data(),"RECREATE");      if ( dolong ){
504      printf(" Create file %s \n",file->GetName());        lfile = new TFile(matrixfile.Data(),"RECREATE");
505          printf(" Create file %s \n",lfile->GetName());
506        };
507        if ( dofull ){
508          ffile = new TFile(matrixfile.Data(),"RECREATE");
509          printf(" Create file %s \n",ffile->GetName());
510        };
511    } else {    } else {
512        lfile->Close();
513      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());
514      return(false);      return(false);
515    };    };
# Line 278  Bool_t CaloFranzini::CreateMatrixFile(TS Line 520  Bool_t CaloFranzini::CreateMatrixFile(TS
520    
521  Bool_t CaloFranzini::UpdateMatrixFile(TString matrixfile){    Bool_t CaloFranzini::UpdateMatrixFile(TString matrixfile){  
522    //    //
523    file = new TFile(matrixfile.Data(),"UPDATE");    if ( dolong ){
524        lfile = new TFile(matrixfile.Data(),"UPDATE");
525        //
526        if ( !lfile || lfile->IsZombie() ){
527          printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
528          return(false);
529        };
530      };
531    //    //
532    if ( !file || file->IsZombie() ){    if ( dofull ){
533      printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());      ffile = new TFile(matrixfile.Data(),"UPDATE");
534      return(false);      //
535        if ( !ffile || ffile->IsZombie() ){
536          printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
537          return(false);
538        };
539    };    };
540    //    //
541    return(true);    return(true);
# Line 290  Bool_t CaloFranzini::UpdateMatrixFile(TS Line 543  Bool_t CaloFranzini::UpdateMatrixFile(TS
543  }  }
544    
545  void CaloFranzini::WriteNumBin(Int_t numbin){  void CaloFranzini::WriteNumBin(Int_t numbin){
546    file->cd();    if ( dolong ){
547    TArrayI nbi(1, &numbin);      lfile->cd();
548    file->WriteObject(&nbi, "nbinenergy");      TArrayI nbi(1, &numbin);
549        lfile->WriteObject(&nbi, "nbinenergy");
550      };
551      if ( dofull ){
552        ffile->cd();
553        TArrayI nbi(1, &numbin);
554        ffile->WriteObject(&nbi, "nbinenergy");
555      };
556  }  }
557    
558  void CaloFranzini::WriteRigBin(TArrayF *rigbin){  void CaloFranzini::WriteRigBin(TArrayF *rigbin){
559    file->cd();    if ( dolong ){
560    //  rigbin->Write("binrig");      lfile->cd();
561    file->WriteObject(&(*rigbin), "binrig");      //  rigbin->Write("binrig");
562        lfile->WriteObject(&(*rigbin), "binrig");
563      };
564      if ( dofull ){
565        ffile->cd();
566        //  rigbin->Write("binrig");
567        ffile->WriteObject(&(*rigbin), "binrig");
568      };
569  }  }
570    
571  void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){  void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){
572    file->cd();    lfile->cd();
573    TString name = Form("qplmeann%i",bin);    TString name = Form("qplmeann%i",bin);
574    file->WriteObject(&(*qpl),name.Data());    lfile->WriteObject(&(*qpl),name.Data());
575  }  }
576    
577  void CaloFranzini::WriteFullMean(TMatrixD *qpl, Int_t bin){  void CaloFranzini::WriteFullMean(TMatrixD *qpl, Int_t bin){
578    file->cd();    ffile->cd();
579    TString name = Form("fqplmeann%i",bin);    TString name = Form("fqplmeann%i",bin);
580    file->WriteObject(&(*qpl),name.Data());    ffile->WriteObject(&(*qpl),name.Data());
581    file->Purge();    ffile->Purge();
582  }  }
583    
584  void CaloFranzini::WriteFullNMean(TMatrixD *qpl, Int_t bin){  void CaloFranzini::WriteFullNMean(TMatrixD *qpl, Int_t bin){
585    file->cd();    ffile->cd();
586    TString name = Form("fnqplmeann%i",bin);    TString name = Form("fnqplmeann%i",bin);
587    file->WriteObject(&(*qpl),name.Data());    ffile->WriteObject(&(*qpl),name.Data());
588    file->Purge();    ffile->Purge();
589  }  }
590    
591  void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){  void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){
592    file->cd();    lfile->cd();
593    TString name = Form("matrixn%i",bin);    TString name = Form("matrixn%i",bin);
594    //  mat.Write(name.Data());    //  mat.Write(name.Data());
595    file->WriteObject(&mat,name.Data());    lfile->WriteObject(&mat,name.Data());
596  }  }
597    
598  void CaloFranzini::WriteInvertedFullMatrix(TMatrixD mat, Int_t bin){  void CaloFranzini::WriteInvertedFullMatrix(TMatrixD mat, Int_t bin){
599    file->cd();    ffile->cd();
600    TString name = Form("fmatrixn%i",bin);    TString name = Form("fmatrixn%i",bin);
601    //  mat.Write(name.Data());    //  mat.Write(name.Data());
602    file->WriteObject(&mat,name.Data());    ffile->WriteObject(&mat,name.Data());
603  }  }
604    
605  void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){  void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
606    file->cd();    lfile->cd();
607    TString name = Form("origmatrixn%i",bin);    TString name = Form("origmatrixn%i",bin);
608    //  mat.Write(name.Data());    //  mat.Write(name.Data());
609    file->WriteObject(&(*mat),name.Data());    lfile->WriteObject(&(*mat),name.Data());
610  }  }
611    
612  void CaloFranzini::WriteFullMatrix(TMatrixD *mat, Int_t bin){  void CaloFranzini::WriteFullMatrix(TMatrixD *mat, Int_t bin){
613    file->cd();    ffile->cd();
614    TString name = Form("origfmatrixn%i",bin);    TString name = Form("origfmatrixn%i",bin);
615    //  mat.Write(name.Data());    //  mat.Write(name.Data());
616    file->WriteObject(&(*mat),name.Data());    ffile->WriteObject(&(*mat),name.Data());
617    file->Purge();    ffile->Purge();
618  }  }
619    
620  void CaloFranzini::WriteFullNMatrix(TMatrixF *mat, Int_t bin){  void CaloFranzini::WriteFullNMatrix(TMatrixF *mat, Int_t bin){
621    file->cd();    ffile->cd();
622    TString name = Form("origfnmatrixn%i",bin);    TString name = Form("origfnmatrixn%i",bin);
623    //  mat.Write(name.Data());    //  mat.Write(name.Data());
624    file->WriteObject(&(*mat),name.Data());    ffile->WriteObject(&(*mat),name.Data());
625    file->Purge();    ffile->Purge();
626  }  }
627    
628  void CaloFranzini::CloseMatrixFile(){  void CaloFranzini::CloseMatrixFile(){
629    file->cd();    if ( dolong ){
630    file->Close();      lfile->cd();
631        lfile->Close();
632      };
633      if ( dofull ){
634        ffile->cd();
635        ffile->Close();
636      };
637  }  }
638    
   
639  Bool_t CaloFranzini::Open(TString matrixfile){    Bool_t CaloFranzini::Open(TString matrixfile){  
640      return(this->Open(matrixfile,""));
641    }
642    
643    Bool_t CaloFranzini::Open(TString longmatrixfile, TString fullmatrixfile){  
644    //    //
645    // find matrix file    // find matrix file
646    //    //
647    if ( !strcmp(matrixfile.Data(),"") ){    if ( !strcmp(longmatrixfile.Data(),"") ){
648      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";    
     if (dofull) matrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_fullel.root";      
   };  
   //  
   file = new TFile(matrixfile.Data(),"READ");  
   //  
   if ( !file || file->IsZombie() ){  
     printf(" ERROR: cannot open file %s \n",matrixfile.Data());  
     return(false);  
649    };    };
650    //    if ( !strcmp(fullmatrixfile.Data(),"") ){
651    if ( !this->LoadBin() ){      if (dofull) fullmatrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_fullel.root";    
     printf(" %s \n",matrixfile.Data());  
     return(false);  
652    };    };
653    //    //
654    if ( dolong ){    if ( dolong ){
655        //
656        lfile = new TFile(longmatrixfile.Data(),"READ");
657        //
658        if ( !lfile || lfile->IsZombie() ){
659          printf(" ERROR: cannot open file %s \n",longmatrixfile.Data());
660          return(false);
661        };
662        //
663        if ( !this->LoadBin() ){
664          printf(" %s \n",longmatrixfile.Data());
665          return(false);
666        };
667        //
668      if ( !this->LoadLong() ){      if ( !this->LoadLong() ){
669        printf(" %s \n",matrixfile.Data());        printf(" %s \n",longmatrixfile.Data());
670        return(false);        return(false);
671      };      };
672      //      //
673      if ( !this->LoadMatrices() ){      if ( !this->LoadMatrices() ){
674        printf(" %s \n",matrixfile.Data());        printf(" %s \n",longmatrixfile.Data());
675        return(false);        return(false);
676      };      };
677    };    };
678    //    //
679    if ( dofull ){    if ( dofull ){
680        //
681        ffile = new TFile(fullmatrixfile.Data(),"READ");
682        //
683        if ( !ffile || ffile->IsZombie() ){
684          printf(" ERROR: cannot open file %s \n",fullmatrixfile.Data());
685          return(false);
686        };
687        //
688        if ( !dolong ){
689          //
690          if ( !this->LoadBin(true) ){
691            printf(" %s \n",fullmatrixfile.Data());
692            return(false);
693          };      
694          //
695        };
696      if ( !this->LoadFull() ){      if ( !this->LoadFull() ){
697        printf(" %s \n",matrixfile.Data());        printf(" %s \n",fullmatrixfile.Data());
698        return(false);        return(false);
699      };      };
700      //      //
701      if ( !this->LoadFullMatrices() ){      if ( !this->LoadFullMatrices() ){
702        printf(" %s \n",matrixfile.Data());        printf(" %s \n",fullmatrixfile.Data());
703        return(false);        return(false);
704      };      };
705    };    };
# Line 414  Bool_t CaloFranzini::Open(TString matrix Line 709  Bool_t CaloFranzini::Open(TString matrix
709    //    //
710  }  }
711    
712    
713  Bool_t CaloFranzini::LoadBin(){  Bool_t CaloFranzini::LoadBin(){
714      return(this->LoadBin(false));
715    }
716    
717    Bool_t CaloFranzini::LoadBin(Bool_t full){
718    //    //
719    TArrayI *numbin = (TArrayI*)file->Get("nbinenergy");    if ( !full ) {
720    if ( !numbin ){      //
721      printf(" ERROR: cannot read number of bins from file ");      TArrayI *numbin = (TArrayI*)lfile->Get("nbinenergy");
722      return(false);      if ( !numbin ){
723    };        printf(" ERROR: cannot read number of bins from file ");
724    nbin = (Int_t)numbin->At(0);        return(false);
725    if ( nbin <= 0 ){      };
726      printf(" ERROR: cannot work with 0 energy bins from file ");      nbin = (Int_t)numbin->At(0);
727      return(false);      if ( nbin <= 0 ){
728    };        printf(" ERROR: cannot work with 0 energy bins from file ");
729    //        return(false);
730    brig = (TArrayF*)file->Get("binrig");      };
731    if ( !brig ){      //
732      printf(" ERROR: cannot read rigidity binning from file ");      brig = (TArrayF*)lfile->Get("binrig");
733      return(false);      if ( !brig ){
734    };        printf(" ERROR: cannot read rigidity binning from file ");
735    //        return(false);
736    brigm=(TArrayF*)file->Get("binrigmean");      };
737    if ( !brigm ){      //
738      printf(" ERROR: cannot read mean rigidity binning from file ");      brigm=(TArrayF*)lfile->Get("binrigmean");
739      return(false);      if ( !brigm ){
740          printf(" ERROR: cannot read mean rigidity binning from file ");
741          return(false);
742        };
743        //
744      } else {
745        //
746        TArrayI *numbin = (TArrayI*)ffile->Get("nbinenergy");
747        if ( !numbin ){
748          printf(" ERROR: cannot read number of bins from file ");
749          return(false);
750        };
751        nbin = (Int_t)numbin->At(0);
752        if ( nbin <= 0 ){
753          printf(" ERROR: cannot work with 0 energy bins from file ");
754          return(false);
755        };
756        //
757        brig = (TArrayF*)ffile->Get("binrig");
758        if ( !brig ){
759          printf(" ERROR: cannot read rigidity binning from file ");
760          return(false);
761        };
762        //
763        brigm=(TArrayF*)ffile->Get("binrigmean");
764        if ( !brigm ){
765          printf(" ERROR: cannot read mean rigidity binning from file ");
766          return(false);
767        };
768    };    };
769    //    //
770    return(true);    return(true);
# Line 446  Bool_t CaloFranzini::LoadLong(){ Line 774  Bool_t CaloFranzini::LoadLong(){
774    //    //
775    for (Int_t i=0;i<17;i++){    for (Int_t i=0;i<17;i++){
776      TString name = Form("qplmeann%i",i);      TString name = Form("qplmeann%i",i);
777      qplmean[i] = (TArrayF*)file->Get(name.Data());      qplmean[i] = (TArrayF*)lfile->Get(name.Data());
778      if ( !qplmean[i] ){      if ( !qplmean[i] ){
779        printf(" ERROR: cannot read average from file ");        printf(" ERROR: cannot read average from file ");
780        return(false);        return(false);
# Line 460  Bool_t CaloFranzini::LoadFull(){ Line 788  Bool_t CaloFranzini::LoadFull(){
788    //    //
789    for (Int_t i=0;i<17;i++){    for (Int_t i=0;i<17;i++){
790      TString name = Form("fqplmeann%i",i);      TString name = Form("fqplmeann%i",i);
791      fqplmean[i] = (TMatrixD*)file->Get(name.Data());      fqplmean[i] = (TMatrixD*)ffile->Get(name.Data());
792      if ( !fqplmean[i] ){      if ( !fqplmean[i] ){
793        printf(" ERROR: cannot read average from file ");        printf(" ERROR: cannot read average from file ");
794        return(false);        return(false);
# Line 474  Bool_t CaloFranzini::LoadMatrices(){ Line 802  Bool_t CaloFranzini::LoadMatrices(){
802    //    //
803    for (Int_t i=0;i<17;i++){    for (Int_t i=0;i<17;i++){
804          TString name1 = Form("matrixn%i",i);          TString name1 = Form("matrixn%i",i);
805          hmat[i] = (TMatrixD*)file->Get(name1.Data());          hmat[i] = (TMatrixD*)lfile->Get(name1.Data());
806    };    };
807    //    //
808    return(true);    return(true);
# Line 484  Bool_t CaloFranzini::LoadFullMatrices(){ Line 812  Bool_t CaloFranzini::LoadFullMatrices(){
812    //    //
813    for (Int_t i=0;i<17;i++){    for (Int_t i=0;i<17;i++){
814          TString name1 = Form("fmatrixn%i",i);          TString name1 = Form("fmatrixn%i",i);
815          hfmat[i] = (TMatrixF*)file->Get(name1.Data());          hfmat[i] = (TMatrixD*)ffile->Get(name1.Data());
816    };    };
817    //    //
818    return(true);    return(true);
# Line 515  TMatrixD *CaloFranzini::LoadCovarianceMa Line 843  TMatrixD *CaloFranzini::LoadCovarianceMa
843  TMatrixD *CaloFranzini::LoadFullAverage(Int_t rigbin){  TMatrixD *CaloFranzini::LoadFullAverage(Int_t rigbin){
844    //    //
845    TString name = Form("fqplmeann%i",rigbin);    TString name = Form("fqplmeann%i",rigbin);
846    TMatrixD *fmean=(TMatrixD*)file->Get(name.Data());    TMatrixD *fmean=(TMatrixD*)ffile->Get(name.Data());
847    //    //
848    return(fmean);    return(fmean);
849    //    //
# Line 524  TMatrixD *CaloFranzini::LoadFullAverage( Line 852  TMatrixD *CaloFranzini::LoadFullAverage(
852  TMatrixF *CaloFranzini::LoadFullMatrix(Int_t rigbin){  TMatrixF *CaloFranzini::LoadFullMatrix(Int_t rigbin){
853    //    //
854    TString name = Form("origfmatrixn%i",rigbin);    TString name = Form("origfmatrixn%i",rigbin);
855    TMatrixF *fmatri=(TMatrixF*)file->Get(name.Data());    TMatrixF *fmatri=(TMatrixF*)ffile->Get(name.Data());
856    //    //
857    return(fmatri);    return(fmatri);
858    //    //
# Line 533  TMatrixF *CaloFranzini::LoadFullMatrix(I Line 861  TMatrixF *CaloFranzini::LoadFullMatrix(I
861  void CaloFranzini::LoadFullMatrix(Int_t rigbin, TMatrixF *&fmatri){  void CaloFranzini::LoadFullMatrix(Int_t rigbin, TMatrixF *&fmatri){
862    //    //
863    TString name = Form("origfmatrixn%i",rigbin);    TString name = Form("origfmatrixn%i",rigbin);
864    fmatri=(TMatrixF*)file->Get(name.Data());    fmatri=(TMatrixF*)ffile->Get(name.Data());
865    //    //
866  }  }
867    
868  void CaloFranzini::UnLoadFullAverage(Int_t rigbin){  void CaloFranzini::UnLoadFullAverage(Int_t rigbin){
869    //    //
870    TString name = Form("fqplmeann%i",rigbin);    TString name = Form("fqplmeann%i",rigbin);
871    file->Delete(name.Data());    ffile->Delete(name.Data());
872    //    //
873  }  }
874    
875  void CaloFranzini::UnLoadFullMatrix(Int_t rigbin){  void CaloFranzini::UnLoadFullMatrix(Int_t rigbin){
876    //    //
877    TString name = Form("origfmatrixn%i",rigbin);    TString name = Form("origfmatrixn%i",rigbin);
878    file->Delete(name.Data());    ffile->Delete(name.Data());
879    //    //
880  }  }
881    
882  TMatrixF *CaloFranzini::LoadFullNMatrix(Int_t rigbin){  TMatrixF *CaloFranzini::LoadFullNMatrix(Int_t rigbin){
883    //    //
884    TString name = Form("origfnmatrixn%i",rigbin);    TString name = Form("origfnmatrixn%i",rigbin);
885    TMatrixF *fnmatri=(TMatrixF*)file->Get(name.Data());    TMatrixF *fnmatri=(TMatrixF*)ffile->Get(name.Data());
886    //    //
887    return(fnmatri);    return(fnmatri);
888    //    //
# Line 563  TMatrixF *CaloFranzini::LoadFullNMatrix( Line 891  TMatrixF *CaloFranzini::LoadFullNMatrix(
891  void CaloFranzini::UnLoadFullNMatrix(Int_t rigbin){  void CaloFranzini::UnLoadFullNMatrix(Int_t rigbin){
892    //    //
893    TString name = Form("origfnmatrixn%i",rigbin);    TString name = Form("origfnmatrixn%i",rigbin);
894    file->Delete(name.Data());    ffile->Delete(name.Data());
895    //    //
896  }  }
897    
898  TMatrixD *CaloFranzini::LoadFullNAverage(Int_t rigbin){  TMatrixD *CaloFranzini::LoadFullNAverage(Int_t rigbin){
899    //    //
900    TString name = Form("fnqplmeann%i",rigbin);    TString name = Form("fnqplmeann%i",rigbin);
901    TMatrixD *fnmean=(TMatrixD*)file->Get(name.Data());    TMatrixD *fnmean=(TMatrixD*)ffile->Get(name.Data());
902    //    //
903    return(fnmean);    return(fnmean);
904    //    //
# Line 579  TMatrixD *CaloFranzini::LoadFullNAverage Line 907  TMatrixD *CaloFranzini::LoadFullNAverage
907  void CaloFranzini::UnLoadFullNAverage(Int_t rigbin){  void CaloFranzini::UnLoadFullNAverage(Int_t rigbin){
908    //    //
909    TString name = Form("fnqplmeann%i",rigbin);    TString name = Form("fnqplmeann%i",rigbin);
910    file->Delete(name.Data());    ffile->Delete(name.Data());
911    //    //
912  }  }
913    
# Line 706  Float_t CaloFranzini::GetFullAverageAt(I Line 1034  Float_t CaloFranzini::GetFullAverageAt(I
1034    if ( a == 0. || minene == 0. || ave != ave ){    if ( a == 0. || minene == 0. || ave != ave ){
1035      Float_t m = (maxene-minene)/(maxrig-minrig);      Float_t m = (maxene-minene)/(maxrig-minrig);
1036      Float_t q = minene - m * minrig;      Float_t q = minene - m * minrig;
1037      ave = rig * m + q;      ave = rig * m + q;
1038    };    };
1039    //      //  
1040    //  ave += (44.-plane)*strip;    //  ave += (44.-plane)*strip;
# Line 731  Float_t CaloFranzini::GetHmatrixAt(Int_t Line 1059  Float_t CaloFranzini::GetHmatrixAt(Int_t
1059    //    //
1060    if ( rig < brigm->At(0) ){    if ( rig < brigm->At(0) ){
1061      if ( rig < brig->At(0) ){      if ( rig < brig->At(0) ){
1062  //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));        //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
1063      };      };
1064      therigb = 0;      therigb = 0;
1065    };    };
1066  //  if ( rig >= brigm->At(nbin-4) ){ // -2    //  if ( rig >= brigm->At(nbin-4) ){ // -2
1067    if ( rig >= brigm->At(nbin-2) ){    if ( rig >= brigm->At(nbin-2) ){
1068      if ( rig >= brig->At(nbin-2) ) {      if ( rig >= brig->At(nbin-2) ) {
1069  //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));        //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
1070      };      };
1071  //    therigb = nbin - 5;// -3      //    therigb = nbin - 5;// -3
1072      therigb = nbin - 3;      therigb = nbin - 3;
1073    };    };
1074    //    //
# Line 785  Float_t CaloFranzini::GetFullHmatrixAt(I Line 1113  Float_t CaloFranzini::GetFullHmatrixAt(I
1113    //    //
1114    if ( rig < brigm->At(0) ){    if ( rig < brigm->At(0) ){
1115      if ( rig < brig->At(0) ){      if ( rig < brig->At(0) ){
1116  //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));        //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
1117      };      };
1118      therigb = 0;      therigb = 0;
1119    };    };
1120  //  if ( rig >= brigm->At(nbin-4) ){ // -2  //  if ( rig >= brigm->At(nbin-4) ){ // -2
1121    if ( rig >= brigm->At(nbin-2) ){    if ( rig >= brigm->At(nbin-2) ){
1122        //if ( rig >= brigm->At(nbin-5) ){
1123      if ( rig >= brig->At(nbin-2) ) {      if ( rig >= brig->At(nbin-2) ) {
1124  //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));        //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
1125      };      };
1126  //    therigb = nbin - 5;// -3      //    therigb = nbin - 5;// -3
1127        //  therigb = nbin - 5;// -3
1128      therigb = nbin - 3;      therigb = nbin - 3;
1129    };    };
1130    //    //
1131    if ( therigb < 2 ) therigb = 2;    if ( therigb < 2 ) therigb = 2;
1132      if ( therigb > 13 ) therigb = 13;
1133    //    //
1134    return(this->GetFullHmatrixAt(iindex,jindex,rig,therigb));    return(this->GetFullHmatrixAt(iindex,jindex,rig,therigb));
1135    //    //
1136  }  }
1137    
1138    
1139  Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig, Int_t therigb){  Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig, Int_t therigb){
1140    //    //
1141      return(this->GetFullHmatrixAt(iindex,jindex,rig,therigb,0));
1142      //
1143    }
1144    
1145    Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig, Int_t therigb, Int_t mtherig){  
1146      //
1147      //  Int_t lofit = 12;
1148      Int_t lofit = 3;
1149      //
1150      Float_t lowrig = 0.;
1151    Float_t minrig;    Float_t minrig;
1152    Float_t maxrig;    Float_t maxrig;
1153    //    //
1154      if ( mtherig > lofit ) lowrig = brigm->At(therigb-1);
1155    minrig = brigm->At(therigb);    minrig = brigm->At(therigb);
1156    maxrig = brigm->At(therigb+1);    maxrig = brigm->At(therigb+1);
1157    //  printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);    //if ( therigb > 10 )  printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);
1158    //    //
1159      Float_t lowene = 0.;
1160      if ( mtherig > lofit ) lowene = (*hfmat[therigb-1])[iindex][jindex];
1161    Float_t minene = (*hfmat[therigb])[iindex][jindex];    Float_t minene = (*hfmat[therigb])[iindex][jindex];
1162    Float_t maxene = (*hfmat[therigb+1])[iindex][jindex];    Float_t maxene = (*hfmat[therigb+1])[iindex][jindex];
1163    //  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);    //  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);
1164    //    //
1165  //  Float_t a = 0.;    Float_t ave = 0.;
 //  Float_t b = 0.;  
 //  Float_t ave = 0.;  
 //  if ( minene == 0. ){  
 //        
 //  } else {  
 //      b = log(maxene/minene)/(maxrig-minrig);  
 //      a = minene/exp(minrig*b);  
 //      ave = a*exp(b*rig);  
 //  };  
1166    //    //
1167    Float_t m = (maxene-minene)/(maxrig-minrig);    //  Float_t a = 0.;
1168    Float_t q = minene - m * minrig;    //  Float_t b = 0.;
1169    Float_t ave = rig * m + q;    //  Float_t ave = 0.;
1170    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);    //  if ( minene == 0. ){
1171      //    
1172      //  } else {
1173      //    b = log(maxene/minene)/(maxrig-minrig);
1174      //    a = minene/exp(minrig*b);
1175      //    ave = a*exp(b*rig);
1176      //  };
1177      //
1178      if ( mtherig > lofit ){
1179        //
1180        // FIT
1181        //
1182    //     Float_t x[3]={lowrig,minrig,maxrig};
1183    //     Float_t y[3]={lowene,minene,maxene};
1184    //     //
1185    //     TGraph *tg= new TGraph(3,x,y);
1186    //     //
1187    //     //     gStyle->SetLabelSize(0.04);
1188    //     //     gStyle->SetNdivisions(510,"XY");
1189    //     //     TCanvas *c = new TCanvas();
1190    //     //     c->Draw();
1191    //     //     c->cd();
1192    //     //     tg->Draw("AP");
1193    //     //
1194    //     TF1 *fun = new TF1("fun","pol2");
1195    //     tg->Fit("fun","QNC");
1196    //     //
1197    //     ave = fun->GetParameter(0) + rig * fun->GetParameter(1) + rig * rig * fun->GetParameter(2);
1198    //     //
1199    //     //    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);
1200    //     //
1201    //     //
1202    //     //     c->Modified();
1203    //     //     c->Update();
1204    //     //     gSystem->ProcessEvents();
1205    //     //     gSystem->Sleep(6000);
1206    //     //     c->Close();
1207    //     //     gStyle->SetLabelSize(0);
1208    //     //     gStyle->SetNdivisions(1,"XY");
1209    //     //
1210    //    tg->Delete();
1211    //    fun->Delete();
1212        Float_t mrigl = (lowrig+minrig)/2.;
1213        Float_t mrigu = (minrig+maxrig)/2.;
1214        Float_t menel = (lowene+minene)/2.;
1215        Float_t meneu = (minene+maxene)/2.;
1216        Float_t m = (meneu-menel)/(mrigu-mrigl);
1217        Float_t q = menel - m * mrigl;
1218        ave = rig * m + q;
1219        //
1220      } else {
1221        Float_t m = (maxene-minene)/(maxrig-minrig);
1222        Float_t q = minene - m * minrig;
1223        ave = rig * m + q;
1224        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);
1225      };
1226    //    //
1227    //    //
1228    return(ave);    return(ave);
# Line 927  void CaloFranzini::DrawLongAverage(Int_t Line 1318  void CaloFranzini::DrawLongAverage(Int_t
1318  Int_t CaloFranzini::ConvertStrip(Int_t mstrip){    Int_t CaloFranzini::ConvertStrip(Int_t mstrip){  
1319    //    //
1320    Int_t lastrip = 0;    Int_t lastrip = 0;
1321    //  //   //
1322    if ( mstrip < 50 ) lastrip = 0;  //   if ( mstrip < 50 ) lastrip = 0;
1323    //  //   //
1324    if ( mstrip >= 50 && mstrip < 64 ) lastrip = 1;  //   if ( mstrip >= 50 && mstrip < 64 ) lastrip = 1;
1325    //  //   //
1326    if ( mstrip >= 64 && mstrip < 70 ) lastrip = 2;  //   if ( mstrip >= 64 && mstrip < 70 ) lastrip = 2;
1327    //  //   //
1328    if ( mstrip >= 70 && mstrip < 75 ) lastrip = 3;  //   if ( mstrip >= 70 && mstrip < 75 ) lastrip = 3;
1329    //  //   //
1330    if ( mstrip >= 75 && mstrip < 84 ){      //   if ( mstrip >= 75 && mstrip < 84 ){    
1331      lastrip = (int)trunc((mstrip - 75)/3.) + 4;  //     lastrip = (int)trunc((mstrip - 75)/3.) + 4;
1332    };  //   };
1333    //  //   //
1334    if ( mstrip >= 84 && mstrip < 90 ){      //   if ( mstrip >= 84 && mstrip < 90 ){    
1335      lastrip = (int)trunc((mstrip - 84)/2.) + 7;  //     lastrip = (int)trunc((mstrip - 84)/2.) + 7;
1336    };  //   };
1337    //  //   //
1338    if ( mstrip >= 90 && mstrip < 101 ){      //   if ( mstrip >= 90 && mstrip < 101 ){    
1339      lastrip = mstrip - 90 + 10;  //     lastrip = mstrip - 90 + 10;
1340    //   };
1341    //   //
1342    //   if ( mstrip >= 101 && mstrip < 107 ){    
1343    //     lastrip = (int)trunc((mstrip - 101)/2.) + 21;
1344    //   };
1345    //   //
1346    //   //
1347    //   if ( mstrip >= 107 && mstrip < 116 ){    
1348    //     lastrip = (int)trunc((mstrip - 107)/3.) + 24;
1349    //   };
1350    //   //
1351    //   if ( mstrip >= 116 && mstrip < 121 ) lastrip = 27;
1352    //   //
1353    //   if ( mstrip >= 121 && mstrip < 127 ) lastrip = 28;
1354    //   //
1355    //   if ( mstrip >= 127 && mstrip < 141 ) lastrip = 29;
1356    //   //
1357    //   if ( mstrip >= 141 ) lastrip = 30;
1358    //   //
1359      //
1360      if ( mstrip < 83 ) lastrip = 0;
1361      //
1362      if ( mstrip >= 83 && mstrip < 90 ) lastrip = 1;
1363      //
1364      if ( mstrip >= 90 && mstrip < 93 ) lastrip = 2;
1365      //
1366      if ( mstrip >= 93 && mstrip < 98 ){    
1367        lastrip = mstrip - 93 + 3;
1368    };    };
1369    //    //
1370    if ( mstrip >= 101 && mstrip < 107 ){        if ( mstrip >= 98 && mstrip < 101 ) lastrip = 8;
     lastrip = (int)trunc((mstrip - 101)/2.) + 21;  
   };  
   //  
   //  
   if ( mstrip >= 107 && mstrip < 116 ){      
     lastrip = (int)trunc((mstrip - 107)/3.) + 24;  
   };  
   //  
   if ( mstrip >= 116 && mstrip < 121 ) lastrip = 27;  
   //  
   if ( mstrip >= 121 && mstrip < 127 ) lastrip = 28;  
1371    //    //
1372    if ( mstrip >= 127 && mstrip < 141 ) lastrip = 29;    if ( mstrip >= 101 && mstrip < 107 ) lastrip = 9;
1373    //    //
1374    if ( mstrip >= 141 ) lastrip = 30;    if ( mstrip >= 107 ) lastrip = 10;
1375    //    //
1376    return(lastrip);    return(lastrip);
1377  }  }
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.23