1 #ifndef KITUTIL_LDC00_h
   2 #define KITUTIL_LDC00_h 1
   3 
   4 #include "lcio.h"
   5 #include "EVENT/LCIO.h"
   6 #include <EVENT/SimCalorimeterHit.h>
   7 #include <EVENT/CalorimeterHit.h>
   8 #include <vector>
   9 #include <list>
  10 #include <stack>
  11 #include <string>
  12 #include <UTIL/CellIDDecoder.h>
  13 #include "myPhys_Geom_Database.h" 
  14 #include <gsl/gsl_vector.h>
  15 #include <gsl/gsl_matrix.h>
  16 #include <gsl/gsl_blas.h>
  17 #include <gsl/gsl_linalg.h>
  18 #include <gsl/gsl_eigen.h>
  19 #include <gsl/gsl_multifit_nlin.h>
  20 #include <gsl/gsl_rng.h>
  21 #include <gsl/gsl_sf_gamma.h>
  22 #include <gsl/gsl_integration.h>
  23 #include <gsl/gsl_sf_pow_int.h>
  24 #include <ANN/ANN.h>
  25 #include "mySuperhit.h"
  26 
  27 using namespace lcio;
  28 using namespace std;
  29 
  30 class Tmpcl2;
  31 /**
  32  *    Basic hit class for reconstruction, contains the calorimeter hit plus     
  33  *    additional parameters
  34  *    @authors P.Krstonosic (DESY)
  35  */
  36 class Superhit2
  37 {
  38  public:
  39 
  40   /**
  41    *  Constructor , needs absolute MIP energy (calibrated ) , and pointer to 
  42    *  calorimeter hit
  43    */
  44   Superhit2(float E,CalorimeterHit* chit_in);  
  45   /**
  46    *    Destructor
  47    */
  48   ~Superhit2();
  49   /**
  50    *   Returns position of the hit i=1 transformed position, i=2 true position
  51    */
  52   const float* getPosition(int i);
  53   /**
  54    *   Pointer to the LCIO calorimeter hit
  55    */
  56   CalorimeterHit* chit;
  57   
  58   bool connect;
  59   /**
  60    *   Transformed position of the hit
  61    */
  62   float point[3];  // transform position
  63   /**
  64    *   Real coordinate of the hit, copy of calorimeter hit data for direct access
  65    */
  66   float pointt[3]; // true position
  67   /**
  68    *   Energy of hit i terms of MIP
  69    */ 
  70   float mip;
  71   /**
  72    *   MIP value [GeV]
  73    */
  74   float mipE;
  75   /**
  76    *   Number of neighbours
  77    */
  78   int top;         // topological parameter
  79   /**
  80    *    Is hit assigned to a cluster or not 
  81    */
  82   bool is_assigned;  
  83   /**
  84    *    Vector of pointers to the neighbouring hits of Superhit2 type
  85    */
  86   vector <Superhit2*> neighbours;
  87   /**
  88    *   Pointer to the cluster to wich hit belong 
  89    */
  90   Tmpcl2 * cl;
  91   /**
  92    *  stave as coded by Mokka 
  93    */
  94   int S;
  95   /**
  96    *  layer as coded by Mokka 
  97    */
  98   int K;
  99   /**
 100    *  module as coded by Mokka
 101    */
 102   int M;
 103   /**
 104    *  0 by constructor  1 for ecal 2 for hcal 
 105    */ 
 106   int isCalorimeter; 
 107  
 108 };  
 109 /**
 110  *    Basic cluster class for reconstruction
 111  *    @authors P.Krstonosic (DESY)
 112  */
 113 class Tmpcl2{
 114 
 115  public : 
 116   /**
 117     * Constructor.
 118     */
 119    Tmpcl2();
 120   /**
 121     * Destructor. 
 122     */
 123   ~Tmpcl2();
 124   /**
 125    *  Calculates center of inerta for the objects in hits vector
 126    */
 127   void calcCenter();  
 128   /**
 129    *  returnes double[3] for the center as calculated with calcCenter method
 130    */ 
 131   double* getCenter();
 132   /**
 133    *  returnes energy of cluster in GeV.
 134    */
 135   double getEnergy();
 136   /**
 137    *  calculates eigenvalues and eigenvectors of inertia tensor
 138    */
 139   void findInertia();
 140 
 141   // vector<Superhit2*> konektori; 
 142   /**
 143    * hit in the cluster
 144    */
 145    vector<Superhit2*> hits;
 146    /**
 147     * pointers to clusters that are contained in this cluster
 148     */
 149   vector<Tmpcl2*> daughters;
 150   /**
 151    * pointers to clusters that contain this cluster
 152    */
 153   vector<Tmpcl2*> parents;
 154   /**
 155    * energy  in GeV (sum over hits in cluster) 
 156    */
 157   double energy;  
 158   /**
 159    *  position of cluster center (x,y,z)
 160    */
 161   double center[3];
 162   /**
 163    *  principal axes of inertia tensor , calculated in calcInertia
 164    */
 165   double direction[3];
 166   /**
 167    * normalized eigenvalues of inertia tensor
 168    */
 169   double inteigen[3]; 
 170   /**
 171    * eigenvectors of inertia tensor
 172    */       
 173   double inteigenvec[9];     
 174 
 175   /**
 176    *  Internal type of the cluster
 177    */
 178   int type ;      
 179 };
 180 
 181 class Photon2
 182 {
 183 
 184  public:  
 185   /**
 186    * Constructor.
 187    */
 188   Photon2(double Ein,double* pravac, double* pocetak);
 189    /**
 190     * Destructor.
 191     */
 192   ~Photon2();
 193   
 194   void Prob(CalorimeterHit* ch,double cut,double* out);
 195   
 196   // data- stvari koje se racunaju jednom i gotovo 
 197   double z1;
 198   double z2;
 199   double k1;
 200   double k2;
 201   double k3;
 202   double k4;
 203   double p1;
 204   double p2;
 205   double p3;
 206   double y;
 207   double eprime;
 208   double Z;
 209   double x0eff;
 210   double sampling;
 211   double Eceff;
 212   double Rm;
 213   double Fs;
 214   double Thom;
 215   double Tsam;
 216 
 217   double alfahom;
 218   double alfasam; 
 219   double betasam;
 220   //  
 221   double Ee;
 222   double dir[3];
 223   double start[3];
 224 };
 225 typedef vector<Superhit2*>  Shitvec2;
 226 typedef vector<Tmpcl2*>     Tmpclvec2;
 227 
 228 /**
 229  * container for holding the numbers needed for energy estimation
 230  */
 231 typedef struct{
 232   /**
 233    * level of cluster
 234    */
 235   int level;
 236   /**
 237    * nominal i.e. input energy in GeV
 238    */
 239   double Enom;
 240   /**
 241    *   Eestimate = a+b*Ecore ,  a parameter of this funcition
 242    */
 243   double a;
 244   /**
 245    *   Eestimate = a+b*Ecore ,  b parameter of this funcition
 246    */
 247   double b;
 248   /**
 249    *   lower validity range of energy estimation funciton in GeV
 250    */
 251   double Emin;
 252  /**
 253    *   upper validity range of energy estimation funciton in GeV
 254    */
 255   double Emax;
 256 }  CoreCalib2;
 257 
 258 /**
 259  * container for storing the EM shower core candidates
 260  */
 261 typedef struct{
 262   /**
 263    * pointer to the cluster
 264    */
 265   Tmpcl2* cl;
 266   /**
 267    *  level of the cluster
 268    */
 269   int level;
 270   /**
 271    *  center of the candidate
 272    */
 273   double X[3];
 274   /**
 275    *  flag to activate deactivate 
 276    */
 277   bool active;
 278 }PROTSEED2;
 279 
 280 /**
 281  *  container for keeping the parameters of the core fineder together
 282  */
 283 typedef struct{
 284   /**
 285    * fluctuation suppresion cut 
 286    */
 287   double rCut;
 288   /**
 289    *  distance cut for core merging 
 290    */
 291   double distCut;
 292   /**
 293    *  angular cut for core merging (value of the cosine is stored not the angle!)
 294    */
 295   double cosCut;
 296   /**
 297    * minimal number of hits needed for 0-th level cluster to be accepted as a core candidate
 298    */
 299   unsigned int MinHit0;
 300   /**
 301    * minimal number of hits in i-th level cluster to be accepted for splitting  of the core
 302    */
 303   unsigned int MinHitSplit;
 304 }CoreCut;
 305 
 306 
 307 /**
 308  *  Creation of superhits, input is ECAL collection ,it's decoded and pointer to resulting container for
 309  *  superhits
 310  */
 311 void CreateAllShits2(LCCollection* colt,CellIDDecoder<CalorimeterHit>& id,vector<Superhit2*>* calo, int debug);
 312 
 313 /**
 314  *  Global precalculation function , iput is vector of superhits, ECAL decoder, number of hits, and 
 315  *  number of neighbors  cut for hit separation
 316  */
 317 void TotalPrecalc2(vector<Superhit2*>* calo,CellIDDecoder<CalorimeterHit>& id,unsigned int nelem, int Ncut, int debug);
 318 /**
 319  *  Precalculation function used internaly by TotalPrecalc 
 320  */
 321 void Precalc2(vector< Superhit2* >& shvec,double r, double z, double cell, double dist,bool isCalorimeter,int stave,int module,CellIDDecoder<CalorimeterHit>& id, int debug);
 322 
 323 /**
 324  *  Basic function for transformation of hit coordinates
 325  */
 326 void GridTransform2( CalorimeterHit* clh,float& radius, float& halfz, float& cellsize,float*X,bool isCalorimeter,int stave,int module,CellIDDecoder<CalorimeterHit>& id);
 327 
 328 /**
 329  *  Global EM core finding function , iput is vector of superhits, array of Tmpcl vectors bbb for internal
 330  *  computation, vector of EM core candidates prs2, and parameters of the algorithm
 331  */
 332 void FindCores2(Shitvec2* secal1, Tmpclvec2* bbb , vector <PROTSEED2> * prs2,
 333 		unsigned int N, vector<float> mipstep, CoreCut Ccut, int debug);
 334 /**
 335  *  NN clustering 
 336  */
 337 void cluster5( vector<Superhit2*>* shv, vector<Tmpcl2*>* clv, int debug);
 338 /**
 339  * returns energy estimate for a given core , input  int level, double core energy in GeV, and 
 340  * calibration data table
 341  */
 342 double giveMeEEstimate2(int nivo,double Ecore, vector<CoreCalib2> cc);
 343 /**
 344  * Example function for creation of the energy estimaiton table 
 345  */
 346 void CreateCalibration_LDC00(vector<CoreCalib2>* cc);
 347 void CreateCalibration_LDCPrime02Sc(vector<CoreCalib2>* cc);
 348 
 349 void LineCaloIntersectD2( double* X1, double* dir,double&d,double&zmax, double*X, int debug);
 350 void LineCaloIntersect2(double* X1, double* X2,double&d,double&zmax,  double*X, int debug);
 351 double LinePointDistance2( double* X1, double* X2, double* X0);
 352 void PointOnLine3(const double* X1,const double* X2,const float* X0,double* Xline);
 353 void PointOnLine22(const double* Xstart,const double* dir,const float* X0,double* Xline);
 354 void ModuleNormal2(double* X1,double& zmax, double* X0);
 355 void ClusterInCluster2(Tmpcl2* cl, vector<Tmpcl2*>& clv, int debug);
 356 double D_cl_cl2(Tmpcl2* cl1,Tmpcl2* cl2) ;
 357 inline double Dot2(double* X1,double* X2);
 358 void ClusterInCluster2(Tmpcl2* cl, vector<Tmpcl2*>& clv,vector<Tmpcl2*>& clout, int debug);
 359 
 360 #endif
 361