ClusterShapes.h

#ifndef ClusterShapes_h
#define ClusterShapes_h


#include <stdio.h>
#include <iostream>
#include <iomanip>
#include <string>
#include <sstream>
#include <cstdlib>
#include "HelixClass.h"
#include <math.h>


class ClusterShapes {

public:

  ClusterShapes(int nhits, float* a, float* x, float* y, float* z);


  ~ClusterShapes();

  void setErrors(float *ex, float* ey, float *ez);

  void setHitTypes(int *ityp);


  int getNumberOfHits();

  float getTotalAmplitude();

  float* getCentreOfGravity();
  // this is (for now) a pure dummy to allow MarlinPandora development!
  float* getCentreOfGravityErrors();

  inline float* getCenterOfGravity() { return getCentreOfGravity() ; }
  // this is (for now) a pure dummy to allow MarlinPandora development!
  inline float* getCenterOfGravityErrors() { return getCentreOfGravityErrors() ; }
  
  float* getEigenValInertia();
  // this is (for now) a pure dummy to allow MarlinPandora development!
  float* getEigenValInertiaErrors();

  float* getEigenVecInertia();
  // this is (for now) a pure dummy to allow MarlinPandora development!
  float* getEigenVecInertiaErrors();

  float getWidth();

  int getEigenSytemCoordinates(float* xlong, float* xtrans);

  int getEigenSytemCoordinates(float* xlong, float* xtrans, float* a);

  int fit3DProfile(float& chi2, float& a, float& b, float& c, float& d, float& xl0, 
                   float * xStart, int& index_xStart, float* X0, float* Rm);

  float getChi2Fit3DProfileSimple(float a, float b, float c, float d, float* X0, 
                                  float* Rm);

  float getChi2Fit3DProfileAdvanced(float E0, float a, float b, float d, float t0, 
                                    float* X0, float* Rm);

  int FitHelix(int max_iter, int status_out, int parametrisation,
               double* parameter, double* dparameter, double& chi2, double& distmax, int direction=1);


  int FitHelix(int max_iter, int status_out, int parametrisation,
               float* parameter, float* dparameter, float& chi2, float& distmax, int direction=1);

  //here add my functions(variables estimated with detector base)
  //maximum deposit energy of hits
  float getEmax(float* xStart, int& index_xStart, float* X0, float* Rm);

  //shower max of the hits from the shower start hit
  float getsmax(float* xStart, int& index_xStart, float* X0, float* Rm);

  //radius where 90% of the cluster energy exists
  float getxt90(float* xStart, int& index_xStart, float* X0, float* Rm);

  //length where less than 20% of the cluster energy exists
  float getxl20(float* xStart, int& index_xStart, float* X0, float* Rm);

  //for cluster study
  void gethits(float* xStart, int& index_xStart, float* X0, float* Rm, float *okxl, float *okxt, float *oke);
  inline float radius() { return _radius; }

  inline float getElipsoid_r1() { return _r1; }

  inline float getElipsoid_r2() { return _r2; }

  inline float getElipsoid_r3() { return _r3; }

  inline float getElipsoid_vol() { return _vol; }

  inline float getElipsoid_r_ave() { return _r_ave; }

  inline float getElipsoid_density() { return _density; }

  inline float getElipsoid_eccentricity() { return _eccentricity; }

  inline float getElipsoid_r_forw() { return _r1_forw; }

  inline float getElipsoid_r_back() { return _r1_back; }

  //Mean of the radius of the hits
  float getRhitMean(float* xStart, int& index_xStart, float* X0, float* Rm);

  //RMS of the radius of the hits
  float getRhitRMS(float* xStart, int& index_xStart, float* X0, float* Rm);



private:

  int _nHits;

  float* _aHit;
  float* _xHit;
  float* _yHit;
  float* _zHit;
  float* _exHit;
  float* _eyHit;
  float* _ezHit;
  int* _types;
  float* _xl;
  float* _xt;
  float* _t;
  float* _s;

  int   _ifNotGravity;
  float _totAmpl;
  float _radius;
  float _xgr;
  float _ygr;
  float _zgr;
  float _analogGravity[3];

  int   _ifNotWidth;
  float _analogWidth;

  int   _ifNotInertia;
  float _ValAnalogInertia[3];
  float _VecAnalogInertia[9];

  int _ifNotEigensystem;

  int   _ifNotElipsoid;
  float _r1           ;  // Cluster spatial axis length -- the largest
  float _r2           ;  // Cluster spatial axis length -- less
  float _r3           ;  // Cluster spatial axis length -- less
  float _vol          ;  // Cluster ellipsoid volume
  float _r_ave        ;  // Cluster average radius  (qubic root)
  float _density      ;  // Cluster density
  float _eccentricity ;  // Cluster Eccentricity
  float _r1_forw      ;
  float _r1_back      ;

  void  findElipsoid();
  void  findGravity();
  void  findInertia();
  void  findWidth();
  float findDistance(int i);
  float vecProduct(float * x1, float * x2);
  float vecProject(float * x, float * axis);
  double DistanceHelix(double x, double y, double z, double X0, double Y0, double R0, double bz,
                       double phi0, double * distRPhiZ);
  int transformToEigensystem(float* xStart, int& index_xStart, float* X0, float* Xm);
  float calculateChi2Fit3DProfileSimple(float a, float b, float c, float d);
  float calculateChi2Fit3DProfileAdvanced(float E0, float a, float b, float d, float t0);
  int fit3DProfileSimple(float& chi2, float& a, float& b, float& c, float& d);
  int fit3DProfileAdvanced(float& chi2, double* par_init, double* par, int npar,
                           float* t, float* s, float* E, float E0);

  // private methods for non-linear, multidim. fitting (helix)
  // static int functParametrisation1(const gsl_vector* par, void* data, gsl_vector* f);
  // static int dfunctParametrisation1(const gsl_vector* par, void* d, gsl_matrix* J);
  // static int fdfParametrisation1(const gsl_vector* par, void* d, gsl_vector* f, gsl_matrix* J);


};


#endif