KITutil.h

#ifndef KITUTIL_h
#define KITUTIL_h 1

#include "lcio.h"
#include "EVENT/LCIO.h"
#include <EVENT/SimCalorimeterHit.h>
#include <EVENT/CalorimeterHit.h>
#include <vector>
#include <list>
#include <stack>
#include <string>
#include <UTIL/CellIDDecoder.h>
#include "Phys_Geom_Database.h" 
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_pow_int.h>
#include <ANN/ANN.h>

using namespace lcio;
using namespace std;

class Tmpcl2;
class Superhit2
{
 public:

  Superhit2(const Superhit2&) = delete;
  Superhit2& operator=(const Superhit2&) = delete;
  Superhit2(float E,CalorimeterHit* chitin);  
  ~Superhit2();
  const float* getPosition(int i);
  CalorimeterHit* chit{};
  
  bool connect{};
  float point[3]{};  // transform position
  float pointt[3]{}; // true position
  float mip{};
  float mipE{};
  int top{};         // topological parameter
  bool is_assigned{};
  vector <Superhit2*> neighbours{};
  Tmpcl2 * cl{};
  int S{};
  int K{};
  int M{};
  int tip{};
 
};  
class Tmpcl2{

 public : 
   Tmpcl2();
  ~Tmpcl2();
  void calcCenter();  
  double* getCenter();
  double getEnergy();
  void findInertia();

  // vector<Superhit2*> konektori; 
   vector<Superhit2*> hits{};
  vector<Tmpcl2*> daughters{};
  vector<Tmpcl2*> parents{};
  double energy{};
  double center[3]{};
  double direction[3]{};
  double inteigen[3]{}; 
  double inteigenvec[9]{};     

  int type{};      
};

class Photon2
{

 public:  
  Photon2(double Ein,double* pravac, double* pocetak);
  ~Photon2();
  
  void Prob(CalorimeterHit* ch,double cut,double* out);
  
  // data- stvari koje se racunaju jednom i gotovo 
  double z1{};
  double z2{};
  double k1{};
  double k2{};
  double k3{};
  double k4{};
  double p1{};
  double p2{};
  double p3{};
  double y{};
  double eprime{};
  double Z{};
  double x0eff{};
  double sampling{};
  double Eceff{};
  double Rm{};
  double Fs{};
  double Thom{};
  double Tsam{};

  double alfahom{};
  double alfasam{};
  double betasam{};
  //  
  double Ee{};
  double dir[3]{};
  double start[3]{};
};
typedef vector<Superhit2*>  Shitvec2;
typedef vector<Tmpcl2*>     Tmpclvec2;

typedef struct{
  int level;
  double Enom;
  double a;
  double b;
  double Emin;
  double Emax;
}  CoreCalib2;

typedef struct{
  Tmpcl2* cl;
  int level;
  double X[3];
  bool active;
}PROTSEED2;

typedef struct{
  double Rcut;
  double Distcut;
  double Coscut;
  unsigned int MinHit0;
  unsigned int MinHitSplit;
}CoreCut2;


void CreateAllShits2(LCCollection* colt,CellIDDecoder<CalorimeterHit>& id,vector<Superhit2*>* calo);

void TotalPrecalc2(vector<Superhit2*>* calo,CellIDDecoder<CalorimeterHit>& id,
                  unsigned int nelem, int Ncut);
void Precalc2(vector< Superhit2* >& shvec,double r, double z, double cell, double dist,bool tip,int stove,int module,CellIDDecoder<CalorimeterHit>& id);

void GridTransform2( CalorimeterHit* clh,float& radius, float& halfz, float& cellsize,float*X,
                    bool tip,int stove,int module,CellIDDecoder<CalorimeterHit>& id);

void FindCores2(Shitvec2* secal1, vector< Tmpclvec2>& bbb , vector <PROTSEED2> * prs2,
                unsigned int N, vector<float> miipstep, CoreCut2 Ccut);
void cluster5( vector<Superhit2*>* shv, vector<Tmpcl2*>* clv);
double  giveMeEEstimate2(int nivo,double Ecore, vector<CoreCalib2> cc);
void CreateCalibrationLDC00(vector<CoreCalib2>* cc);


void LineCaloIntersectD2( double* X1, double* dir,double&d,double&zmax, double*X);
void LineCaloIntersect2(double* X1, double* X2,double&d,double&zmax,  double*X);
double LinePointDistance2( double* X1, double* X2, double* X0);
void PointOnLine3(const double* X1,const double* X2,const float* X0,double* Xline);
void PointOnLine22(const double* Xstart,const double* dir,const float* X0,double* Xline);
void ModuleNormal2(double* X1,double& zmax, double* X0);
void ClusterInCluster2(Tmpcl2* cl, vector<Tmpcl2*>& clv);
double D_cl_cl2(Tmpcl2* cl1,Tmpcl2* cl2) ;
inline double Dot2(double* X1,double* X2);
void ClusterInCluster2(Tmpcl2* cl, vector<Tmpcl2*>& clv,vector<Tmpcl2*>& clout);

#endif