--- calo/flight/CaloFranzini/src/CaloFranzini.cpp	2007/12/04 12:42:54	1.2
+++ calo/flight/CaloFranzini/src/CaloFranzini.cpp	2008/01/21 10:24:11	1.7
@@ -18,6 +18,7 @@
   //
   file = NULL;
   brig = NULL;
+  brigm = NULL;
   nbin = 0;
   //
   L2 = l2p;
@@ -30,6 +31,17 @@
   dolong = true;
   dofull = false;
   sntr = 0;
+  OBT = 0;
+  PKT = 0;
+  atime = 0;
+  //
+  crig = false;
+  sel = true;
+  cont = false;
+  N = 0;
+  NC = 43;
+  //
+  mask18b = -1;
   //
   Clear();
   //
@@ -37,14 +49,11 @@
 
 void CaloFranzini::Clear(){
   //
-  OBT = 0;
-  PKT = 0;
-  atime = 0;
-  longtzeta = -100.;
-  fulltzeta = -100.;
+  longtzeta = 0.;
+  fulltzeta = 0.;
   degfre = 0;
   memset(estrip, 0, 2*22*96*sizeof(Float_t));
-  memset(qplane, 0, 2*22*sizeof(Float_t));
+  memset(qplane, 0, 43*sizeof(Float_t));
   //
 }
 
@@ -55,8 +64,11 @@
   printf("========================================================================\n");
   printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
   printf(" debug                [debug flag]:.. %i\n",debug);
+  printf(" degree of freedom                :.. %i\n",this->GetDegreeOfFreedom());  
   printf(" longtzeta                        :.. %f\n",longtzeta);  
-  printf(" fulltzeta                        :.. %f\n",fulltzeta);  
+  printf(" longtzeta normalized             :.. %f\n",this->GetNormLongTZeta());  
+  //  printf(" fulltzeta                        :.. %f\n",fulltzeta);  
+  //  printf(" longtzeta normalized             :.. %f\n",this->GetNormFullTZeta());  
   printf("========================================================================\n");
   //
 }
@@ -69,6 +81,32 @@
   //
 }
 
+void CaloFranzini::SetNoWpreSampler(Int_t n){
+  Int_t nc2 = NC/2;
+  if ( NC >= 37 ) nc2 = (NC+1)/2;
+  if ( nc2+n < 23 ){
+    N = n;
+  } else {
+    printf(" ERROR! Calorimeter is made of 22 W planes\n");
+    printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
+    printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
+    NC = 43;
+    N = 0;
+  };
+}
+
+void CaloFranzini::SetNoWcalo(Int_t n){
+  if ( N+n < 23 ){
+    NC = n*2;
+    if ( NC >37 ) NC--;
+  } else {
+    printf(" ERROR! Calorimeter is made of 22 W planes\n");
+    printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
+    printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
+    NC = 43;
+    N = 0;
+  };
+}
 
 void CaloFranzini::Process(){
   this->Process(0);
@@ -111,6 +149,7 @@
   this->Clear();
   //
   Float_t rig = L2->GetTrack(itr)->GetTrkTrack()->GetRigidity();
+  if ( crig ) rig = L2->GetCaloLevel2()->qtot/260.;
   //
   // Fill the estrip matrix
   //
@@ -124,8 +163,12 @@
     mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
     //
     estrip[view][plane][strip] = mip;
-    nplane = 1 - view + 2 * plane;
-    qplane[nplane] += 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;
     //
   };
   //
@@ -133,27 +176,64 @@
   //
   if ( dolong ){
     //
-    degfre = 44;
+    if ( cont ){
+      for (Int_t i=0; i<22; i++){
+ 	if ( i == (18+N) ){
+	  mask18b =  1 + 2 * (i - N);
+	  break;
+	};
+      };
+    };  
+    //  
+    if ( sel ){
+      for (Int_t i=0; i<22; i++){
+ 	if ( i == (18-N) ){
+	  mask18b =  1 + 2 * (i - N);
+	  break;
+	};
+      };
+    };
+    //
+    if ( mask18b == 37 ) mask18b = -1;
+    //
+    Int_t dgf = 43;
     //
     for (Int_t i=0; i < 22; i++){
       if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){
-	degfre = 2 * i;
+	dgf = 2 * i;
 	break;
       };
       if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){
-      degfre = 1 + 2 * i;
+      dgf = 1 + 2 * i;
       break;
       };
     };
     //
+    if ( dgf < 43 && dgf > 37 ) dgf--; 
+    // 
+    degfre = TMath::Min(dgf,NC);
+    //
+    Float_t longzdiag = 0.;
+    Float_t longzout = 0.;
+    //
     if ( degfre > 0 ){
-      TArrayF *qplmean = this->LoadLongAverage(rig);
-      TMatrixD *lmatrix = this->LoadCovarianceMatrix(rig);
       for (Int_t i = 0; i < degfre; i++){
-	for (Int_t j = 0; j < degfre; j++){	
-	  longtzeta += (qplane[i]-qplmean->At(i)) * (*lmatrix)[i][j] * (qplane[j]-qplmean->At(j));
+	if ( i != mask18b ){
+	  for (Int_t j = 0; j < degfre; j++){	
+	    if ( j != mask18b ){
+	      if ( i == j ){
+		longzdiag += (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));
+	      } else {
+		longzout += (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));
+	      };
+	      longtzeta += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
+	    };
+	  };
 	};
       };
+      if ( debug ) printf(" diagonal total %f out of diagonal total %f \n",longzdiag,longzout);
     };
   };
   if ( dofull ){
@@ -161,7 +241,7 @@
     printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
   };  
   //
-  if ( debug ) this->Print();
+//  if ( debug ) this->Print();
   if ( debug ) printf(" exit \n");
 }
 
@@ -186,6 +266,7 @@
   //
   if ( !file || file->IsZombie() ){
     file = new TFile(matrixfile.Data(),"RECREATE"); 
+    printf(" Create file %s \n",file->GetName());
   } else {
     printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
     return(false);
@@ -195,26 +276,90 @@
   //
 }
 
+Bool_t CaloFranzini::UpdateMatrixFile(TString matrixfile){  
+  //
+  file = new TFile(matrixfile.Data(),"UPDATE");
+  //
+  if ( !file || file->IsZombie() ){
+    printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
+    return(false);
+  };
+  //
+  return(true);
+  //
+}
+
 void CaloFranzini::WriteNumBin(Int_t numbin){
+  file->cd();
   TArrayI nbi(1, &numbin);
   file->WriteObject(&nbi, "nbinenergy");
 }
 
 void CaloFranzini::WriteRigBin(TArrayF *rigbin){
-  file->WriteObject(&rigbin, "binrig");
+  file->cd();
+  //  rigbin->Write("binrig");
+  file->WriteObject(&(*rigbin), "binrig");
 }
 
 void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){
+  file->cd();
   TString name = Form("qplmeann%i",bin);
-  file->WriteObject(&qpl,name.Data());
+  file->WriteObject(&(*qpl),name.Data());
 }
 
-void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
+void CaloFranzini::WriteFullMean(TMatrixD *qpl, Int_t bin){
+  file->cd();
+  TString name = Form("fqplmeann%i",bin);
+  file->WriteObject(&(*qpl),name.Data());
+  file->Purge();
+}
+
+void CaloFranzini::WriteFullNMean(TMatrixD *qpl, Int_t bin){
+  file->cd();
+  TString name = Form("fnqplmeann%i",bin);
+  file->WriteObject(&(*qpl),name.Data());
+  file->Purge();
+}
+
+void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){
+  file->cd();
   TString name = Form("matrixn%i",bin);
+  //  mat.Write(name.Data());
+  file->WriteObject(&mat,name.Data());
+}
+
+void CaloFranzini::WriteInvertedFullMatrix(TMatrixD mat, Int_t bin){
+  file->cd();
+  TString name = Form("fmatrixn%i",bin);
+  //  mat.Write(name.Data());
+  file->WriteObject(&mat,name.Data());
+}
+
+void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
+  file->cd();
+  TString name = Form("origmatrixn%i",bin);
+  //  mat.Write(name.Data());
+  file->WriteObject(&(*mat),name.Data());
+}
+
+void CaloFranzini::WriteFullMatrix(TMatrixD *mat, Int_t bin){
+  file->cd();
+  TString name = Form("origfmatrixn%i",bin);
+  //  mat.Write(name.Data());
+  file->WriteObject(&(*mat),name.Data());
+  file->Purge();
+}
+
+void CaloFranzini::WriteFullNMatrix(TMatrixF *mat, Int_t bin){
+  file->cd();
+  TString name = Form("origfnmatrixn%i",bin);
+  //  mat.Write(name.Data());
   file->WriteObject(&(*mat),name.Data());
+  file->Purge();
 }
 
 void CaloFranzini::CloseMatrixFile(){
+  file->cd();
   file->Close();
 }
 
@@ -235,76 +380,590 @@
     return(false);
   };
   //
+  if ( !this->LoadBin() ){
+    printf(" %s \n",matrixfile.Data());
+    return(false);
+  };
+  //
+  if ( dolong ){
+    if ( !this->LoadLong() ){
+      printf(" %s \n",matrixfile.Data());
+      return(false);
+    };
+    //
+    if ( !this->LoadMatrices() ){
+      printf(" %s \n",matrixfile.Data());
+      return(false);
+    };
+  };
+  //
+  if ( dofull ){
+    if ( !this->LoadFull() ){
+      printf(" %s \n",matrixfile.Data());
+      return(false);
+    };
+    //
+    if ( !this->LoadFullMatrices() ){
+      printf(" %s \n",matrixfile.Data());
+      return(false);
+    };
+  };
+  //
+  //
+  return(true);
+  //
+}
+
+Bool_t CaloFranzini::LoadBin(){
+  //
   TArrayI *numbin = (TArrayI*)file->Get("nbinenergy");
   if ( !numbin ){
-    printf(" ERROR: cannot read number of bins from file %s \n",matrixfile.Data());
+    printf(" ERROR: cannot read number of bins from file ");
     return(false);
   };
   nbin = (Int_t)numbin->At(0);
   if ( nbin <= 0 ){
-    printf(" ERROR: cannot work with 0 energy bins (from file %s) \n",matrixfile.Data());
+    printf(" ERROR: cannot work with 0 energy bins from file ");
     return(false);
   };
   //
   brig = (TArrayF*)file->Get("binrig");
   if ( !brig ){
-    printf(" ERROR: cannot read rigidity binning from file %s \n",matrixfile.Data());
+    printf(" ERROR: cannot read rigidity binning from file ");
     return(false);
   };
   //
+  brigm=(TArrayF*)file->Get("binrigmean");
+  if ( !brigm ){
+    printf(" ERROR: cannot read mean rigidity binning from file ");
+    return(false);
+  };
+  //
+  return(true);
+}
+
+Bool_t CaloFranzini::LoadLong(){
+  //
+  for (Int_t i=0;i<17;i++){
+    TString name = Form("qplmeann%i",i);
+    qplmean[i] = (TArrayF*)file->Get(name.Data());
+    if ( !qplmean[i] ){
+      printf(" ERROR: cannot read average from file ");
+      return(false);
+    };
+  };
+  //
+  return(true);
+}
+
+Bool_t CaloFranzini::LoadFull(){
+  //
+  for (Int_t i=0;i<17;i++){
+    TString name = Form("fqplmeann%i",i);
+    fqplmean[i] = (TMatrixD*)file->Get(name.Data());
+    if ( !fqplmean[i] ){
+      printf(" ERROR: cannot read average from file ");
+      return(false);
+    };
+  };
+  //
   return(true);
+}
+
+Bool_t CaloFranzini::LoadMatrices(){
+  //
+  for (Int_t i=0;i<17;i++){
+        TString name1 = Form("matrixn%i",i);
+        hmat[i] = (TMatrixD*)file->Get(name1.Data());
+  };
   //
+  return(true);
+}
+
+Bool_t CaloFranzini::LoadFullMatrices(){
+  //
+  for (Int_t i=0;i<17;i++){
+        TString name1 = Form("fmatrixn%i",i);
+        hfmat[i] = (TMatrixF*)file->Get(name1.Data());
+  };
+  //
+  return(true);
 }
 
 TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){
   //
-  TString name;
+  Int_t mv = 0;
   for (Int_t i = 0; i<nbin-1; i++){
     if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
-      name = Form("matrixn%i",i);
+      mv = i;
       break;
     };
   };
   if ( 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));
-    name = "matrixn0";
-    printf(" Using matrix %s \n",name.Data());
+    mv = 0;
   };
-  if ( 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));
-    name = Form("matrixn%i",nbin-1);
-    printf(" Using matrix %s \n",name.Data());
+  if ( rig >= brig->At(nbin-1) ){
+    printf(" WARNING: Event with rigidity higher than the last covariance matrix bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
+    mv = nbin-2;
   };
   //
-  TMatrixD *matrix = (TMatrixD*)file->Get(name.Data());
+  return(hmat[mv]);
+  //
+}
+
+TMatrixD *CaloFranzini::LoadFullAverage(Int_t rigbin){
+  //
+  TString name = Form("fqplmeann%i",rigbin);
+  TMatrixD *fmean=(TMatrixD*)file->Get(name.Data());
+  //
+  return(fmean);
+  //
+}
+
+TMatrixF *CaloFranzini::LoadFullMatrix(Int_t rigbin){
+  //
+  TString name = Form("origfmatrixn%i",rigbin);
+  TMatrixF *fmatri=(TMatrixF*)file->Get(name.Data());
+  //
+  return(fmatri);
+  //
+}
+
+void CaloFranzini::LoadFullMatrix(Int_t rigbin, TMatrixF *&fmatri){
+  //
+  TString name = Form("origfmatrixn%i",rigbin);
+  fmatri=(TMatrixF*)file->Get(name.Data());
+  //
+}
+
+void CaloFranzini::UnLoadFullAverage(Int_t rigbin){
+  //
+  TString name = Form("fqplmeann%i",rigbin);
+  file->Delete(name.Data());
+  //
+}
+
+void CaloFranzini::UnLoadFullMatrix(Int_t rigbin){
+  //
+  TString name = Form("origfmatrixn%i",rigbin);
+  file->Delete(name.Data());
+  //
+}
+
+TMatrixF *CaloFranzini::LoadFullNMatrix(Int_t rigbin){
+  //
+  TString name = Form("origfnmatrixn%i",rigbin);
+  TMatrixF *fnmatri=(TMatrixF*)file->Get(name.Data());
+  //
+  return(fnmatri);
+  //
+}
+
+void CaloFranzini::UnLoadFullNMatrix(Int_t rigbin){
+  //
+  TString name = Form("origfnmatrixn%i",rigbin);
+  file->Delete(name.Data());
+  //
+}
+
+TMatrixD *CaloFranzini::LoadFullNAverage(Int_t rigbin){
+  //
+  TString name = Form("fnqplmeann%i",rigbin);
+  TMatrixD *fnmean=(TMatrixD*)file->Get(name.Data());
+  //
+  return(fnmean);
+  //
+}
+
+void CaloFranzini::UnLoadFullNAverage(Int_t rigbin){
   //
-  return(matrix);
+  TString name = Form("fnqplmeann%i",rigbin);
+  file->Delete(name.Data());
   //
 }
 
 
 TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){
   //
-  TString name;
+  Int_t mv=0;
   for (Int_t i = 0; i<nbin-1; i++){
     if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
-      name = Form("qplmeann%i",i);
+	mv = i;
       break;
     };
   };
   if ( 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));
-    name = "qplmeann0";
-    printf(" Using qplmean %s \n",name.Data());
+    mv = 0;
+  };
+  if ( rig >= brig->At(nbin-1) ){
+    printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
+    mv=nbin-2;
+  };
+  //
+  return(qplmean[mv]);
+  //
+}
+
+Float_t CaloFranzini::GetAverageAt(Int_t plane, Float_t rig){
+  //
+  Int_t therigb = 0;
+  for (Int_t i = 0; i<nbin-2; i++){
+    if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
+      therigb = i;
+      break;
+    };
+  };
+  //
+  Float_t minrig;
+  Float_t maxrig;
+  //
+  //
+  if ( rig < brigm->At(0) ){
+    if ( rig < brig->At(0) ){
+//      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
+    };
+    therigb = 0;
+  };
+  if ( rig >= brigm->At(nbin-2) ){
+    if ( rig >= brig->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));
+    };
+    therigb = nbin - 3;
+  };
+  //
+  minrig = brigm->At(therigb);
+  maxrig = brigm->At(therigb+1);
+  //
+  Float_t minene = (*qplmean[therigb])[plane];
+  Float_t maxene = (*qplmean[therigb+1])[plane];
+  //
+  if ( maxrig == minrig ){
+   printf("Unrecoverable ERROR! Matrix will be screwed up... \n");
+   return(0.);
+  };
+  Float_t b = log(maxene/minene)/(maxrig-minrig);
+  Float_t a = minene/exp(minrig*b);
+  Float_t ave = a*exp(b*rig);
+  if ( a == 0. || minene == 0. || ave != ave ){
+    //  if ( a == 0. || minene == 0. ){
+    Float_t m = (maxene-minene)/(maxrig-minrig);
+    Float_t q = minene - m * minrig;
+    ave = rig * m + q;
+  };
+  //
+  return(ave);
+  //
+}
+
+Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig){
+  //
+  Int_t therigb = 0;
+  for (Int_t i = 0; i<nbin-2; i++){
+    if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
+      therigb = i;
+      break;
+    };
   };
-  if ( rig >= brig->At(nbin) ){
-    printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin));
-    name = Form("qplmeann%i",nbin-1);
-    printf(" Using qplmean %s \n",name.Data());
+  //
+  //
+  if ( rig < brigm->At(0) ){
+    if ( rig < brig->At(0) ){
+//      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
+    };
+    therigb = 0;
+  };
+  if ( rig >= brigm->At(nbin-2) ){
+    if ( rig >= brig->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));
+    };
+    therigb = nbin - 3;
+  };
+  //
+  return(this->GetFullAverageAt(plane,strip,rig,therigb));
+  //
+}
+
+Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig, Int_t therigb){
+  //
+  Float_t minrig;
+  Float_t maxrig;
+  //
+  minrig = brigm->At(therigb);
+  maxrig = brigm->At(therigb+1);
+  //
+  Float_t minene = (*fqplmean[therigb])[plane][strip];
+  Float_t maxene = (*fqplmean[therigb+1])[plane][strip];
+  //
+  if ( maxrig == minrig ){
+   printf("Unrecoverable ERROR! Matrix will be screwed up... \n");
+   return(0.);
+  };
+  Float_t b = log(maxene/minene)/(maxrig-minrig);
+  Float_t a = minene/exp(minrig*b);
+  Float_t ave = a*exp(b*rig);
+  if ( a == 0. || minene == 0. || ave != ave ){
+    Float_t m = (maxene-minene)/(maxrig-minrig);
+    Float_t q = minene - m * minrig;
+    ave = rig * m + q;
+  };
+  //  
+  //  ave += (44.-plane)*strip;
+  //if ( a == 0. ) ave = 0.;
+  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);
+  //
+  return(ave);
+  //
+}
+
+Float_t CaloFranzini::GetHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
+  Int_t therigb = 0;
+  for (Int_t i = 0; i<nbin-2; i++){
+    if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
+      therigb = i;
+      break;
+    };
   };
   //
-  TArrayF *qplmean = (TArrayF*)file->Get(name.Data());
+  Float_t minrig;
+  Float_t maxrig;
+  //
+  if ( rig < brigm->At(0) ){
+    if ( rig < brig->At(0) ){
+//      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
+    };
+    therigb = 0;
+  };
+//  if ( rig >= brigm->At(nbin-4) ){ // -2
+  if ( rig >= brigm->At(nbin-2) ){ 
+    if ( rig >= brig->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));
+    };
+//    therigb = nbin - 5;// -3
+    therigb = nbin - 3;
+  };
+  //
+  if ( therigb < 2 ) therigb = 2;
+  minrig = brigm->At(therigb);
+  maxrig = brigm->At(therigb+1);
+  //  printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);
+  //
+  Float_t minene = (*hmat[therigb])[iindex][jindex];
+  Float_t maxene = (*hmat[therigb+1])[iindex][jindex];
+  //  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);
+  //
+//  Float_t a = 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);
+//  };
+  //
+  Float_t m = (maxene-minene)/(maxrig-minrig);
+  Float_t q = minene - m * minrig;
+  Float_t ave = rig * m + q;
+  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);
+  //
+  //
+  return(ave);
+  //
+}
+
+Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
+  Int_t therigb = 0;
+  for (Int_t i = 0; i<nbin-2; i++){
+    if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
+      therigb = i;
+      break;
+    };
+  };
+  //
+  if ( rig < brigm->At(0) ){
+    if ( rig < brig->At(0) ){
+//      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
+    };
+    therigb = 0;
+  };
+//  if ( rig >= brigm->At(nbin-4) ){ // -2
+  if ( rig >= brigm->At(nbin-2) ){ 
+    if ( rig >= brig->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));
+    };
+//    therigb = nbin - 5;// -3
+    therigb = nbin - 3;
+  };
+  //
+  if ( therigb < 2 ) therigb = 2;
+  //
+  return(this->GetFullHmatrixAt(iindex,jindex,rig,therigb));
+  //
+}
+
+Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig, Int_t therigb){
+  //
+  Float_t minrig;
+  Float_t maxrig;
+  //
+  minrig = brigm->At(therigb);
+  maxrig = brigm->At(therigb+1);
+  //  printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);
+  //
+  Float_t minene = (*hfmat[therigb])[iindex][jindex];
+  Float_t maxene = (*hfmat[therigb+1])[iindex][jindex];
+  //  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);
+  //
+//  Float_t a = 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);
+//  };
+  //
+  Float_t m = (maxene-minene)/(maxrig-minrig);
+  Float_t q = minene - m * minrig;
+  Float_t ave = rig * m + q;
+  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);
+  //
+  //
+  return(ave);
+  //
+}
+
+void CaloFranzini::DrawLongAverage(Float_t rig){
+  //
+  TArrayF *ll = this->LoadLongAverage(rig);
+  //
+  gStyle->SetLabelSize(0.04);
+  gStyle->SetNdivisions(510,"XY");
+  //
+  TString hid = Form("cfralongvyvx");	
+  TCanvas *tcf  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
+  if ( tcf ){
+    tcf->cd();
+  } else {
+    tcf = new TCanvas(hid,hid);
+  };
+  //
+  TString thid = Form("hfralongvyvx");	
+  TH1F *thf  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
+  if ( thf ) thf->Delete();
+  thf = new TH1F(thid,thid,44,-0.5,43.5);
+  tcf->cd();
+  Float_t qpl[43];
+  for (Int_t st=0;st<43;st++){
+    qpl[st] = ll->At(st);
+    printf("st %i qpl %f\n",st,qpl[st]);
+  };
+  for (Int_t st=0;st<44;st++){
+    if ( st == 37 ){
+      thf->Fill(st,0.);
+    } else {
+      Int_t ss = st;
+      if ( st > 37 ) ss--;
+      thf->Fill(st,qpl[ss]);
+    };
+  };
+  thf->Draw();
+  tcf->Modified();
+  tcf->Update();
+  //
+  gStyle->SetLabelSize(0);
+  gStyle->SetNdivisions(1,"XY");
+  //
+};
+
+void CaloFranzini::DrawLongAverage(Int_t bin){
+  //
+  TArrayF *ll = this->LoadLongAverage(brigm->At(bin));
+  //
+  gStyle->SetLabelSize(0.04);
+  gStyle->SetNdivisions(510,"XY");
+  //
+  TString hid = Form("cfralongvyvx");	
+  TCanvas *tcf  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
+  if ( tcf ){
+    tcf->cd();
+  } else {
+    tcf = new TCanvas(hid,hid);
+  };
+  //
+  TString thid = Form("hfralongvyvx");	
+  TH1F *thf  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
+  if ( thf ) thf->Delete();
+  thf = new TH1F(thid,thid,44,-0.5,43.5);
+  tcf->cd();
+  Float_t qpl[43];
+  for (Int_t st=0;st<43;st++){
+    qpl[st] = ll->At(st);
+    printf("st %i qpl %f\n",st,qpl[st]);
+  };
+  for (Int_t st=0;st<44;st++){
+    if ( st == 37 ){
+      thf->Fill(st,0.);
+    } else {
+      Int_t ss = st;
+      if ( st > 37 ) ss--;
+      thf->Fill(st,qpl[ss]);
+    };
+  };
+  thf->Draw();
+  tcf->Modified();
+  tcf->Update();
+  //
+  gStyle->SetLabelSize(0);
+  gStyle->SetNdivisions(1,"XY");
+  //
+};
+
+Int_t CaloFranzini::ConvertStrip(Int_t mstrip){  
+  //
+  Int_t lastrip = 0;
+  //
+  if ( mstrip < 50 ) lastrip = 0;
+  //
+  if ( mstrip >= 50 && mstrip < 64 ) lastrip = 1;
+  //
+  if ( mstrip >= 64 && mstrip < 70 ) lastrip = 2;
+  //
+  if ( mstrip >= 70 && mstrip < 75 ) lastrip = 3;
+  //
+  if ( mstrip >= 75 && mstrip < 84 ){    
+    lastrip = (int)trunc((mstrip - 75)/3.) + 4;
+  };
+  //
+  if ( mstrip >= 84 && mstrip < 90 ){    
+    lastrip = (int)trunc((mstrip - 84)/2.) + 7;
+  };
+  //
+  if ( mstrip >= 90 && mstrip < 101 ){    
+    lastrip = mstrip - 90 + 10;
+  };
+  //
+  if ( mstrip >= 101 && mstrip < 107 ){    
+    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;
+  //
+  if ( mstrip >= 127 && mstrip < 141 ) lastrip = 29;
   //
-  return(qplmean);
+  if ( mstrip >= 141 ) lastrip = 30;
   //
+  return(lastrip);
 }