ClusterShapes.cc

  gsl_matrix* covar = gsl_matrix_alloc(npar,npar);   // covariance matrix

  data DataSet;
  DataSet.n = _nHits;
  DataSet.x = &t[0];
  DataSet.y = &s[0];
  DataSet.z = &E[0];  // _aHit[0]; // ???? normalise per volume ????


  fitfunct.f = &ShapeFitFunct;
  fitfunct.df = &dShapeFitFunct;
  fitfunct.fdf = &fdfShapeFitFunct;
  fitfunct.n = _nHits;
  fitfunct.p = npar;
  fitfunct.params = &DataSet;

  gsl_vector_view pinit = gsl_vector_view_array(par_init,npar);

  //std::cout << Solver << " " << &fitfunct << " " << &pinit.vector << std::endl;
  gsl_multifit_fdfsolver_set(Solver,&fitfunct,&pinit.vector);

  gsl_set_error_handler_off();

  // perform fit
  do {
    iter++;
    //std::cout << "Multidimensional Fit Iteration started ... ... ";
    status = gsl_multifit_fdfsolver_iterate(Solver);
    //std::cout << "--- DONE ---" << std::endl;

    if (status) break;
    status = gsl_multifit_test_delta (Solver->dx,Solver->x,abs_error,rel_error);

    // debug
    /*
    //    E0  = (float)gsl_vector_get(Solver->x,0);
    par[0] = (float)gsl_vector_get(Solver->x,0);
    par[1] = (float)gsl_vector_get(Solver->x,1);
    par[2] = (float)gsl_vector_get(Solver->x,2);
    par[3] = (float)gsl_vector_get(Solver->x,3);

    std::cout << "Status of multidimensional fit : " << status << "  "
    << "Iterations : " << iter << std::endl;
    std::cout << "E0 : " <<  "FIXED" << "\t" << "A : " << par[0] << "\t" << "B : " 
    <<  par[1] << "\t" << "D : " << par[2] << "\t" << "t0 : "  << par[3] 
    << std::endl << std::endl;
    */

  } while ( status==GSL_CONTINUE && iter < max_iter);

  //fg: jacobian has been dropped from gsl_multifit_fdfsolver in gsl 2:
  gsl_matrix * J = gsl_matrix_alloc(Solver->fdf->n, Solver->fdf->p);
  gsl_multifit_fdfsolver_jac(Solver, J);
  gsl_multifit_covar( J, rel_error, covar );
  //  gsl_multifit_covar (Solver->J,rel_error,covar);

  //  E0  = (float)gsl_vector_get(Solver->x,0);
  par[0] = (float)gsl_vector_get(Solver->x,0); // A
  par[1] = (float)gsl_vector_get(Solver->x,1); // B
  par[2] = (float)gsl_vector_get(Solver->x,2); // D
  par[3] = (float)gsl_vector_get(Solver->x,3); // t0

  gsl_multifit_fdfsolver_free(Solver);
  gsl_matrix_free(covar);
  
  chi2 = calculateChi2Fit3DProfileAdvanced(E0,par[0],par[1],par[2],par[3]);
  if (status) chi2 = -1.0;

  int result = 0;  // no error handling at the moment
  return result;

}

//=============================================================================

float ClusterShapes::getEmax(float* xStart, int& index_xStart, float* X0, float* Rm){

  if (_ifNotEigensystem == 1){
    transformToEigensystem(xStart,index_xStart,X0,Rm);
  }

  float E_max=0.0;
  //int i_max=0;
  for (int i = 0; i < _nHits; ++i) {
    if (E_max < _aHit[i]) {
      E_max = _aHit[i]; 
      //i_max = i;
    }
  }
  
  return E_max;
}

float ClusterShapes::getsmax(float* xStart, int& index_xStart, float* X0, float* Rm){

  if (_ifNotEigensystem == 1){
    transformToEigensystem(xStart,index_xStart,X0,Rm);
  }

  float E_max=0.0,xl_max=0.0;
  float xl_start=1.0e+50;
  //int i_max=0;
  for (int i = 0; i < _nHits; ++i) {
    //check the position of maximum energy deposit
    if (E_max < _aHit[i]) {
      E_max = _aHit[i]; 
      xl_max = _xl[i];
      //i_max = i;
    }
    //check the position of shower start
    if (xl_start > _xl[i]) {
      xl_start = _xl[i];
    }
  }
  
  return fabs(xl_max-xl_start);
}

float ClusterShapes::getxl20(float* xStart, int& index_xStart, float* X0, float* Rm){
  
  if (_ifNotEigensystem == 1){
    transformToEigensystem(xStart,index_xStart,X0,Rm);
  }
  
  float *xl_res = new float[_nHits];
  float *E_res = new float[_nHits];
  float E_tot=0.0;
  for (int i = 0; i < _nHits; ++i) {
    E_res[i]=_aHit[i];
    xl_res[i]=_xl[i];

    E_tot+=_aHit[i];
  }
  
  //sort deposite energy to xs ascending order
  float tmpe=0.0, tmpxl=0.0;
  for (int i = 0; i < _nHits; ++i) {
    for (int j = i + 1; j < _nHits; ++j) {
      if(xl_res[i]>xl_res[j]){
	tmpe=E_res[i];
	E_res[i]=E_res[j];
	E_res[j]=tmpe;

	tmpxl=xl_res[i];
	xl_res[i]=xl_res[j];
	xl_res[j]=tmpxl;
      }
    }
  }

  // std::cout << "check xt: " << xt_res[0] << " " << xt_res[1] << " " << xt_res[2]
  //           <<std::endl;
  float E20=0.0,okxl=0.0;
  int k=0;
  while(E20/E_tot<0.2){
    E20+=E_res[k];
    k++;
  }
  
  //final hit is located in outer radius
  okxl=xl_res[k-2];

  delete[] xl_res;
  delete[] E_res;

  return okxl;
}

float ClusterShapes::getxt90(float* xStart, int& index_xStart, float* X0, float* Rm){
  
  if (_ifNotEigensystem == 1){
    transformToEigensystem(xStart,index_xStart,X0,Rm);
  }
  
  float *xt_res = new float[_nHits];
  float *E_res = new float[_nHits];
  float E_tot=0.0;
  for (int i = 0; i < _nHits; ++i) {
    E_res[i]=_aHit[i];
    xt_res[i]=_xt[i];

    E_tot+=_aHit[i];
  }
  
  //sort deposite energy to xt ascending order
  float tmpe=0.0, tmpxt=0.0;
  for (int i = 0; i < _nHits; ++i) {
    for (int j = i + 1; j < _nHits; ++j) {
      if(xt_res[i]>xt_res[j]){
	tmpe=E_res[i];
	E_res[i]=E_res[j];
	E_res[j]=tmpe;

	tmpxt=xt_res[i];
	xt_res[i]=xt_res[j];
	xt_res[j]=tmpxt;
      }
    }
  }

  // std::cout << "check xt: " << xt_res[0] << " " << xt_res[1] << " " << xt_res[2]
  //           <<std::endl;
  float E90=0.0,okxt=0.0;
  int k=0;
  while(E90/E_tot<0.9){
    E90+=E_res[k];
    k++;
  }
  
  //final hit is located in outer radius
  okxt=xt_res[k-2];

  delete[] xt_res;
  delete[] E_res;

  return okxt;
}

//for test
void ClusterShapes::gethits(float* xStart, int& index_xStart, float* X0, float* Rm, float *okxl, float *okxt, float *oke){
  
  if (_ifNotEigensystem == 1){
    transformToEigensystem(xStart,index_xStart,X0,Rm);
  }
  
  for (int i = 0; i < _nHits; ++i) {
    okxl[i]=_xl[i];
    okxt[i]=_xt[i];
    oke[i]=_aHit[i];
  }
  
  return;
}


float ClusterShapes::getRhitMean(float* xStart, int& index_xStart, float* X0, float* Rm){

  if (_ifNotEigensystem == 1){
    transformToEigensystem(xStart,index_xStart,X0,Rm);
  }

  float MainCentre[3];

  MainCentre[0] = _analogGravity[0];
  MainCentre[1] = _analogGravity[1];
  MainCentre[2] = _analogGravity[2];

  // float Rhit[_nHits]={0};
  float Rhit=0;

  float Rhitsum=0;
  float Rhitmean=0;

  for (int i = 0; i < _nHits; ++i) {
    Rhit = sqrt(pow((_xHit[i]-MainCentre[0]),2) + pow((_yHit[i]-MainCentre[1]),2));
    Rhitsum += Rhit;
  }
  // cogx-- xx[0]  
  
  Rhitmean = Rhitsum/_nHits;

  return Rhitmean;
}

float ClusterShapes::getRhitRMS(float* xStart, int& index_xStart, float* X0, float* Rm){

  if (_ifNotEigensystem == 1){
    transformToEigensystem(xStart,index_xStart,X0,Rm);
  }

  float MainCentre[3];

  MainCentre[0] = _analogGravity[0];
  MainCentre[1] = _analogGravity[1];
  MainCentre[2] = _analogGravity[2];

  // float Rhit[_nHits]={0};
  float Rhit=0;

  float Rhit2sum=0;
  float Rhitrms=0;

  for (int i = 0; i < _nHits; ++i) {
    Rhit = sqrt(pow((_xHit[i]-MainCentre[0]),2) + pow((_yHit[i]-MainCentre[1]),2));
    Rhit2sum += pow(Rhit,2);
  }
  // cogx-- xx[0]  
  
  Rhitrms = sqrt(Rhit2sum/_nHits);

  return Rhitrms;
}