--- calo/flight/CaloFranzini/src/CaloFranzini.cpp	2007/12/18 09:55:07	1.4
+++ calo/flight/CaloFranzini/src/CaloFranzini.cpp	2008/01/03 10:02:28	1.5
@@ -31,7 +31,11 @@
   dolong = true;
   dofull = false;
   sntr = 0;
+  OBT = 0;
+  PKT = 0;
+  atime = 0;
   //
+  crig = false;
   sel = true;
   cont = false;
   N = 0;
@@ -45,11 +49,8 @@
 
 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, 43*sizeof(Float_t));
@@ -81,7 +82,9 @@
 }
 
 void CaloFranzini::SetNoWpreSampler(Int_t n){
-  if ( NC+n < 23 ){
+  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");
@@ -146,6 +149,7 @@
   this->Clear();
   //
   Float_t rig = L2->GetTrack(itr)->GetTrkTrack()->GetRigidity();
+  if ( crig ) rig = L2->GetCaloLevel2()->qtot/260.;
   //
   // Fill the estrip matrix
   //
@@ -163,7 +167,7 @@
     nplane = 1 - view + 2 * (plane - N);
     if ( nplane > (37-(2*N)) ) nplane--; 
     //
-    if ( plane == (18+N) ) mip = 0.;
+//    if ( plane == (18+N) ) mip = 0.;
     if ( nplane > -1 ) qplane[nplane] += mip;
     //
   };
@@ -209,16 +213,27 @@
     // 
     degfre = TMath::Min(dgf,NC);
     //
+    Float_t longzdiag = 0.;
+    Float_t longzout = 0.;
+    //
     if ( degfre > 0 ){
       for (Int_t i = 0; i < degfre; i++){
 	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 ){
@@ -226,7 +241,7 @@
     printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
   };  
   //
-  if ( debug ) this->Print();
+//  if ( debug ) this->Print();
   if ( debug ) printf(" exit \n");
 }
 
@@ -292,6 +307,12 @@
   file->WriteObject(&(*qpl),name.Data());
 }
 
+void CaloFranzini::WriteFullMean(TMatrixD *qpl, Int_t bin){
+  file->cd();
+  TString name = Form("fqplmeann%i",bin);
+  file->WriteObject(&(*qpl),name.Data());
+}
+
 void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){
   file->cd();
   TString name = Form("matrixn%i",bin);
@@ -299,6 +320,13 @@
   file->WriteObject(&mat,name.Data());
 }
 
+void CaloFranzini::WriteInvertedFullMatrix(TMatrixF 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);
@@ -306,6 +334,13 @@
   file->WriteObject(&(*mat),name.Data());
 }
 
+void CaloFranzini::WriteFullMatrix(TMatrixF *mat, Int_t bin){
+  file->cd();
+  TString name = Form("origfmatrixn%i",bin);
+  //  mat.Write(name.Data());
+  file->WriteObject(&(*mat),name.Data());
+}
+
 void CaloFranzini::CloseMatrixFile(){
   file->cd();
   file->Close();
@@ -332,9 +367,29 @@
     printf(" %s \n",matrixfile.Data());
     return(false);
   };
-  if ( !this->LoadMatrices() ){
-    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);
+    };
   };
   //
   //
@@ -367,13 +422,32 @@
     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);
-	};
+    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);
@@ -389,6 +463,16 @@
   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){
   //
   Int_t mv = 0;
@@ -450,13 +534,13 @@
   //
   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));
+//      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));
+//      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;
   };
@@ -478,8 +562,8 @@
   return(ave);
   //
 }
-
-Float_t CaloFranzini::GetHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
+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) ){
@@ -491,15 +575,16 @@
   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));
+//      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));
+//      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;
   };
@@ -507,9 +592,13 @@
   minrig = brigm->At(therigb);
   maxrig = brigm->At(therigb+1);
   //
-  Float_t minene = (*hmat[therigb])[iindex][jindex];
-  Float_t maxene = (*hmat[therigb+1])[iindex][jindex];
+  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);
@@ -518,6 +607,120 @@
   //
 }
 
+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;
+    };
+  };
+  //
+  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;
+    };
+  };
+  //
+  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 = (*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);
@@ -560,3 +763,46 @@
   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");
+  //
+};