CaloFranzini/src/CaloFranzini.cpp

/**
 * \file CaloFranzini.cpp
 * \author Emiliano Mocchiutti (2007/07/18)
 */
//
// headers
//
#include <CaloFranzini.h>
//--------------------------------------
/**
 * Default constructor 
 */
CaloFranzini::CaloFranzini(){
  Clear();
}

CaloFranzini::CaloFranzini(PamLevel2 *l2p){  
  //
  file = NULL;
  brig = NULL;
  brigm = NULL;
  nbin = 0;
  //
  L2 = l2p;
  //
  if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
  //
  // Default variables
  //
  debug = false;
  dolong = true;
  dofull = false;
  sntr = 0;
  //
  sel = true;
  cont = false;
  N = 0;
  NC = 43;
  //
  mask18b = -1;
  //
  Clear();
  //
}

void CaloFranzini::Clear(){
  //
  OBT = 0;
  PKT = 0;
  atime = 0;
  longtzeta = -100.;
  fulltzeta = -100.;
  degfre = 0;
  memset(estrip, 0, 2*22*96*sizeof(Float_t));
  memset(qplane, 0, 43*sizeof(Float_t));
  //
}

void CaloFranzini::Print(){
  //
  Process();
  //
  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(" 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 ){
      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]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
            };
          };
        };
      };
    };
  };
  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){  
  //
  file = new TFile(matrixfile.Data(),"READ");
  //
  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);
  };
  //
  return(true);
  //
}

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->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());
}

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();
}


Bool_t CaloFranzini::Open(TString matrixfile){  
  //
  // find matrix file
  //
  if ( !strcmp(matrixfile.Data(),"") ){
    if (dolong) matrixfile = (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);
  };
  //
  if ( !this->LoadBin() ){
    printf(" %s \n",matrixfile.Data());
    return(false);
  };
  if ( !this->LoadMatrices() ){
    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 ");
    return(false);
  };
  nbin = (Int_t)numbin->At(0);
  if ( nbin <= 0 ){
    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 ");
    return(false);
  };
  //
  brigm=(TArrayF*)file->Get("binrigmean");
  if ( !brigm ){
    printf(" ERROR: cannot read mean rigidity binning from file ");
    return(false);
  };
  //
  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::LoadMatrices(){
  //
  for (Int_t i=0;i<17;i++){
        TString name1 = Form("matrixn%i",i);
        hmat[i] = (TMatrixD*)file->Get(name1.Data());
  };
  //
  return(true);
}

TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){
  //
  Int_t mv = 0;
  for (Int_t i = 0; i<nbin-1; i++){
    if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
      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));
    mv = 0;
  };
  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;
  };
  //
  return(hmat[mv]);
  //
}


TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){
  //
  Int_t mv=0;
  for (Int_t i = 0; i<nbin-1; i++){
    if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
        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));
    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);
  //
  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;
    };
  };
  //
  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 = (*hmat[therigb])[iindex][jindex];
  Float_t maxene = (*hmat[therigb+1])[iindex][jindex];
  //
  Float_t b = log(maxene/minene)/(maxrig-minrig);
  Float_t a = minene/exp(minrig*b);
  Float_t ave = a*exp(b*rig);
  //
  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");
  //
};
1	mocchiut	1.1	/**
2			* \file CaloFranzini.cpp
3			* \author Emiliano Mocchiutti (2007/07/18)
4			*/
5			//
6			// headers
7			//
8			#include <CaloFranzini.h>
9			//--------------------------------------
10			/**
11			* Default constructor
12			*/
13			CaloFranzini::CaloFranzini(){
14			Clear();
15			}
16
17			CaloFranzini::CaloFranzini(PamLevel2 *l2p){
18			//
19			file = NULL;
20			brig = NULL;
21	mocchiut	1.4	brigm = NULL;
22	mocchiut	1.1	nbin = 0;
23			//
24			L2 = l2p;
25			//
26			if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
27			//
28			// Default variables
29			//
30			debug = false;
31			dolong = true;
32			dofull = false;
33	mocchiut	1.2	sntr = 0;
34	mocchiut	1.1	//
35	mocchiut	1.3	sel = true;
36			cont = false;
37			N = 0;
38			NC = 43;
39			//
40			mask18b = -1;
41			//
42	mocchiut	1.1	Clear();
43			//
44			}
45
46			void CaloFranzini::Clear(){
47			//
48			OBT = 0;
49			PKT = 0;
50			atime = 0;
51	mocchiut	1.2	longtzeta = -100.;
52			fulltzeta = -100.;
53			degfre = 0;
54	mocchiut	1.1	memset(estrip, 0, 22296*sizeof(Float_t));
55	mocchiut	1.3	memset(qplane, 0, 43*sizeof(Float_t));
56	mocchiut	1.1	//
57			}
58
59			void CaloFranzini::Print(){
60			//
61			Process();
62			//
63			printf("========================================================================\n");
64			printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
65			printf(" debug [debug flag]:.. %i\n",debug);
66	mocchiut	1.3	printf(" degree of freedom :.. %i\n",this->GetDegreeOfFreedom());
67	mocchiut	1.1	printf(" longtzeta :.. %f\n",longtzeta);
68	mocchiut	1.3	printf(" longtzeta normalized :.. %f\n",this->GetNormLongTZeta());
69			// printf(" fulltzeta :.. %f\n",fulltzeta);
70			// printf(" longtzeta normalized :.. %f\n",this->GetNormFullTZeta());
71	mocchiut	1.1	printf("========================================================================\n");
72			//
73			}
74
75			void CaloFranzini::Delete(){
76	mocchiut	1.2	//
77	mocchiut	1.1	if ( file ) file->Close();
78	mocchiut	1.2	//
79	mocchiut	1.1	Clear();
80	mocchiut	1.2	//
81			}
82
83	mocchiut	1.3	void CaloFranzini::SetNoWpreSampler(Int_t n){
84			if ( NC+n < 23 ){
85			N = n;
86			} else {
87			printf(" ERROR! Calorimeter is made of 22 W planes\n");
88			printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
89			printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
90			NC = 43;
91			N = 0;
92			};
93			}
94
95			void CaloFranzini::SetNoWcalo(Int_t n){
96			if ( N+n < 23 ){
97			NC = n*2;
98			if ( NC >37 ) NC--;
99			} else {
100			printf(" ERROR! Calorimeter is made of 22 W planes\n");
101			printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
102			printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
103			NC = 43;
104			N = 0;
105			};
106			}
107	mocchiut	1.2
108			void CaloFranzini::Process(){
109			this->Process(0);
110			}
111
112			void CaloFranzini::Process(Int_t itr){
113			//
114			if ( !L2 ){
115			printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
116			printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
117			return;
118			};
119			//
120			if ( !nbin \|\| !file \|\| !brig ){
121			printf(" ERROR: it seems covariance matrix file has not been opened (CaloFranzini::Open()) \n");
122			printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
123			return;
124			};
125			//
126			Bool_t newentry = false;
127			//
128			if ( L2->IsORB() ){
129			if ( L2->GetOrbitalInfo()->pkt_num != PKT \|\| L2->GetOrbitalInfo()->OBT != OBT \|\| L2->GetOrbitalInfo()->absTime != atime \|\| itr != sntr ){
130			newentry = true;
131			OBT = L2->GetOrbitalInfo()->OBT;
132			PKT = L2->GetOrbitalInfo()->pkt_num;
133			atime = L2->GetOrbitalInfo()->absTime;
134			sntr = itr;
135			};
136			} else {
137			newentry = true;
138			};
139			//
140			if ( !newentry ) return;
141			//
142			// Some variables
143			//
144			if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
145			//
146			this->Clear();
147			//
148			Float_t rig = L2->GetTrack(itr)->GetTrkTrack()->GetRigidity();
149			//
150			// Fill the estrip matrix
151			//
152			Int_t nplane = 0;
153			Int_t view = 0;
154			Int_t plane = 0;
155			Int_t strip = 0;
156			Float_t mip = 0.;
157			for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
158			//
159			mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
160			//
161			estrip[view][plane][strip] = mip;
162	mocchiut	1.3	//
163			nplane = 1 - view + 2 * (plane - N);
164			if ( nplane > (37-(2*N)) ) nplane--;
165			//
166			if ( plane == (18+N) ) mip = 0.;
167			if ( nplane > -1 ) qplane[nplane] += mip;
168	mocchiut	1.2	//
169			};
170			//
171			//
172			//
173			if ( dolong ){
174			//
175	mocchiut	1.3	if ( cont ){
176			for (Int_t i=0; i<22; i++){
177			if ( i == (18+N) ){
178			mask18b = 1 + 2 * (i - N);
179			break;
180			};
181			};
182			};
183			//
184			if ( sel ){
185			for (Int_t i=0; i<22; i++){
186			if ( i == (18-N) ){
187			mask18b = 1 + 2 * (i - N);
188			break;
189			};
190			};
191			};
192			//
193			if ( mask18b == 37 ) mask18b = -1;
194			//
195			Int_t dgf = 43;
196	mocchiut	1.2	//
197			for (Int_t i=0; i < 22; i++){
198			if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){
199	mocchiut	1.3	dgf = 2 * i;
200	mocchiut	1.2	break;
201			};
202			if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){
203	mocchiut	1.3	dgf = 1 + 2 * i;
204	mocchiut	1.2	break;
205			};
206			};
207			//
208	mocchiut	1.3	if ( dgf < 43 && dgf > 37 ) dgf--;
209			//
210			degfre = TMath::Min(dgf,NC);
211			//
212	mocchiut	1.2	if ( degfre > 0 ){
213			for (Int_t i = 0; i < degfre; i++){
214	mocchiut	1.3	if ( i != mask18b ){
215			for (Int_t j = 0; j < degfre; j++){
216			if ( j != mask18b ){
217	mocchiut	1.4	longtzeta += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
218	mocchiut	1.3	};
219			};
220	mocchiut	1.2	};
221			};
222			};
223			};
224			if ( dofull ){
225			printf(" ERROR: NOT IMPLEMENTED YET\n");
226	mocchiut	1.1	printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
227			};
228			//
229			if ( debug ) this->Print();
230			if ( debug ) printf(" exit \n");
231	mocchiut	1.2	}
232
233
234			Float_t CaloFranzini::GetNormLongTZeta(){
235			Process();
236			Float_t normz = 0.;
237			if ( degfre > 0 ) normz = longtzeta/(Float_t)degfre;
238			return normz;
239			}
240
241			Float_t CaloFranzini::GetNormFullTZeta(){
242			Process();
243			Float_t normz = 0.;
244			if ( degfre > 0 ) normz = fulltzeta/(Float_t)degfre;
245			return normz;
246	mocchiut	1.1	}
247
248			Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){
249			//
250			file = new TFile(matrixfile.Data(),"READ");
251			//
252			if ( !file \|\| file->IsZombie() ){
253			file = new TFile(matrixfile.Data(),"RECREATE");
254	mocchiut	1.3	printf(" Create file %s \n",file->GetName());
255	mocchiut	1.1	} else {
256			printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
257			return(false);
258			};
259			//
260			return(true);
261			//
262			}
263
264	mocchiut	1.4	Bool_t CaloFranzini::UpdateMatrixFile(TString matrixfile){
265			//
266			file = new TFile(matrixfile.Data(),"UPDATE");
267			//
268			if ( !file \|\| file->IsZombie() ){
269			printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
270			return(false);
271			};
272			//
273			return(true);
274			//
275			}
276
277	mocchiut	1.1	void CaloFranzini::WriteNumBin(Int_t numbin){
278	mocchiut	1.3	file->cd();
279	mocchiut	1.1	TArrayI nbi(1, &numbin);
280			file->WriteObject(&nbi, "nbinenergy");
281			}
282
283			void CaloFranzini::WriteRigBin(TArrayF *rigbin){
284	mocchiut	1.3	file->cd();
285			// rigbin->Write("binrig");
286			file->WriteObject(&(*rigbin), "binrig");
287	mocchiut	1.1	}
288
289			void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){
290	mocchiut	1.3	file->cd();
291	mocchiut	1.1	TString name = Form("qplmeann%i",bin);
292	mocchiut	1.3	file->WriteObject(&(*qpl),name.Data());
293			}
294
295			void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){
296			file->cd();
297			TString name = Form("matrixn%i",bin);
298			// mat.Write(name.Data());
299			file->WriteObject(&mat,name.Data());
300	mocchiut	1.1	}
301
302			void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
303	mocchiut	1.3	file->cd();
304			TString name = Form("origmatrixn%i",bin);
305			// mat.Write(name.Data());
306	mocchiut	1.1	file->WriteObject(&(*mat),name.Data());
307			}
308
309			void CaloFranzini::CloseMatrixFile(){
310	mocchiut	1.3	file->cd();
311	mocchiut	1.1	file->Close();
312			}
313
314
315			Bool_t CaloFranzini::Open(TString matrixfile){
316			//
317			// find matrix file
318			//
319			if ( !strcmp(matrixfile.Data(),"") ){
320			if (dolong) matrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_longel.root";
321			if (dofull) matrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_fullel.root";
322			};
323			//
324			file = new TFile(matrixfile.Data(),"READ");
325			//
326			if ( !file \|\| file->IsZombie() ){
327			printf(" ERROR: cannot open file %s \n",matrixfile.Data());
328			return(false);
329			};
330			//
331	mocchiut	1.4	if ( !this->LoadBin() ){
332			printf(" %s \n",matrixfile.Data());
333			return(false);
334			};
335			if ( !this->LoadMatrices() ){
336			printf(" %s \n",matrixfile.Data());
337			return(false);
338			};
339			//
340			//
341			return(true);
342			//
343			}
344
345			Bool_t CaloFranzini::LoadBin(){
346			//
347	mocchiut	1.1	TArrayI numbin = (TArrayI)file->Get("nbinenergy");
348			if ( !numbin ){
349	mocchiut	1.4	printf(" ERROR: cannot read number of bins from file ");
350	mocchiut	1.1	return(false);
351			};
352			nbin = (Int_t)numbin->At(0);
353			if ( nbin <= 0 ){
354	mocchiut	1.4	printf(" ERROR: cannot work with 0 energy bins from file ");
355	mocchiut	1.1	return(false);
356			};
357			//
358			brig = (TArrayF*)file->Get("binrig");
359			if ( !brig ){
360	mocchiut	1.4	printf(" ERROR: cannot read rigidity binning from file ");
361	mocchiut	1.1	return(false);
362			};
363			//
364	mocchiut	1.4	brigm=(TArrayF*)file->Get("binrigmean");
365			if ( !brigm ){
366			printf(" ERROR: cannot read mean rigidity binning from file ");
367			return(false);
368			};
369			//
370			for (Int_t i=0;i<17;i++){
371			TString name = Form("qplmeann%i",i);
372			qplmean[i] = (TArrayF*)file->Get(name.Data());
373			if ( !qplmean[i] ){
374			printf(" ERROR: cannot read average from file ");
375			return(false);
376			};
377			};
378			//
379	mocchiut	1.1	return(true);
380	mocchiut	1.4	}
381
382			Bool_t CaloFranzini::LoadMatrices(){
383	mocchiut	1.1	//
384	mocchiut	1.4	for (Int_t i=0;i<17;i++){
385			TString name1 = Form("matrixn%i",i);
386			hmat[i] = (TMatrixD*)file->Get(name1.Data());
387			};
388			//
389			return(true);
390	mocchiut	1.1	}
391
392			TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){
393			//
394	mocchiut	1.4	Int_t mv = 0;
395	mocchiut	1.1	for (Int_t i = 0; i<nbin-1; i++){
396			if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
397	mocchiut	1.4	mv = i;
398	mocchiut	1.1	break;
399			};
400			};
401			if ( rig < brig->At(0) ){
402			printf(" WARNING: Event with rigidity lower than the first covariance matrix bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
403	mocchiut	1.4	mv = 0;
404	mocchiut	1.1	};
405	mocchiut	1.3	if ( rig >= brig->At(nbin-1) ){
406			printf(" WARNING: Event with rigidity higher than the last covariance matrix bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
407	mocchiut	1.4	mv = nbin-2;
408	mocchiut	1.1	};
409			//
410	mocchiut	1.4	return(hmat[mv]);
411	mocchiut	1.1	//
412			}
413
414
415			TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){
416			//
417	mocchiut	1.4	Int_t mv=0;
418	mocchiut	1.1	for (Int_t i = 0; i<nbin-1; i++){
419			if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
420	mocchiut	1.4	mv = i;
421	mocchiut	1.1	break;
422			};
423			};
424			if ( rig < brig->At(0) ){
425			printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
426	mocchiut	1.4	mv = 0;
427	mocchiut	1.1	};
428	mocchiut	1.3	if ( rig >= brig->At(nbin-1) ){
429			printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
430	mocchiut	1.4	mv=nbin-2;
431	mocchiut	1.1	};
432			//
433	mocchiut	1.4	return(qplmean[mv]);
434	mocchiut	1.1	//
435	mocchiut	1.4	}
436
437			Float_t CaloFranzini::GetAverageAt(Int_t plane, Float_t rig){
438			//
439			Int_t therigb = 0;
440			for (Int_t i = 0; i<nbin-2; i++){
441			if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
442			therigb = i;
443			break;
444			};
445			};
446			//
447			Float_t minrig;
448			Float_t maxrig;
449			//
450			//
451			if ( rig < brigm->At(0) ){
452			if ( rig < brig->At(0) ){
453			printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
454			};
455			therigb = 0;
456			};
457			if ( rig >= brigm->At(nbin-2) ){
458			if ( rig >= brig->At(nbin-2) ) {
459			printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
460			};
461			therigb = nbin - 3;
462			};
463			//
464			minrig = brigm->At(therigb);
465			maxrig = brigm->At(therigb+1);
466			//
467			Float_t minene = (*qplmean[therigb])[plane];
468			Float_t maxene = (*qplmean[therigb+1])[plane];
469			//
470			if ( maxrig == minrig ){
471			printf("Unrecoverable ERROR! Matrix will be screwed up... \n");
472			return(0.);
473			};
474			Float_t b = log(maxene/minene)/(maxrig-minrig);
475			Float_t a = minene/exp(minrig*b);
476			Float_t ave = aexp(brig);
477			//
478			return(ave);
479	mocchiut	1.1	//
480			}
481	mocchiut	1.3
482	mocchiut	1.4	Float_t CaloFranzini::GetHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
483			Int_t therigb = 0;
484			for (Int_t i = 0; i<nbin-2; i++){
485			if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
486			therigb = i;
487			break;
488			};
489			};
490			//
491			Float_t minrig;
492			Float_t maxrig;
493			//
494			if ( rig < brigm->At(0) ){
495			if ( rig < brig->At(0) ){
496			printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
497			};
498			therigb = 0;
499			};
500			if ( rig >= brigm->At(nbin-2) ){
501			if ( rig >= brig->At(nbin-2) ) {
502			printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
503			};
504			therigb = nbin - 3;
505			};
506			//
507			minrig = brigm->At(therigb);
508			maxrig = brigm->At(therigb+1);
509			//
510			Float_t minene = (*hmat[therigb])[iindex][jindex];
511			Float_t maxene = (*hmat[therigb+1])[iindex][jindex];
512			//
513			Float_t b = log(maxene/minene)/(maxrig-minrig);
514			Float_t a = minene/exp(minrig*b);
515			Float_t ave = aexp(brig);
516			//
517			return(ave);
518			//
519			}
520	mocchiut	1.3
521			void CaloFranzini::DrawLongAverage(Float_t rig){
522			//
523			TArrayF *ll = this->LoadLongAverage(rig);
524			//
525			gStyle->SetLabelSize(0.04);
526			gStyle->SetNdivisions(510,"XY");
527			//
528			TString hid = Form("cfralongvyvx");
529			TCanvas tcf = dynamic_cast<TCanvas>(gDirectory->FindObject(hid));
530			if ( tcf ){
531			tcf->cd();
532			} else {
533			tcf = new TCanvas(hid,hid);
534			};
535			//
536			TString thid = Form("hfralongvyvx");
537			TH1F thf = dynamic_cast<TH1F>(gDirectory->FindObject(thid));
538			if ( thf ) thf->Delete();
539			thf = new TH1F(thid,thid,44,-0.5,43.5);
540			tcf->cd();
541			Float_t qpl[43];
542			for (Int_t st=0;st<43;st++){
543			qpl[st] = ll->At(st);
544	mocchiut	1.4	printf("st %i qpl %f\n",st,qpl[st]);
545	mocchiut	1.3	};
546			for (Int_t st=0;st<44;st++){
547			if ( st == 37 ){
548			thf->Fill(st,0.);
549			} else {
550			Int_t ss = st;
551			if ( st > 37 ) ss--;
552			thf->Fill(st,qpl[ss]);
553			};
554			};
555			thf->Draw();
556			tcf->Modified();
557			tcf->Update();
558			//
559			gStyle->SetLabelSize(0);
560			gStyle->SetNdivisions(1,"XY");
561			//
562			};