CaloFranzini/src/Calib.cpp

#include <stdlib.h>
#include <iostream>
#include <iomanip>
//
#include <TString.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TMatrixD.h>
#include <TArrayF.h>
#include <TArrayI.h>
//
#include <PamLevel2.h>
#include <CaloFranzini.h>
//
using namespace std;
//
extern Bool_t MATRIX;
CaloFranzini *cf;
Int_t nbin;
Float_t rig[18];
Float_t rmean[17];
//Float_t qqplane[17][43];
//Int_t nnqplane[17][43];

TArrayF *qplane[17];
TArrayI *nqplane[17];
TMatrixD *matrix[17];
TMatrixD *nmat[17];

//===============================================================================
bool Select( PamLevel2* event ){

  //---------------------------------------------------------
  // single track
  //---------------------------------------------------------
  if( event->GetTrkLevel2()->GetNTracks()!=1 ) return false; 
  PamTrack *track = event->GetTrack(0);
  if(!track)return false;

  //------------------------------------------------------------------
  // tracker pre-selection
  //------------------------------------------------------------------
  TrkTrack *trk = track->GetTrkTrack();

  bool TRACK__OK = false;
  if( 
     trk->chi2 >0     &&
     trk->GetNX()>=4  &&
     trk->GetNY()>=3  &&
     trk->GetLeverArmX()>=5 &&
     true ) TRACK__OK = true;

  if( !TRACK__OK )return false;

  //------------------------------------------------------------------
  // TOF pre-selection
  //------------------------------------------------------------------
  bool TOF__OK = false;
  if(
     event->GetToFLevel2()->GetNHitPaddles(0) == 1 &&
     event->GetToFLevel2()->GetNHitPaddles(1) == 1 &&
     event->GetToFLevel2()->GetNHitPaddles(2) == 1 &&
     event->GetToFLevel2()->GetNHitPaddles(3) == 1 &&
     event->GetToFLevel2()->GetNHitPaddles(4) >= 1 &&
     event->GetToFLevel2()->GetNHitPaddles(5) >= 1 &&
     event->GetToFLevel2()->npmt() <= 18           &&
     !event->GetAcLevel2()->CARDhit() &&
     !event->GetAcLevel2()->CAThit() &&
     true ) TOF__OK = true;
  if( !TOF__OK )return false;
  //------------------------------------------------------
  // no albedo
  //------------------------------------------------------
  if(    track->GetToFTrack()->beta[12]<=0.2 ||
         track->GetToFTrack()->beta[12] >= 1.5 ) return false;

  //------------------------------------------------------
  bool CUT1 = false;
  if( 
     trk->nstep<100   &&
     trk->GetRigidity()<400.   &&
     trk->GetRigidity()>0.1   &&
     trk->GetDeflection()<0.   &&
     trk->resx[0]<0.001 && 
     trk->resx[5]<0.001 && 
     track->IsSolved() &&
     trk->IsInsideCavity() &&
     true ) CUT1 = true;
  //------------------------------------------------------
  if( !CUT1 )return false;
  Float_t defl = trk->GetDeflection();
  float p0 = 1.111588e+00;
  float p1 = 1.707656e+00;
  float p2 = 1.489693e-01;
  float chi2m025 = p0 + fabs(defl)*p1 + defl*defl*p2;
  
  float def_0      = 0.07;
  float chi2m025_0 = p0 + fabs(def_0)*p1 + def_0*def_0*p2;
  
  //  int   nchi2cut=5;
  float chi2cut=3.;
  float chi2m = pow( chi2m025-chi2m025_0+pow(chi2cut,0.25), 4.);
  bool CUT2 = false;
  if(
     track->GetTrkTrack()->chi2 < chi2m &&
     true ) CUT2 = true;  
  if ( !CUT2 ) return false;
  //------------------------------------------------------
  //
  // energy momentum match
  //
  Float_t qtotimp = event->GetCaloLevel2()->qtot / trk->GetRigidity();
  Float_t qcut2 = (-0.5 * trk->GetRigidity() + 150.) * 1.1;
  if ( qcut2 < 55. ) qcut2 = 55.;
  if ( qtotimp <= qcut2 ) return false;
  //
  for (Int_t i=0; i < 22; i++){
    if ( track->GetCaloTrack()->tibar[i][1] < 0 || track->GetCaloTrack()->tibar[i][0] < 0 ){
        return false;
    };
  };
  //
  if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false;
  //
  return true;
}
//===============================================================
// Create histograms
//
//
//
//
//
//===============================================================
void CreateHistos( PamLevel2* event , TString file){

  cf = new CaloFranzini(event);
  //
  if ( MATRIX ){
    cf->UpdateMatrixFile(file.Data());
    cf->LoadBin();
  } else {
    cf->CreateMatrixFile(file.Data());
  };
  //
  //
  nbin = 18;
  rig[0] = 0.1;
  rig[1] = 0.5;
  rig[2] = 1.;
  rig[3] = 1.5;
  rig[4] = 2.2;
  rig[5] = 3.;
  rig[6] = 4.;
  rig[7] = 5.;
  rig[8] = 6.;
  rig[9] = 8.;
  rig[10] = 10.;
  rig[11] = 15.;
  rig[12] = 25.;
  rig[13] = 35.;
  rig[14] = 50.;
  rig[15] = 100.;
  rig[16] = 200.;
  rig[17] = 400.;
  //
  memset(rmean, 0, 17*sizeof(Float_t));
  //
  for (Int_t i=0; i < 17 ; i++){
    matrix[i] = new TMatrixD(43,43);
    qplane[i] = new TArrayF(43);
    nqplane[i] = new TArrayI(43);
    nmat[i] = new TMatrixD(43,43);
  };
  //
}

//=============================================================== 
void FindAverage( PamLevel2* L2, int iev ){
  //
  Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
  //
  Int_t rbi = 0;
  for (Int_t i = 0; i<nbin-1; i++){ 
    if ( erig>=rig[i] && erig < rig[i+1] ){
      rbi = i;
      break;
    };
  };
  //
  if ( erig < rig[0] ) return;
  if ( erig >= rig[nbin-1] ) return;
  //
  rmean[rbi] += erig;
  //
  Int_t dgf = 43;
  //
  for (Int_t i=0; i<dgf; i++){
    (*nqplane[rbi])[i]++;
  };
  //
  // 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);
    //
    nplane = 1 - view + 2 * plane;
    if ( nplane > 37 ) nplane--; 
    if ( nplane < dgf ){
      (*qplane[rbi])[nplane] += mip;
    };
    //
  };
}

void CalculateAverage(){
  for (Int_t i=0; i<nbin-1; i++){
    if ( (*nqplane[i])[0] > 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0];
    for (Int_t j=0; j<43 ; j++){
      if ( (*nqplane[i])[j] > 0 ){  
        (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j];
      } else {
        (*qplane[i])[j] = 0.;
      };
      printf(" BIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(*qplane[i])[j],(*nqplane[i])[j]);
    };
  };
  
  cf->WriteNumBin(nbin);

  TArrayF *rigbin = new TArrayF(18, rig);
  cf->WriteRigBin(rigbin);

  TArrayF *rmeanbin = new TArrayF(17, rmean);
  TFile *file =  cf->GetFile();
  file->cd();
  file->WriteObject(&(*rmeanbin), "binrigmean");
  //
  for (Int_t i=0; i<nbin-1; i++){
    //
    cf->WriteLongMean(qplane[i], i);
    //
  };
  //
}

//=============================================================== 
void FindMatrix( PamLevel2* L2, int iev ){
  //
  Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
  //
  Int_t rbi = 0;
  for (Int_t i = 0; i<nbin-1; i++){
    if ( erig>=rig[i] && erig < rig[i+1] ){
      rbi = i;
      break;
    };
  };
  //
  if ( erig < rig[0] ) return;
  if ( erig >= rig[nbin-1] ) return;
  //
  Int_t dgf = 43;
//  printf("1qui \n");
  //
  for (Int_t i=0; i<dgf; i++){
    for (Int_t j=0; j<dgf; j++){
      (*nmat[rbi])[i][j] += 1.;
    };
  };
//  printf("2qui \n");
  //
  // Fill the estrip matrix
  //
  Int_t nplane = 0;
  Int_t view = 0;
  Int_t plane = 0;
  Int_t strip = 0;
  Float_t mip = 0.;
  Float_t hpl[43];
  memset(hpl,0,43*sizeof(Float_t));
  for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
    //
    mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
    //
    nplane = 1 - view + 2 * plane;
    if ( nplane > 37 ) nplane--; 
    if ( nplane < dgf ){
      hpl[nplane] += mip;
    };
    //
  };
  //
//  printf("3qui \n");
  //
  for (Int_t i=0; i<dgf; i++){
    for (Int_t j=0; j<dgf; j++){
      (*matrix[rbi])[i][j] += (hpl[i] - cf->GetAverageAt(i,erig)) * (hpl[j] - cf->GetAverageAt(j,erig));
    };
  };
//  printf("4qui \n");  
  //  gObjectTable->Print();
}

//===============================================================
// Save histograms
//
//
//
//
//
//===============================================================
void SaveHistos(){
  //
  if ( MATRIX ){
    //
    printf("Finished, calculating average and inverting matrices\n");
    //
    for (Int_t i=0; i<nbin-1; i++){
      //
      // determine the average matrix
      //    
      for (Int_t ii=0; ii<43; ii++){
        for (Int_t j=0; j<43; j++){
          if ( (*nmat[i])[ii][j] > 0. ){
            (*matrix[i])[ii][j] /= (*nmat[i])[ii][j];
          } else {
            (*matrix[i])[ii][j] = 0.;
          };
        };
      };
      //
      cf->WriteLongMatrix(matrix[i],i);
      //
      if ( matrix[i]->Determinant() == 0. ){
        printf("\n");
        for (Int_t ii=0; ii<43; ii++){
          for (Int_t j=0; j<43; j++){
            printf(" %.f",(*matrix[i])[ii][j]);   
          };
          printf("\n");
        };
        printf("\n");
        printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
      } else {
        Double_t det = 0.;
        TMatrixD invmatrix = (TMatrixD)(matrix[i]->Invert(&det));
        printf(" Bin %i determinant is %f \n",i,det);
        cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i);
      };
    };
    //
    printf(" done, closing file and exiting\n");
    //
  };
  //
  cf->CloseMatrixFile();
  //
  cf->Delete();
  //
}
1	mocchiut	1.1	#include <stdlib.h>
2			#include <iostream>
3			#include <iomanip>
4			//
5			#include <TString.h>
6			#include <TH1F.h>
7			#include <TH2F.h>
8			#include <TMatrixD.h>
9			#include <TArrayF.h>
10			#include <TArrayI.h>
11			//
12			#include <PamLevel2.h>
13			#include <CaloFranzini.h>
14			//
15			using namespace std;
16			//
17	mocchiut	1.2	extern Bool_t MATRIX;
18	mocchiut	1.1	CaloFranzini *cf;
19			Int_t nbin;
20			Float_t rig[18];
21	mocchiut	1.2	Float_t rmean[17];
22			//Float_t qqplane[17][43];
23			//Int_t nnqplane[17][43];
24
25			TArrayF *qplane[17];
26			TArrayI *nqplane[17];
27			TMatrixD *matrix[17];
28			TMatrixD *nmat[17];
29	mocchiut	1.1
30			//===============================================================================
31			bool Select( PamLevel2* event ){
32
33			//---------------------------------------------------------
34			// single track
35			//---------------------------------------------------------
36			if( event->GetTrkLevel2()->GetNTracks()!=1 ) return false;
37			PamTrack *track = event->GetTrack(0);
38			if(!track)return false;
39
40			//------------------------------------------------------------------
41			// tracker pre-selection
42			//------------------------------------------------------------------
43			TrkTrack *trk = track->GetTrkTrack();
44
45			bool TRACK__OK = false;
46			if(
47			trk->chi2 >0 &&
48			trk->GetNX()>=4 &&
49			trk->GetNY()>=3 &&
50			trk->GetLeverArmX()>=5 &&
51			true ) TRACK__OK = true;
52
53			if( !TRACK__OK )return false;
54
55			//------------------------------------------------------------------
56			// TOF pre-selection
57			//------------------------------------------------------------------
58			bool TOF__OK = false;
59			if(
60			event->GetToFLevel2()->GetNHitPaddles(0) == 1 &&
61			event->GetToFLevel2()->GetNHitPaddles(1) == 1 &&
62			event->GetToFLevel2()->GetNHitPaddles(2) == 1 &&
63			event->GetToFLevel2()->GetNHitPaddles(3) == 1 &&
64			event->GetToFLevel2()->GetNHitPaddles(4) >= 1 &&
65			event->GetToFLevel2()->GetNHitPaddles(5) >= 1 &&
66			event->GetToFLevel2()->npmt() <= 18 &&
67			!event->GetAcLevel2()->CARDhit() &&
68			!event->GetAcLevel2()->CAThit() &&
69			true ) TOF__OK = true;
70			if( !TOF__OK )return false;
71			//------------------------------------------------------
72			// no albedo
73			//------------------------------------------------------
74			if( track->GetToFTrack()->beta[12]<=0.2 \|\|
75			track->GetToFTrack()->beta[12] >= 1.5 ) return false;
76
77			//------------------------------------------------------
78			bool CUT1 = false;
79			if(
80			trk->nstep<100 &&
81			trk->GetRigidity()<400. &&
82			trk->GetRigidity()>0.1 &&
83			trk->GetDeflection()<0. &&
84			trk->resx[0]<0.001 &&
85			trk->resx[5]<0.001 &&
86			track->IsSolved() &&
87			trk->IsInsideCavity() &&
88			true ) CUT1 = true;
89			//------------------------------------------------------
90			if( !CUT1 )return false;
91			Float_t defl = trk->GetDeflection();
92			float p0 = 1.111588e+00;
93			float p1 = 1.707656e+00;
94			float p2 = 1.489693e-01;
95			float chi2m025 = p0 + fabs(defl)p1 + defldefl*p2;
96
97			float def_0 = 0.07;
98			float chi2m025_0 = p0 + fabs(def_0)p1 + def_0def_0*p2;
99
100			// int nchi2cut=5;
101			float chi2cut=3.;
102			float chi2m = pow( chi2m025-chi2m025_0+pow(chi2cut,0.25), 4.);
103			bool CUT2 = false;
104			if(
105			track->GetTrkTrack()->chi2 < chi2m &&
106			true ) CUT2 = true;
107			if ( !CUT2 ) return false;
108			//------------------------------------------------------
109			//
110			// energy momentum match
111			//
112			Float_t qtotimp = event->GetCaloLevel2()->qtot / trk->GetRigidity();
113			Float_t qcut2 = (-0.5 * trk->GetRigidity() + 150.) * 1.1;
114			if ( qcut2 < 55. ) qcut2 = 55.;
115			if ( qtotimp <= qcut2 ) return false;
116			//
117			for (Int_t i=0; i < 22; i++){
118			if ( track->GetCaloTrack()->tibar[i][1] < 0 \|\| track->GetCaloTrack()->tibar[i][0] < 0 ){
119			return false;
120			};
121			};
122			//
123			if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false;
124			//
125			return true;
126			}
127			//===============================================================
128			// Create histograms
129			//
130			//
131			//
132			//
133			//
134			//===============================================================
135			void CreateHistos( PamLevel2* event , TString file){
136
137			cf = new CaloFranzini(event);
138	mocchiut	1.2	//
139			if ( MATRIX ){
140			cf->UpdateMatrixFile(file.Data());
141			cf->LoadBin();
142			} else {
143			cf->CreateMatrixFile(file.Data());
144			};
145			//
146	mocchiut	1.1	//
147			nbin = 18;
148			rig[0] = 0.1;
149			rig[1] = 0.5;
150			rig[2] = 1.;
151			rig[3] = 1.5;
152			rig[4] = 2.2;
153			rig[5] = 3.;
154			rig[6] = 4.;
155			rig[7] = 5.;
156			rig[8] = 6.;
157			rig[9] = 8.;
158			rig[10] = 10.;
159			rig[11] = 15.;
160			rig[12] = 25.;
161			rig[13] = 35.;
162			rig[14] = 50.;
163			rig[15] = 100.;
164			rig[16] = 200.;
165			rig[17] = 400.;
166			//
167	mocchiut	1.2	memset(rmean, 0, 17*sizeof(Float_t));
168	mocchiut	1.1	//
169	mocchiut	1.2	for (Int_t i=0; i < 17 ; i++){
170	mocchiut	1.1	matrix[i] = new TMatrixD(43,43);
171			qplane[i] = new TArrayF(43);
172			nqplane[i] = new TArrayI(43);
173			nmat[i] = new TMatrixD(43,43);
174			};
175			//
176			}
177
178			//===============================================================
179			void FindAverage( PamLevel2* L2, int iev ){
180			//
181			Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
182			//
183			Int_t rbi = 0;
184	mocchiut	1.2	for (Int_t i = 0; i<nbin-1; i++){
185	mocchiut	1.1	if ( erig>=rig[i] && erig < rig[i+1] ){
186			rbi = i;
187			break;
188			};
189			};
190			//
191			if ( erig < rig[0] ) return;
192			if ( erig >= rig[nbin-1] ) return;
193			//
194			rmean[rbi] += erig;
195			//
196			Int_t dgf = 43;
197			//
198			for (Int_t i=0; i<dgf; i++){
199			(*nqplane[rbi])[i]++;
200			};
201			//
202			// Fill the estrip matrix
203			//
204			Int_t nplane = 0;
205			Int_t view = 0;
206			Int_t plane = 0;
207			Int_t strip = 0;
208			Float_t mip = 0.;
209			//
210			for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
211			//
212			mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
213			//
214			nplane = 1 - view + 2 * plane;
215			if ( nplane > 37 ) nplane--;
216			if ( nplane < dgf ){
217			(*qplane[rbi])[nplane] += mip;
218			};
219			//
220			};
221			}
222
223	mocchiut	1.2	void CalculateAverage(){
224			for (Int_t i=0; i<nbin-1; i++){
225			if ( (nqplane[i])[0] > 0 ) rmean[i] /= (Float_t)(nqplane[i])[0];
226			for (Int_t j=0; j<43 ; j++){
227			if ( (*nqplane[i])[j] > 0 ){
228			(qplane[i])[j] /= (Float_t)(nqplane[i])[j];
229			} else {
230			(*qplane[i])[j] = 0.;
231			};
232			printf(" BIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(qplane[i])[j],(nqplane[i])[j]);
233			};
234			};
235
236			cf->WriteNumBin(nbin);
237
238			TArrayF *rigbin = new TArrayF(18, rig);
239			cf->WriteRigBin(rigbin);
240
241			TArrayF *rmeanbin = new TArrayF(17, rmean);
242			TFile *file = cf->GetFile();
243			file->cd();
244			file->WriteObject(&(*rmeanbin), "binrigmean");
245			//
246			for (Int_t i=0; i<nbin-1; i++){
247			//
248			cf->WriteLongMean(qplane[i], i);
249			//
250			};
251			//
252			}
253
254	mocchiut	1.1	//===============================================================
255			void FindMatrix( PamLevel2* L2, int iev ){
256			//
257			Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
258			//
259			Int_t rbi = 0;
260			for (Int_t i = 0; i<nbin-1; i++){
261			if ( erig>=rig[i] && erig < rig[i+1] ){
262			rbi = i;
263			break;
264			};
265			};
266			//
267			if ( erig < rig[0] ) return;
268			if ( erig >= rig[nbin-1] ) return;
269			//
270			Int_t dgf = 43;
271	mocchiut	1.2	// printf("1qui \n");
272	mocchiut	1.1	//
273			for (Int_t i=0; i<dgf; i++){
274			for (Int_t j=0; j<dgf; j++){
275			(*nmat[rbi])[i][j] += 1.;
276			};
277			};
278	mocchiut	1.2	// printf("2qui \n");
279	mocchiut	1.1	//
280			// Fill the estrip matrix
281			//
282			Int_t nplane = 0;
283			Int_t view = 0;
284			Int_t plane = 0;
285			Int_t strip = 0;
286			Float_t mip = 0.;
287			Float_t hpl[43];
288			memset(hpl,0,43*sizeof(Float_t));
289			for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
290			//
291			mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
292			//
293			nplane = 1 - view + 2 * plane;
294			if ( nplane > 37 ) nplane--;
295			if ( nplane < dgf ){
296			hpl[nplane] += mip;
297			};
298			//
299			};
300			//
301	mocchiut	1.2	// printf("3qui \n");
302			//
303	mocchiut	1.1	for (Int_t i=0; i<dgf; i++){
304			for (Int_t j=0; j<dgf; j++){
305	mocchiut	1.2	(matrix[rbi])[i][j] += (hpl[i] - cf->GetAverageAt(i,erig)) (hpl[j] - cf->GetAverageAt(j,erig));
306	mocchiut	1.1	};
307			};
308	mocchiut	1.2	// printf("4qui \n");
309			// gObjectTable->Print();
310	mocchiut	1.1	}
311
312			//===============================================================
313			// Save histograms
314			//
315			//
316			//
317			//
318			//
319			//===============================================================
320			void SaveHistos(){
321	mocchiut	1.2	//
322			if ( MATRIX ){
323			//
324			printf("Finished, calculating average and inverting matrices\n");
325	mocchiut	1.1	//
326	mocchiut	1.2	for (Int_t i=0; i<nbin-1; i++){
327			//
328			// determine the average matrix
329			//
330			for (Int_t ii=0; ii<43; ii++){
331			for (Int_t j=0; j<43; j++){
332			if ( (*nmat[i])[ii][j] > 0. ){
333			(matrix[i])[ii][j] /= (nmat[i])[ii][j];
334			} else {
335			(*matrix[i])[ii][j] = 0.;
336			};
337	mocchiut	1.1	};
338			};
339	mocchiut	1.2	//
340			cf->WriteLongMatrix(matrix[i],i);
341			//
342			if ( matrix[i]->Determinant() == 0. ){
343			printf("\n");
344			for (Int_t ii=0; ii<43; ii++){
345			for (Int_t j=0; j<43; j++){
346			printf(" %.f",(*matrix[i])[ii][j]);
347			};
348			printf("\n");
349			};
350			printf("\n");
351			printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
352			} else {
353			Double_t det = 0.;
354			TMatrixD invmatrix = (TMatrixD)(matrix[i]->Invert(&det));
355			printf(" Bin %i determinant is %f \n",i,det);
356			cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i);
357	mocchiut	1.1	};
358			};
359	mocchiut	1.2	//
360			printf(" done, closing file and exiting\n");
361			//
362	mocchiut	1.1	};
363	mocchiut	1.2	//
364	mocchiut	1.1	cf->CloseMatrixFile();
365	mocchiut	1.2	//
366	mocchiut	1.1	cf->Delete();
367	mocchiut	1.2	//
368	mocchiut	1.1	}