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;
  //
  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, 2*22*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(" longtzeta                        :.. %f\n",longtzeta);  
  printf(" fulltzeta                        :.. %f\n",fulltzeta);  
  printf("========================================================================\n");
  //
}

void CaloFranzini::Delete(){
  //
  if ( file ) file->Close();
  //
  Clear();
  //
}


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;
    qplane[nplane] += mip;
    //
  };
  //
  //
  //
  if ( dolong ){
    //
    degfre = 44;
    //
    for (Int_t i=0; i < 22; i++){
      if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){
        degfre = 2 * i;
        break;
      };
      if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){
      degfre = 1 + 2 * i;
      break;
      };
    };
    //
    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 ( 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"); 
  } 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){
  TArrayI nbi(1, &numbin);
  file->WriteObject(&nbi, "nbinenergy");
}

void CaloFranzini::WriteRigBin(TArrayF *rigbin){
  file->WriteObject(&rigbin, "binrig");
}

void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){
  TString name = Form("qplmeann%i",bin);
  file->WriteObject(&qpl,name.Data());
}

void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
  TString name = Form("matrixn%i",bin);
  file->WriteObject(&(*mat),name.Data());
}

void CaloFranzini::CloseMatrixFile(){
  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) ){
    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());
  };
  //
  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) ){
    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());
  };
  //
  TArrayF *qplmean = (TArrayF*)file->Get(name.Data());
  //
  return(qplmean);
  //
}
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			Clear();
35			//
36			}
37
38			void CaloFranzini::Clear(){
39			//
40			OBT = 0;
41			PKT = 0;
42			atime = 0;
43	mocchiut	1.2	longtzeta = -100.;
44			fulltzeta = -100.;
45			degfre = 0;
46	mocchiut	1.1	memset(estrip, 0, 22296*sizeof(Float_t));
47			memset(qplane, 0, 222sizeof(Float_t));
48			//
49			}
50
51			void CaloFranzini::Print(){
52			//
53			Process();
54			//
55			printf("========================================================================\n");
56			printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
57			printf(" debug [debug flag]:.. %i\n",debug);
58			printf(" longtzeta :.. %f\n",longtzeta);
59			printf(" fulltzeta :.. %f\n",fulltzeta);
60			printf("========================================================================\n");
61			//
62			}
63
64			void CaloFranzini::Delete(){
65	mocchiut	1.2	//
66	mocchiut	1.1	if ( file ) file->Close();
67	mocchiut	1.2	//
68	mocchiut	1.1	Clear();
69	mocchiut	1.2	//
70			}
71
72
73			void CaloFranzini::Process(){
74			this->Process(0);
75			}
76
77			void CaloFranzini::Process(Int_t itr){
78			//
79			if ( !L2 ){
80			printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
81			printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
82			return;
83			};
84			//
85			if ( !nbin \|\| !file \|\| !brig ){
86			printf(" ERROR: it seems covariance matrix file has not been opened (CaloFranzini::Open()) \n");
87			printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
88			return;
89			};
90			//
91			Bool_t newentry = false;
92			//
93			if ( L2->IsORB() ){
94			if ( L2->GetOrbitalInfo()->pkt_num != PKT \|\| L2->GetOrbitalInfo()->OBT != OBT \|\| L2->GetOrbitalInfo()->absTime != atime \|\| itr != sntr ){
95			newentry = true;
96			OBT = L2->GetOrbitalInfo()->OBT;
97			PKT = L2->GetOrbitalInfo()->pkt_num;
98			atime = L2->GetOrbitalInfo()->absTime;
99			sntr = itr;
100			};
101			} else {
102			newentry = true;
103			};
104			//
105			if ( !newentry ) return;
106			//
107			// Some variables
108			//
109			if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
110			//
111			this->Clear();
112			//
113			Float_t rig = L2->GetTrack(itr)->GetTrkTrack()->GetRigidity();
114			//
115			// Fill the estrip matrix
116			//
117			Int_t nplane = 0;
118			Int_t view = 0;
119			Int_t plane = 0;
120			Int_t strip = 0;
121			Float_t mip = 0.;
122			for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
123			//
124			mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
125			//
126			estrip[view][plane][strip] = mip;
127			nplane = 1 - view + 2 * plane;
128			qplane[nplane] += mip;
129			//
130			};
131			//
132			//
133			//
134			if ( dolong ){
135			//
136			degfre = 44;
137			//
138			for (Int_t i=0; i < 22; i++){
139			if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){
140			degfre = 2 * i;
141			break;
142			};
143			if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){
144			degfre = 1 + 2 * i;
145			break;
146			};
147			};
148			//
149			if ( degfre > 0 ){
150			TArrayF *qplmean = this->LoadLongAverage(rig);
151			TMatrixD *lmatrix = this->LoadCovarianceMatrix(rig);
152			for (Int_t i = 0; i < degfre; i++){
153			for (Int_t j = 0; j < degfre; j++){
154			longtzeta += (qplane[i]-qplmean->At(i)) * (lmatrix)[i][j] (qplane[j]-qplmean->At(j));
155			};
156			};
157			};
158			};
159			if ( dofull ){
160			printf(" ERROR: NOT IMPLEMENTED YET\n");
161	mocchiut	1.1	printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
162			};
163			//
164			if ( debug ) this->Print();
165			if ( debug ) printf(" exit \n");
166	mocchiut	1.2	}
167
168
169			Float_t CaloFranzini::GetNormLongTZeta(){
170			Process();
171			Float_t normz = 0.;
172			if ( degfre > 0 ) normz = longtzeta/(Float_t)degfre;
173			return normz;
174			}
175
176			Float_t CaloFranzini::GetNormFullTZeta(){
177			Process();
178			Float_t normz = 0.;
179			if ( degfre > 0 ) normz = fulltzeta/(Float_t)degfre;
180			return normz;
181	mocchiut	1.1	}
182
183			Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){
184			//
185			file = new TFile(matrixfile.Data(),"READ");
186			//
187			if ( !file \|\| file->IsZombie() ){
188			file = new TFile(matrixfile.Data(),"RECREATE");
189			} else {
190			printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
191			return(false);
192			};
193			//
194			return(true);
195			//
196			}
197
198			void CaloFranzini::WriteNumBin(Int_t numbin){
199			TArrayI nbi(1, &numbin);
200			file->WriteObject(&nbi, "nbinenergy");
201			}
202
203			void CaloFranzini::WriteRigBin(TArrayF *rigbin){
204			file->WriteObject(&rigbin, "binrig");
205			}
206
207			void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){
208			TString name = Form("qplmeann%i",bin);
209			file->WriteObject(&qpl,name.Data());
210			}
211
212			void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
213			TString name = Form("matrixn%i",bin);
214			file->WriteObject(&(*mat),name.Data());
215			}
216
217			void CaloFranzini::CloseMatrixFile(){
218			file->Close();
219			}
220
221
222			Bool_t CaloFranzini::Open(TString matrixfile){
223			//
224			// find matrix file
225			//
226			if ( !strcmp(matrixfile.Data(),"") ){
227			if (dolong) matrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_longel.root";
228			if (dofull) matrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_fullel.root";
229			};
230			//
231			file = new TFile(matrixfile.Data(),"READ");
232			//
233			if ( !file \|\| file->IsZombie() ){
234			printf(" ERROR: cannot open file %s \n",matrixfile.Data());
235			return(false);
236			};
237			//
238			TArrayI numbin = (TArrayI)file->Get("nbinenergy");
239			if ( !numbin ){
240			printf(" ERROR: cannot read number of bins from file %s \n",matrixfile.Data());
241			return(false);
242			};
243			nbin = (Int_t)numbin->At(0);
244			if ( nbin <= 0 ){
245			printf(" ERROR: cannot work with 0 energy bins (from file %s) \n",matrixfile.Data());
246			return(false);
247			};
248			//
249			brig = (TArrayF*)file->Get("binrig");
250			if ( !brig ){
251			printf(" ERROR: cannot read rigidity binning from file %s \n",matrixfile.Data());
252			return(false);
253			};
254			//
255			return(true);
256			//
257			}
258
259			TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){
260			//
261			TString name;
262			for (Int_t i = 0; i<nbin-1; i++){
263			if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
264			name = Form("matrixn%i",i);
265			break;
266			};
267			};
268			if ( rig < brig->At(0) ){
269			printf(" WARNING: Event with rigidity lower than the first covariance matrix bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
270			name = "matrixn0";
271			printf(" Using matrix %s \n",name.Data());
272			};
273			if ( rig >= brig->At(nbin) ){
274			printf(" WARNING: Event with rigidity higher than the last covariance matrix bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin));
275			name = Form("matrixn%i",nbin-1);
276			printf(" Using matrix %s \n",name.Data());
277			};
278			//
279			TMatrixD matrix = (TMatrixD)file->Get(name.Data());
280			//
281			return(matrix);
282			//
283			}
284
285
286			TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){
287			//
288			TString name;
289			for (Int_t i = 0; i<nbin-1; i++){
290			if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
291			name = Form("qplmeann%i",i);
292			break;
293			};
294			};
295			if ( rig < brig->At(0) ){
296			printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
297			name = "qplmeann0";
298			printf(" Using qplmean %s \n",name.Data());
299			};
300			if ( rig >= brig->At(nbin) ){
301			printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin));
302			name = Form("qplmeann%i",nbin-1);
303			printf(" Using qplmean %s \n",name.Data());
304			};
305			//
306			TArrayF qplmean = (TArrayF)file->Get(name.Data());
307			//
308			return(qplmean);
309			//
310			}