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;
  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 ){
      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){  
  //
  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);
  //
}

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);
  };
  //
  TArrayI *numbin = (TArrayI*)file->Get("nbinenergy");
  if ( !numbin ){
    printf(" ERROR: cannot read number of bins from file %s \n",matrixfile.Data());
    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());
    return(false);
  };
  //
  brig = (TArrayF*)file->Get("binrig");
  if ( !brig ){
    printf(" ERROR: cannot read rigidity binning from file %s \n",matrixfile.Data());
    return(false);
  };
  //
  return(true);
  //
}

TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){
  //
  TString name;
  for (Int_t i = 0; i<nbin-1; i++){
    if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
      name = Form("matrixn%i",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());
  };
  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());
  };
  //
  TMatrixD *matrix = (TMatrixD*)file->Get(name.Data());
  //
  return(matrix);
  //
}


TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){
  //
  TString name;
  for (Int_t i = 0; i<nbin-1; i++){
    if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
      name = Form("qplmeann%i",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());
  };
  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());
  };
  //
  TArrayF *qplmean = (TArrayF*)file->Get(name.Data());
  //
  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");
  //
};
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			nbin = 0;
22			//
23			L2 = l2p;
24			//
25			if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
26			//
27			// Default variables
28			//
29			debug = false;
30			dolong = true;
31			dofull = false;
32	mocchiut	1.2	sntr = 0;
33	mocchiut	1.1	//
34	mocchiut	1.3	sel = true;
35			cont = false;
36			N = 0;
37			NC = 43;
38			//
39			mask18b = -1;
40			//
41	mocchiut	1.1	Clear();
42			//
43			}
44
45			void CaloFranzini::Clear(){
46			//
47			OBT = 0;
48			PKT = 0;
49			atime = 0;
50	mocchiut	1.2	longtzeta = -100.;
51			fulltzeta = -100.;
52			degfre = 0;
53	mocchiut	1.1	memset(estrip, 0, 22296*sizeof(Float_t));
54	mocchiut	1.3	memset(qplane, 0, 43*sizeof(Float_t));
55	mocchiut	1.1	//
56			}
57
58			void CaloFranzini::Print(){
59			//
60			Process();
61			//
62			printf("========================================================================\n");
63			printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
64			printf(" debug [debug flag]:.. %i\n",debug);
65	mocchiut	1.3	printf(" degree of freedom :.. %i\n",this->GetDegreeOfFreedom());
66	mocchiut	1.1	printf(" longtzeta :.. %f\n",longtzeta);
67	mocchiut	1.3	printf(" longtzeta normalized :.. %f\n",this->GetNormLongTZeta());
68			// printf(" fulltzeta :.. %f\n",fulltzeta);
69			// printf(" longtzeta normalized :.. %f\n",this->GetNormFullTZeta());
70	mocchiut	1.1	printf("========================================================================\n");
71			//
72			}
73
74			void CaloFranzini::Delete(){
75	mocchiut	1.2	//
76	mocchiut	1.1	if ( file ) file->Close();
77	mocchiut	1.2	//
78	mocchiut	1.1	Clear();
79	mocchiut	1.2	//
80			}
81
82	mocchiut	1.3	void CaloFranzini::SetNoWpreSampler(Int_t n){
83			if ( NC+n < 23 ){
84			N = n;
85			} else {
86			printf(" ERROR! Calorimeter is made of 22 W planes\n");
87			printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
88			printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
89			NC = 43;
90			N = 0;
91			};
92			}
93
94			void CaloFranzini::SetNoWcalo(Int_t n){
95			if ( N+n < 23 ){
96			NC = n*2;
97			if ( NC >37 ) NC--;
98			} else {
99			printf(" ERROR! Calorimeter is made of 22 W planes\n");
100			printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
101			printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
102			NC = 43;
103			N = 0;
104			};
105			}
106	mocchiut	1.2
107			void CaloFranzini::Process(){
108			this->Process(0);
109			}
110
111			void CaloFranzini::Process(Int_t itr){
112			//
113			if ( !L2 ){
114			printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
115			printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
116			return;
117			};
118			//
119			if ( !nbin \|\| !file \|\| !brig ){
120			printf(" ERROR: it seems covariance matrix file has not been opened (CaloFranzini::Open()) \n");
121			printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
122			return;
123			};
124			//
125			Bool_t newentry = false;
126			//
127			if ( L2->IsORB() ){
128			if ( L2->GetOrbitalInfo()->pkt_num != PKT \|\| L2->GetOrbitalInfo()->OBT != OBT \|\| L2->GetOrbitalInfo()->absTime != atime \|\| itr != sntr ){
129			newentry = true;
130			OBT = L2->GetOrbitalInfo()->OBT;
131			PKT = L2->GetOrbitalInfo()->pkt_num;
132			atime = L2->GetOrbitalInfo()->absTime;
133			sntr = itr;
134			};
135			} else {
136			newentry = true;
137			};
138			//
139			if ( !newentry ) return;
140			//
141			// Some variables
142			//
143			if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
144			//
145			this->Clear();
146			//
147			Float_t rig = L2->GetTrack(itr)->GetTrkTrack()->GetRigidity();
148			//
149			// Fill the estrip matrix
150			//
151			Int_t nplane = 0;
152			Int_t view = 0;
153			Int_t plane = 0;
154			Int_t strip = 0;
155			Float_t mip = 0.;
156			for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
157			//
158			mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
159			//
160			estrip[view][plane][strip] = mip;
161	mocchiut	1.3	//
162			nplane = 1 - view + 2 * (plane - N);
163			if ( nplane > (37-(2*N)) ) nplane--;
164			//
165			if ( plane == (18+N) ) mip = 0.;
166			if ( nplane > -1 ) qplane[nplane] += mip;
167	mocchiut	1.2	//
168			};
169			//
170			//
171			//
172			if ( dolong ){
173			//
174	mocchiut	1.3	if ( cont ){
175			for (Int_t i=0; i<22; i++){
176			if ( i == (18+N) ){
177			mask18b = 1 + 2 * (i - N);
178			break;
179			};
180			};
181			};
182			//
183			if ( sel ){
184			for (Int_t i=0; i<22; i++){
185			if ( i == (18-N) ){
186			mask18b = 1 + 2 * (i - N);
187			break;
188			};
189			};
190			};
191			//
192			if ( mask18b == 37 ) mask18b = -1;
193			//
194			Int_t dgf = 43;
195	mocchiut	1.2	//
196			for (Int_t i=0; i < 22; i++){
197			if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){
198	mocchiut	1.3	dgf = 2 * i;
199	mocchiut	1.2	break;
200			};
201			if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){
202	mocchiut	1.3	dgf = 1 + 2 * i;
203	mocchiut	1.2	break;
204			};
205			};
206			//
207	mocchiut	1.3	if ( dgf < 43 && dgf > 37 ) dgf--;
208			//
209			degfre = TMath::Min(dgf,NC);
210			//
211	mocchiut	1.2	if ( degfre > 0 ){
212			TArrayF *qplmean = this->LoadLongAverage(rig);
213			TMatrixD *lmatrix = this->LoadCovarianceMatrix(rig);
214			for (Int_t i = 0; i < degfre; i++){
215	mocchiut	1.3	if ( i != mask18b ){
216			for (Int_t j = 0; j < degfre; j++){
217			if ( j != mask18b ){
218			longtzeta += (qplane[i]-qplmean->At(i)) * (lmatrix)[i][j] (qplane[j]-qplmean->At(j));
219			};
220			};
221	mocchiut	1.2	};
222			};
223			};
224			};
225			if ( dofull ){
226			printf(" ERROR: NOT IMPLEMENTED YET\n");
227	mocchiut	1.1	printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
228			};
229			//
230			if ( debug ) this->Print();
231			if ( debug ) printf(" exit \n");
232	mocchiut	1.2	}
233
234
235			Float_t CaloFranzini::GetNormLongTZeta(){
236			Process();
237			Float_t normz = 0.;
238			if ( degfre > 0 ) normz = longtzeta/(Float_t)degfre;
239			return normz;
240			}
241
242			Float_t CaloFranzini::GetNormFullTZeta(){
243			Process();
244			Float_t normz = 0.;
245			if ( degfre > 0 ) normz = fulltzeta/(Float_t)degfre;
246			return normz;
247	mocchiut	1.1	}
248
249			Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){
250			//
251			file = new TFile(matrixfile.Data(),"READ");
252			//
253			if ( !file \|\| file->IsZombie() ){
254			file = new TFile(matrixfile.Data(),"RECREATE");
255	mocchiut	1.3	printf(" Create file %s \n",file->GetName());
256	mocchiut	1.1	} else {
257			printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
258			return(false);
259			};
260			//
261			return(true);
262			//
263			}
264
265			void CaloFranzini::WriteNumBin(Int_t numbin){
266	mocchiut	1.3	file->cd();
267	mocchiut	1.1	TArrayI nbi(1, &numbin);
268			file->WriteObject(&nbi, "nbinenergy");
269			}
270
271			void CaloFranzini::WriteRigBin(TArrayF *rigbin){
272	mocchiut	1.3	file->cd();
273			// rigbin->Write("binrig");
274			file->WriteObject(&(*rigbin), "binrig");
275	mocchiut	1.1	}
276
277			void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){
278	mocchiut	1.3	file->cd();
279	mocchiut	1.1	TString name = Form("qplmeann%i",bin);
280	mocchiut	1.3	file->WriteObject(&(*qpl),name.Data());
281			}
282
283			void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){
284			file->cd();
285			TString name = Form("matrixn%i",bin);
286			// mat.Write(name.Data());
287			file->WriteObject(&mat,name.Data());
288	mocchiut	1.1	}
289
290			void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
291	mocchiut	1.3	file->cd();
292			TString name = Form("origmatrixn%i",bin);
293			// mat.Write(name.Data());
294	mocchiut	1.1	file->WriteObject(&(*mat),name.Data());
295			}
296
297			void CaloFranzini::CloseMatrixFile(){
298	mocchiut	1.3	file->cd();
299	mocchiut	1.1	file->Close();
300			}
301
302
303			Bool_t CaloFranzini::Open(TString matrixfile){
304			//
305			// find matrix file
306			//
307			if ( !strcmp(matrixfile.Data(),"") ){
308			if (dolong) matrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_longel.root";
309			if (dofull) matrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_fullel.root";
310			};
311			//
312			file = new TFile(matrixfile.Data(),"READ");
313			//
314			if ( !file \|\| file->IsZombie() ){
315			printf(" ERROR: cannot open file %s \n",matrixfile.Data());
316			return(false);
317			};
318			//
319			TArrayI numbin = (TArrayI)file->Get("nbinenergy");
320			if ( !numbin ){
321			printf(" ERROR: cannot read number of bins from file %s \n",matrixfile.Data());
322			return(false);
323			};
324			nbin = (Int_t)numbin->At(0);
325			if ( nbin <= 0 ){
326			printf(" ERROR: cannot work with 0 energy bins (from file %s) \n",matrixfile.Data());
327			return(false);
328			};
329			//
330			brig = (TArrayF*)file->Get("binrig");
331			if ( !brig ){
332			printf(" ERROR: cannot read rigidity binning from file %s \n",matrixfile.Data());
333			return(false);
334			};
335			//
336			return(true);
337			//
338			}
339
340			TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){
341			//
342			TString name;
343			for (Int_t i = 0; i<nbin-1; i++){
344			if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
345			name = Form("matrixn%i",i);
346			break;
347			};
348			};
349			if ( rig < brig->At(0) ){
350			printf(" WARNING: Event with rigidity lower than the first covariance matrix bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
351			name = "matrixn0";
352			printf(" Using matrix %s \n",name.Data());
353			};
354	mocchiut	1.3	if ( rig >= brig->At(nbin-1) ){
355			printf(" WARNING: Event with rigidity higher than the last covariance matrix bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
356			name = Form("matrixn%i",nbin-2);
357	mocchiut	1.1	printf(" Using matrix %s \n",name.Data());
358			};
359			//
360			TMatrixD matrix = (TMatrixD)file->Get(name.Data());
361			//
362			return(matrix);
363			//
364			}
365
366
367			TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){
368			//
369			TString name;
370			for (Int_t i = 0; i<nbin-1; i++){
371			if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
372			name = Form("qplmeann%i",i);
373			break;
374			};
375			};
376			if ( rig < brig->At(0) ){
377			printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
378			name = "qplmeann0";
379			printf(" Using qplmean %s \n",name.Data());
380			};
381	mocchiut	1.3	if ( rig >= brig->At(nbin-1) ){
382			printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
383			name = Form("qplmeann%i",nbin-2);
384	mocchiut	1.1	printf(" Using qplmean %s \n",name.Data());
385			};
386			//
387			TArrayF qplmean = (TArrayF)file->Get(name.Data());
388			//
389			return(qplmean);
390			//
391			}
392	mocchiut	1.3
393
394
395			void CaloFranzini::DrawLongAverage(Float_t rig){
396			//
397			TArrayF *ll = this->LoadLongAverage(rig);
398			//
399			gStyle->SetLabelSize(0.04);
400			gStyle->SetNdivisions(510,"XY");
401			//
402			TString hid = Form("cfralongvyvx");
403			TCanvas tcf = dynamic_cast<TCanvas>(gDirectory->FindObject(hid));
404			if ( tcf ){
405			tcf->cd();
406			} else {
407			tcf = new TCanvas(hid,hid);
408			};
409			//
410			TString thid = Form("hfralongvyvx");
411			TH1F thf = dynamic_cast<TH1F>(gDirectory->FindObject(thid));
412			if ( thf ) thf->Delete();
413			thf = new TH1F(thid,thid,44,-0.5,43.5);
414			tcf->cd();
415			// Int_t pp=0;
416			Float_t qpl[43];
417			for (Int_t st=0;st<43;st++){
418			qpl[st] = ll->At(st);
419			};
420			for (Int_t st=0;st<44;st++){
421			if ( st == 37 ){
422			thf->Fill(st,0.);
423			} else {
424			Int_t ss = st;
425			if ( st > 37 ) ss--;
426			thf->Fill(st,qpl[ss]);
427			};
428			};
429			thf->Draw();
430			tcf->Modified();
431			tcf->Update();
432			//
433			gStyle->SetLabelSize(0);
434			gStyle->SetNdivisions(1,"XY");
435			//
436			};