--- calo/flight/CaloFranzini/src/CaloFranzini.cpp	2007/12/04 12:05:29	1.1
+++ calo/flight/CaloFranzini/src/CaloFranzini.cpp	2007/12/13 17:08:11	1.3
@@ -29,6 +29,14 @@
   debug = false;
   dolong = true;
   dofull = false;
+  sntr = 0;
+  //
+  sel = true;
+  cont = false;
+  N = 0;
+  NC = 43;
+  //
+  mask18b = -1;
   //
   Clear();
   //
@@ -39,10 +47,11 @@
   OBT = 0;
   PKT = 0;
   atime = 0;
-  longtzeta = 0.;
-  fulltzeta = 0.;
+  longtzeta = -100.;
+  fulltzeta = -100.;
+  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));
   //
 }
 
@@ -53,21 +62,188 @@
   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");
   //
 }
 
 void CaloFranzini::Delete(){
+  //
   if ( file ) file->Close();
+  //
   Clear();
+  //
+}
+
+void CaloFranzini::SetNoWpreSampler(Int_t n){
+  if ( NC+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);
+}
+
+void CaloFranzini::Process(Int_t itr){
+  //  
+  if ( !L2 ){
+    printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
+    printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
+    return;
+  };
+  //
+  if ( !nbin || !file || !brig ){
+    printf(" ERROR: it seems covariance matrix file has not been opened (CaloFranzini::Open()) \n");
+    printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
+    return;
+  };
+  //
+  Bool_t newentry = false;
+  //
+  if ( L2->IsORB() ){
+    if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || itr != sntr ){
+      newentry = true;
+      OBT = L2->GetOrbitalInfo()->OBT;
+      PKT = L2->GetOrbitalInfo()->pkt_num;
+      atime = L2->GetOrbitalInfo()->absTime;
+      sntr = itr;
+    };
+  } else {
+    newentry = true;
+  };
+  //
+  if ( !newentry ) return;
+  //
+  // Some variables
+  //
+  if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
+  //
+  this->Clear();
+  //
+  Float_t rig = L2->GetTrack(itr)->GetTrkTrack()->GetRigidity();
+  //
+  // Fill the estrip matrix
+  //
+  Int_t nplane = 0;
+  Int_t view = 0;
+  Int_t plane = 0;
+  Int_t strip = 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;
+    //
+  };
+  //
+  //
+  //
+  if ( dolong ){
+    //
+    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 ){
+	dgf = 2 * i;
+	break;
+      };
+      if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){
+      dgf = 1 + 2 * i;
+      break;
+      };
+    };
+    //
+    if ( dgf < 43 && dgf > 37 ) dgf--; 
+    // 
+    degfre = TMath::Min(dgf,NC);
+    //
+    if ( degfre > 0 ){
+      TArrayF *qplmean = this->LoadLongAverage(rig);
+      TMatrixD *lmatrix = this->LoadCovarianceMatrix(rig);
+      for (Int_t i = 0; i < degfre; i++){
+	if ( i != mask18b ){
+	  for (Int_t j = 0; j < degfre; j++){	
+	    if ( j != mask18b ){
+	      longtzeta += (qplane[i]-qplmean->At(i)) * (*lmatrix)[i][j] * (qplane[j]-qplmean->At(j));
+	    };
+	  };
+	};
+      };
+    };
+  };
+  if ( dofull ){
+    printf(" ERROR: NOT IMPLEMENTED YET\n");
     printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
   };  
   //
   if ( debug ) this->Print();
   if ( debug ) printf(" exit \n");
-  //
+}
+
+
+Float_t CaloFranzini::GetNormLongTZeta(){
+  Process(); 
+  Float_t normz = 0.;
+  if ( degfre > 0 ) normz = longtzeta/(Float_t)degfre;
+  return normz; 
+}
+
+Float_t CaloFranzini::GetNormFullTZeta(){
+  Process(); 
+  Float_t normz = 0.;
+  if ( degfre > 0 ) normz = fulltzeta/(Float_t)degfre;
+  return normz; 
 }
 
 Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){  
@@ -76,6 +252,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);
@@ -86,25 +263,39 @@
 }
 
 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::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::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::CloseMatrixFile(){
+  file->cd();
   file->Close();
 }
 
@@ -160,9 +351,9 @@
     name = "matrixn0";
     printf(" Using matrix %s \n",name.Data());
   };
-  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);
+  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));
+    name = Form("matrixn%i",nbin-2);
     printf(" Using matrix %s \n",name.Data());
   };
   //
@@ -187,9 +378,9 @@
     name = "qplmeann0";
     printf(" Using qplmean %s \n",name.Data());
   };
-  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);
+  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));
+    name = Form("qplmeann%i",nbin-2);
     printf(" Using qplmean %s \n",name.Data());
   };
   //
@@ -198,3 +389,48 @@
   return(qplmean);
   //
 }
+
+
+
+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();
+  //  Int_t pp=0;
+  Float_t qpl[43];
+  for (Int_t st=0;st<43;st++){
+    qpl[st] = ll->At(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");
+  //
+};