GFlash/TB/myFunctions.cc

0001 // -------------------------------------------

0002 // useful functions

0003 // -------------------------------------------

0004
0005 // return the position [0->49] of the max amplitude crystal

0006 int maxAmplitInMatrix(double myMatrix[])
0007 {
0008   int maxXtal   = 999;
0009   double maxADC = -999.;
0010
0011   for (int icry=0; icry<49; icry++) {
0012     if (myMatrix[icry] > maxADC){
0013       maxADC  = myMatrix[icry];
0014       maxXtal = icry;
0015     }}
0016
0017   return maxXtal;
0018 }
0019
0020 // max amplitude xtal energy - passing the xtal number in the matrix

0021 double ene1x1_xtal(double mymatrix[], int xtal)
0022 {
0023   double E1x1 = 0.;
0024   if (mymatrix[xtal]<-50) { E1x1 = -1000.; }
0025   else { E1x1 = mymatrix[xtal]; }
0026   return E1x1;
0027 }
0028
0029 // 3x3 matrix energy around the max amplitude xtal

0030 double ene3x3_xtal(double mymatrix[], int xtal)
0031 {
0032   double E3x3 = 0.;
0033
0034   if ( (mymatrix[xtal-8]<-50) || (mymatrix[xtal-7]<-50)  || (mymatrix[xtal-6]<-50)  || (mymatrix[xtal-1]<-50) || (mymatrix[xtal]<-50)  || (mymatrix[xtal+1]<-50)  || (mymatrix[xtal+6]<-50) || (mymatrix[xtal+7]<-50)  || (mymatrix[xtal+8]<-50) )
0035     { E3x3 = -1000.; }
0036   else
0037     { E3x3 = mymatrix[xtal-8] + mymatrix[xtal-7] + mymatrix[xtal-6] + mymatrix[xtal-1] + mymatrix[xtal] + mymatrix[xtal+1] + mymatrix[xtal+6] + mymatrix[xtal+7] + mymatrix[xtal+8]; }
0038
0039   return E3x3;
0040 }
0041
0042 // 5x5 matrix energy around the max amplitude xtal

0043 double ene5x5_xtal(double mymatrix[], int xtal)
0044 {
0045   double E5x5 = 0.;
0046
0047   if( (mymatrix[xtal-16]<-50) || (mymatrix[xtal-15]<-50) || (mymatrix[xtal-14]<-50) || (mymatrix[xtal-13]<-50) || (mymatrix[xtal-12]<-50) || (mymatrix[xtal-9]<-50) || (mymatrix[xtal-8]<-50) || (mymatrix[xtal-7]<-50) || (mymatrix[xtal-6]<-50) || (mymatrix[xtal-5]<-50) || (mymatrix[xtal-2]<-50) || (mymatrix[xtal-1]<-50) || (mymatrix[xtal]<-50) || (mymatrix[xtal+1]<-50) || (mymatrix[xtal+2]<-50) || (mymatrix[xtal+5]<-50) || (mymatrix[xtal+6]<-50) || (mymatrix[xtal+7]<-50) || (mymatrix[xtal+8]<-50) || (mymatrix[xtal+9]<-50) || (mymatrix[xtal+12]<-50) || (mymatrix[xtal+13]<-50) || (mymatrix[xtal+14]<-50) || (mymatrix[xtal+15]<-50) || (mymatrix[xtal+16]<-50) )
0048     { E5x5 = -1000.; }
0049   else
0050     { E5x5 = mymatrix[xtal-16] + mymatrix[xtal-15] + mymatrix[xtal-14] + mymatrix[xtal-13] + mymatrix[xtal-12] + mymatrix[xtal-9] + mymatrix[xtal-8] + mymatrix[xtal-7] + mymatrix[xtal-6] + mymatrix[xtal-5] + mymatrix[xtal-2] + mymatrix[xtal-1] + mymatrix[xtal] + mymatrix[xtal+1] + mymatrix[xtal+2] + mymatrix[xtal+5] + mymatrix[xtal+6] + mymatrix[xtal+7] + mymatrix[xtal+8] + mymatrix[xtal+9] + mymatrix[xtal+12] + mymatrix[xtal+13] + mymatrix[xtal+14] + mymatrix[xtal+15] + mymatrix[xtal+16]; }
0051
0052   return E5x5;
0053 }
0054
0055 // lateral energy along ieta - syjun

0056 void energy_ieta(double mymatrix[], double *energyieta)
0057 {
0058   for(int jj = 0 ; jj < 5 ; jj++) energyieta[jj] = 0.0;
0059
0060   for(int jj = 0 ; jj < 5 ; jj++){
0061     for (int ii = 0; ii < 5 ; ii++) energyieta[jj] += mymatrix[(8+jj)+7*ii];
0062   }
0063 }
0064
0065 // lateral energy along iphi - syjun

0066 void energy_iphi(double mymatrix[], double *energyiphi)
0067 {
0068   for(int jj = 0 ; jj < 5 ; jj++) energyiphi[jj] = 0.0;
0069
0070   for(int jj = 0 ; jj < 5 ; jj++){
0071     for (int ii = 0; ii < 5 ; ii++) energyiphi[jj] += mymatrix[(8+7*jj)+ii];
0072   }
0073 }
0074
0075
0076 // 7x7 matrix energy around the max amplitude xtal

0077 double ene7x7_xtal(double mymatrix[])
0078 {
0079   double E7x7 = 0.;
0080
0081   for (int ii=0;ii<49;ii++)
0082     {
0083       if (mymatrix[ii]<-50)
0084     E7x7 = -1000.;
0085       else
0086     E7x7 += mymatrix[ii];
0087     }
0088
0089   return E7x7;
0090 }
0091
0092
0093 // -------------------------------------------

0094 // fitting functions

0095 // -------------------------------------------

0096
0097 // crystal ball fit

0098 double crystalball(double *x, double *par) {
0099   // par[0]:  mean

0100   // par[1]:  sigma

0101   // par[2]:  alpha, crossover point

0102   // par[3]:  n, length of tail

0103   // par[4]:  N, normalization

0104
0105   double cb = 0.0;
0106   double exponent = 0.0;
0107   double bla = 0.0;
0108
0109   if (x[0] > par[0] - par[2]*par[1]) {
0110     exponent = (x[0] - par[0])/par[1];
0111     cb = exp(-exponent*exponent/2.);
0112   } else {
0113     double nenner  = pow(par[3]/par[2], par[3])*exp(-par[2]*par[2]/2.);
0114     double zaehler = (par[0] - x[0])/par[1] + par[3]/par[2] - par[2];
0115     zaehler = pow(zaehler, par[3]);
0116     cb = nenner/zaehler;
0117   }
0118
0119   if (par[4] > 0.) {
0120     cb *= par[4];
0121   }
0122   return cb;
0123 }
0124
0125 double shapeFunction(double *x, double *par) {
0126
0127   // par[0] = constant value

0128   // par[1], par[2], par[3]: a3, a2 x>0 pol3

0129   // par[4], par[5], par[6]: a3, a2 x<=0 pol3

0130
0131   double ret_val;
0132   if(x[0]>0.)
0133     ret_val = par[0] + par[1]*x[0]*x[0] + par[2]*x[0]*x[0]*x[0] + par[3]*x[0]*x[0]*x[0]*x[0];
0134   else
0135     ret_val = par[0] + par[4]*x[0]*x[0] + par[5]*x[0]*x[0]*x[0] + par[6]*x[0]*x[0]*x[0]*x[0];
0136
0137   return ret_val;
0138 }