diff --git a/Utilities/DataHelper/CMakeLists.txt b/Utilities/DataHelper/CMakeLists.txt
index 3d1e61224c731b39cc3527efa526078ba36a7e9c..a0b346f3580fbbf1567c01dda9710193118d43a0 100644
--- a/Utilities/DataHelper/CMakeLists.txt
+++ b/Utilities/DataHelper/CMakeLists.txt
@@ -1,10 +1,11 @@
gaudi_subdir(DataHelper v0r0)
find_package(EDM4HEP REQUIRED)
-set (gsl_include "/cvmfs/cepc.ihep.ac.cn/software/cepcsoft/x86_64-sl6-gcc49/external/GSL/1.14/install/include")
-set (gsl_lib1 "/cvmfs/cepc.ihep.ac.cn/software/cepcsoft/x86_64-sl6-gcc49/external/GSL/1.14/install/lib/libgsl.so")
-set (gsl_lib2 "/cvmfs/cepc.ihep.ac.cn/software/cepcsoft/x86_64-sl6-gcc49/external/GSL/1.14/install/lib/libgslcblas.so")
-INCLUDE_DIRECTORIES(${gsl_include})
+find_package(GSL REQUIRED )
+message("GSL: ${GSL_LIBRARIES} ")
+message("GSL INCLUDE_DIRS: ${GSL_INCLUDE_DIRS} ")
+
+INCLUDE_DIRECTORIES(${GSL_INCLUDE_DIRS})
gaudi_depends_on_subdirs()
@@ -14,5 +15,5 @@ set(DataHelperLib_srcs src/*.cc src/*.cpp)
gaudi_add_library(DataHelperLib ${DataHelperLib_srcs}
PUBLIC_HEADERS DataHelper
- LINK_LIBRARIES EDM4HEP::edm4hep EDM4HEP::edm4hepDict ${gsl_lib1} ${gsl_lib2}
+ LINK_LIBRARIES EDM4HEP::edm4hep EDM4HEP::edm4hepDict ${GSL_LIBRARIES}
)
diff --git a/Utilities/DataHelper/DataHelper/ClusterShapes.h b/Utilities/DataHelper/DataHelper/ClusterShapes.h
index 695dabeed2be2e97cde2217b8340e04fee8f78c6..3bcd5db00b82fdeefc21e70afe0843f376be2ec2 100644
--- a/Utilities/DataHelper/DataHelper/ClusterShapes.h
+++ b/Utilities/DataHelper/DataHelper/ClusterShapes.h
@@ -8,7 +8,7 @@
#include <string>
#include <sstream>
#include <cstdlib>
-#include "DataHelper/HelixClass.h"
+#include "HelixClass.h"
#include <math.h>
@@ -24,7 +24,7 @@
*/
class ClusterShapes {
- public:
+public:
/**
* Constructor
@@ -48,20 +48,20 @@ class ClusterShapes {
/**
* Defining errors for Helix Fit
*/
- void setErrors(float *ex, float* ey, float *ez);
+ void setErrors(float *ex, float* ey, float *ez);
- /**
- * Defining hit types for Helix Fit :
- * type 1 - cyllindrical detector
- * type 2 - Z disk detector
- */
- void setHitTypes(int *ityp);
+ /**
+ * Defining hit types for Helix Fit :
+ * type 1 - cyllindrical detector
+ * type 2 - Z disk detector
+ */
+ void setHitTypes(int *ityp);
/**
* returns the number of elements of the cluster
*/
- int getNumberOfHits();
+ int getNumberOfHits();
/**
* returns the accumulated amplitude for the whole cluster (just the sum of the
@@ -75,15 +75,21 @@ class ClusterShapes {
* of gravity is calculated with the energy of the entries of the cluster.
*/
float* getCentreOfGravity();
+ // this is (for now) a pure dummy to allow MarlinPandora development!
+ float* getCentreOfGravityErrors();
/** US spelling of getCentreOfGravity */
inline float* getCenterOfGravity() { return getCentreOfGravity() ; }
+ // this is (for now) a pure dummy to allow MarlinPandora development!
+ inline float* getCenterOfGravityErrors() { return getCentreOfGravityErrors() ; }
/**
* array of the inertias of mass (i.\ e.\ energy) corresponding to the three main axes
* of inertia. The array is sorted in ascending order.
*/
float* getEigenValInertia();
+ // this is (for now) a pure dummy to allow MarlinPandora development!
+ float* getEigenValInertiaErrors();
/**
* array of the three main axes of inertia (9 entries) starting
@@ -92,6 +98,8 @@ class ClusterShapes {
* of 1.
*/
float* getEigenVecInertia();
+ // this is (for now) a pure dummy to allow MarlinPandora development!
+ float* getEigenVecInertiaErrors();
/**
* 'mean' width of the cluster perpendicular to the main
@@ -144,7 +152,7 @@ class ClusterShapes {
* detector this is meant to be the 'mean' Moliere radius.
*/
int fit3DProfile(float& chi2, float& a, float& b, float& c, float& d, float& xl0,
- float * xStart, int& index_xStart, float X0, float Rm);
+ float * xStart, int& index_xStart, float* X0, float* Rm);
/**
* returns the chi2 of the fit in the method Fit3DProfile (if simple
@@ -155,8 +163,8 @@ class ClusterShapes {
* @param Rm : Moliere radius of the the detector material. For a composite
* detector this is meant to be the 'mean' Moliere radius.
*/
- float getChi2Fit3DProfileSimple(float a, float b, float c, float d, float X0,
- float Rm);
+ float getChi2Fit3DProfileSimple(float a, float b, float c, float d, float* X0,
+ float* Rm);
/**
* returns the chi2 of the fit in the method Fit3DProfile (if advanced
@@ -168,7 +176,7 @@ class ClusterShapes {
* detector this is meant to be the 'mean' Moliere radius.
*/
float getChi2Fit3DProfileAdvanced(float E0, float a, float b, float d, float t0,
- float X0, float Rm);
+ float* X0, float* Rm);
/**
* performs a least square fit on a helix path in space, which
@@ -220,6 +228,21 @@ class ClusterShapes {
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);
/**
* distance to the centre of gravity measured from IP
* (absolut value of the vector to the centre of gravity)
@@ -280,54 +303,58 @@ class ClusterShapes {
*/
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:
+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;
+ std::vector<float> _aHit;
+ std::vector<float> _xHit;
+ std::vector<float> _yHit;
+ std::vector<float> _zHit;
+ std::vector<float> _exHit;
+ std::vector<float> _eyHit;
+ std::vector<float> _ezHit;
+ std::vector<float> _xl;
+ std::vector<float> _xt;
+ std::vector<float> _t;
+ std::vector<float> _s;
+ std::vector<int> _types;
+
+ int _ifNotGravity=1;
+ float _totAmpl=0.0;
+ float _radius=0.0;
+ float _xgr=0.0;
+ float _ygr=0.0;
+ float _zgr=0.0;
+ float _analogGravity[3]={0.0, 0.0,0.0};
+
+ int _ifNotWidth=1;
+ float _analogWidth=0.0;
+
+ int _ifNotInertia=1;
float _ValAnalogInertia[3];
float _VecAnalogInertia[9];
- int _ifNotEigensystem;
+ int _ifNotEigensystem=1;
- 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 ;
+ //int _ifNotElipsoid=1;
+ float _r1 =0.0; // Cluster spatial axis length -- the largest
+ float _r2 =0.0; // Cluster spatial axis length -- less
+ float _r3 =0.0; // Cluster spatial axis length -- less
+ float _vol =0.0; // Cluster ellipsoid volume
+ float _r_ave =0.0; // Cluster average radius (qubic root)
+ float _density =0.0; // Cluster density
+ float _eccentricity =0.0; // Cluster Eccentricity
+ float _r1_forw =0.0;
+ float _r1_back =0.0;
void findElipsoid();
void findGravity();
@@ -337,8 +364,8 @@ class ClusterShapes {
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);
+ 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);
@@ -353,4 +380,5 @@ class ClusterShapes {
};
+
#endif
diff --git a/Utilities/DataHelper/src/ClusterShapes.cc b/Utilities/DataHelper/src/ClusterShapes.cc
index e86bcefac13f546dbea1292567c2e841c3b62cad..a5fbba971fde00a57acd352bc325d9cd3f5d2d96 100644
--- a/Utilities/DataHelper/src/ClusterShapes.cc
+++ b/Utilities/DataHelper/src/ClusterShapes.cc
@@ -1,3 +1,5 @@
+
+
/* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
#include "DataHelper/ClusterShapes.h"
@@ -14,7 +16,6 @@
//#include <gsl/gsl_sf_pow_int.h>
-
// #################################################
// ##### #####
// ##### Additional Structures and Functions #####
@@ -46,7 +47,7 @@ double G(double x) {
// inverse Gammafunction
double invG(double x) {
- return gsl_sf_gammainv(x);
+ return gsl_sf_gammainv(x);
}
@@ -68,24 +69,23 @@ double DinvG(double x) {
double rel_error = 1e-6;
double result = 0.0;
double error = 0.0;
- int status = 0;
gsl_integration_workspace* w = gsl_integration_workspace_alloc(workspace_size);
gsl_function F;
F.function = &Integral_G;
F.params = &x;
- status = gsl_integration_qagiu(&F,0,abs_error,rel_error,workspace_size,w,
- &result,&error);
+ /*int status=*/gsl_integration_qagiu(&F,0,abs_error,rel_error,workspace_size,w,
+ &result,&error);
// debug
/*
- printf ("Numeric Integration : \n");
- printf ("parameter of integration = % .18f\n", x);
- printf ("status of integration = %d \n" , status);
- printf ("result = % .18f\n", result);
- printf ("estimated error = % .18f\n", error);
- printf ("intervals = %d\n\n", w->size);
+ printf ("Numeric Integration : \n");
+ printf ("parameter of integration = % .18f\n", x);
+ printf ("status of integration = %d \n" , status);
+ printf ("result = % .18f\n", result);
+ printf ("estimated error = % .18f\n", error);
+ printf ("intervals = %d\n\n", w->size);
*/
double G2 = pow(gsl_sf_gamma(x),2);
@@ -128,7 +128,7 @@ int ShapeFitFunct(const gsl_vector* par, void* d, gsl_vector* f) {
for (int i(0); i < n; i++) {
fi = /*E0 * */ B * invG(A) * pow(B*(t[i]-t0),A-1) * exp(-B*(t[i]-t0)) * exp(-D*s[i])
- - a[i];
+ - a[i];
gsl_vector_set(f,i,fi);
}
@@ -155,30 +155,30 @@ int dShapeFitFunct(const gsl_vector* par, void* d, gsl_matrix* J) {
for (int i(0); i < n; i++) {
/*
- gsl_matrix_set(J,i,0,B * invG(A) * pow(B*(t[i]-t0),A-1) * exp(-B*(t[i]-t0))
- * exp(-D*s[i]) );
- */
+ gsl_matrix_set(J,i,0,B * invG(A) * pow(B*(t[i]-t0),A-1) * exp(-B*(t[i]-t0))
+ * exp(-D*s[i]) );
+ */
gsl_matrix_set(J,i,0,( /* E0 * */ B * invG(A) * log(B*(t[i]-t0))*pow(B*(t[i]-t0),A-1) *
exp(-B*(t[i]-t0)) + DinvG(A) * /* E0 * */ B * pow(B*(t[i]-t0),A-1) *
- exp(-B*(t[i]-t0))
- ) * exp(-D*s[i]));
+ exp(-B*(t[i]-t0))
+ ) * exp(-D*s[i]));
gsl_matrix_set(J,i,1,( /* E0 * */ invG(A) * pow(B*(t[i]-t0),A-1) * exp(-B*(t[i]-t0)) +
/* E0 * */ invG(A) * (A-1) * B * (t[i]-t0) * pow(B*(t[i]-t0),A-2) *
- exp(-B*(t[i]-t0)) -
+ exp(-B*(t[i]-t0)) -
/* E0 * */ B * invG(A) * (t[i]-t0) * pow(B*(t[i]-t0),A-1) *
- exp(-B*(t[i]-t0))
- ) * exp(-D*s[i]));
+ exp(-B*(t[i]-t0))
+ ) * exp(-D*s[i]));
gsl_matrix_set(J,i,2,-/* E0 * */ B * invG(A) * s[i] * pow(B*(t[i]-t0),A-1) *
- exp(-B*(t[i]-t0)) * exp(-D*s[i]));
+ exp(-B*(t[i]-t0)) * exp(-D*s[i]));
gsl_matrix_set(J,i,3,(-/* E0 * */ pow(B,2) * invG(A) * (A-1) * pow(B*(t[i]-t0),A-2) *
- exp(-B*(t[i]-t0)) +
- /* E0 * */ pow(B,2) * invG(A) * pow(B*(t[i]-t0),A-1) *
- exp(-B*(t[i]-t0))
- ) * exp(-D*s[i]));
+ exp(-B*(t[i]-t0)) +
+ /* E0 * */ pow(B,2) * invG(A) * pow(B*(t[i]-t0),A-1) *
+ exp(-B*(t[i]-t0))
+ ) * exp(-D*s[i]));
}
@@ -229,9 +229,9 @@ int functParametrisation1(const gsl_vector* par, void* d, gsl_vector* f) {
float* x = ((struct data*)d)->x;
float* y = ((struct data*)d)->y;
float* z = ((struct data*)d)->z;
-//float* ex = ((struct data*)d)->ex;
-//float* ey = ((struct data*)d)->ey;
-//float* ez = ((struct data*)d)->ez;
+ //float* ex = ((struct data*)d)->ex;
+ //float* ey = ((struct data*)d)->ey;
+ //float* ez = ((struct data*)d)->ez;
float fi = 0.0;
// first dimension
@@ -263,9 +263,9 @@ int dfunctParametrisation1(const gsl_vector* par, void* d, gsl_matrix* J) {
int n = ((struct data*)d)->n;
float* z = ((struct data*)d)->z;
-//float* ex = ((struct data*)d)->ex;
-//float* ey = ((struct data*)d)->ey;
-//float* ez = ((struct data*)d)->ez;
+ //float* ex = ((struct data*)d)->ex;
+ //float* ey = ((struct data*)d)->ey;
+ //float* ez = ((struct data*)d)->ez;
// calculate Jacobi's matrix J[i][j] = dfi/dparj
@@ -508,7 +508,7 @@ int dfunctParametrisation3(const gsl_vector* par, void* d, gsl_matrix* J) {
gsl_matrix_set(J,i,0,0);
gsl_matrix_set(J,i,1,((1/omega) - d0) * cos(Phi0) - (1/fabs(omega)) * sin(phii) );
gsl_matrix_set(J,i,2,((-1.0)*sin(Phi0))/pow(omega,2) - ( (signum(omega))/(pow(fabs(omega),2)) ) * cos(phii) -
- (1/fabs(omega)) * sin(phii) * ( ( ((-1.0)/sqrt(1+pow(tanL,2)))*si) + (M_PI)/(2*fabs(omega)) - (signum(omega)*omega*M_PI)/(2*pow(fabs(omega),2)) ) );
+ (1/fabs(omega)) * sin(phii) * ( ( ((-1.0)/sqrt(1+pow(tanL,2)))*si) + (M_PI)/(2*fabs(omega)) - (signum(omega)*omega*M_PI)/(2*pow(fabs(omega),2)) ) );
gsl_matrix_set(J,i,3,(-1.0)*sin(Phi0));
gsl_matrix_set(J,i,4,((-1.0)/fabs(omega))*sin(phii) * ( (omega*tanL*si)/sqrt(pow(1+pow(tanL,2),3)) ) );
@@ -522,7 +522,7 @@ int dfunctParametrisation3(const gsl_vector* par, void* d, gsl_matrix* J) {
gsl_matrix_set(J,i+n,0,0);
gsl_matrix_set(J,i+n,1,((1/omega) - d0)*sin(Phi0) + (1/fabs(omega)) * cos(phii) );
gsl_matrix_set(J,i+n,2,cos(Phi0)/pow(omega,2) + ( (signum(omega))/(pow(fabs(omega),2)) ) * sin(phii) +
- (1/fabs(omega))*cos(phii) * ( ( ((-1.0)/sqrt(1+pow(tanL,2)))*si) + (M_PI)/(2*fabs(omega)) - (signum(omega)*omega*M_PI)/(2*pow(fabs(omega),2)) ) );
+ (1/fabs(omega))*cos(phii) * ( ( ((-1.0)/sqrt(1+pow(tanL,2)))*si) + (M_PI)/(2*fabs(omega)) - (signum(omega)*omega*M_PI)/(2*pow(fabs(omega),2)) ) );
gsl_matrix_set(J,i+n,3,cos(Phi0));
gsl_matrix_set(J,i+n,4,(1/fabs(omega))*cos(phii) * ( (omega*tanL*si)/sqrt(pow(1+pow(tanL,2),3)) ) );
@@ -570,39 +570,35 @@ int fdfParametrisation3(const gsl_vector* par, void* d, gsl_vector* f, gsl_matri
//=============================================================================
-ClusterShapes::ClusterShapes(int nhits, float* a, float* x, float* y, float* z){
-
- _nHits = nhits;
- _aHit = new float[_nHits];
- _xHit = new float[_nHits];
- _yHit = new float[_nHits];
- _zHit = new float[_nHits];
- _exHit = new float[_nHits];
- _eyHit = new float[_nHits];
- _ezHit = new float[_nHits];
- _types = new int[_nHits];
-
- _xl = new float[_nHits];
- _xt = new float[_nHits];
-
- _t = new float[_nHits];
- _s = new float[_nHits];
+ClusterShapes::ClusterShapes(int nhits, float* a, float* x, float* y, float* z):
+
+ _nHits(nhits),
+ _aHit (nhits, 0.0),
+ _xHit (nhits, 0.0),
+ _yHit (nhits, 0.0),
+ _zHit (nhits, 0.0),
+ _exHit (nhits, 1.0),
+ _eyHit (nhits, 1.0),
+ _ezHit (nhits, 1.0),
+ _xl (nhits, 0.0),
+ _xt (nhits, 0.0),
+ _t (nhits, 0.0),
+ _s (nhits, 0.0),
+ _types(nhits, 1), // all hits are assumed to be "cylindrical"
+
+ _ifNotGravity (1),
+ _ifNotWidth (1),
+ _ifNotInertia (1),
+ _ifNotEigensystem(1)
+{
for (int i(0); i < nhits; ++i) {
_aHit[i] = a[i];
_xHit[i] = x[i];
_yHit[i] = y[i];
_zHit[i] = z[i];
- _exHit[i] = 1.0;
- _eyHit[i] = 1.0;
- _ezHit[i] = 1.0;
- _types[i] = 1; // all hits are assumed to be "cyllindrical"
}
- _ifNotGravity = 1;
- _ifNotWidth = 1;
- _ifNotInertia = 1;
- _ifNotEigensystem = 1;
}
@@ -611,18 +607,6 @@ ClusterShapes::ClusterShapes(int nhits, float* a, float* x, float* y, float* z){
ClusterShapes::~ClusterShapes() {
- delete[] _aHit;
- delete[] _xHit;
- delete[] _yHit;
- delete[] _zHit;
- delete[] _exHit;
- delete[] _eyHit;
- delete[] _ezHit;
- delete[] _types;
- delete[] _xl;
- delete[] _xt;
- delete[] _t;
- delete[] _s;
}
//=============================================================================
@@ -672,6 +656,11 @@ float* ClusterShapes::getCentreOfGravity() {
if (_ifNotGravity == 1) findGravity() ;
return &_analogGravity[0] ;
}
+float* ClusterShapes::getCentreOfGravityErrors() {
+ // this is a pure dummy to allow MarlinPandora development!
+ if (_ifNotGravity == 1) findGravity() ;
+ return &_analogGravity[0] ;
+}
//=============================================================================
@@ -679,6 +668,11 @@ float* ClusterShapes::getEigenValInertia() {
if (_ifNotInertia == 1) findInertia();
return &_ValAnalogInertia[0] ;
}
+float* ClusterShapes::getEigenValInertiaErrors() {
+ // this is a pure dummy to allow MarlinPandora development!
+ if (_ifNotInertia == 1) findInertia();
+ return &_ValAnalogInertia[0] ;
+}
//=============================================================================
@@ -686,6 +680,11 @@ float* ClusterShapes::getEigenVecInertia() {
if (_ifNotInertia == 1) findInertia();
return &_VecAnalogInertia[0] ;
}
+float* ClusterShapes::getEigenVecInertiaErrors() {
+ // this is a pure dummy to allow MarlinPandora development!
+ if (_ifNotInertia == 1) findInertia();
+ return &_VecAnalogInertia[0] ;
+}
//=============================================================================
@@ -703,8 +702,8 @@ int ClusterShapes::getEigenSytemCoordinates(float* xlong, float* xtrans) {
// NOT SAVE, change to class variables !!!!!
- float X0 = 7.0;
- float Rm = 13.5;
+ float X0[2]={3.50,17.57}; //in mm. //this is the exact value of tungsten and iron
+ float Rm[2]={9.00,17.19}; //in mm. need to change to estimate correctly
if (_ifNotEigensystem == 1) transformToEigensystem(xStart,index_xStart,X0,Rm);
@@ -725,8 +724,8 @@ int ClusterShapes::getEigenSytemCoordinates(float* xlong, float* xtrans, float*
int index_xStart;
// NOT SAVE, change to class variables !!!!!
- float X0 = 7.0;
- float Rm = 13.5;
+ float X0[2]={3.50,17.57}; //in mm. //this is the exact value of tungsten and iron
+ float Rm[2]={9.00,17.19}; //in mm. need to change to estimate correctly
if (_ifNotEigensystem == 1) transformToEigensystem(xStart,index_xStart,X0,Rm);
@@ -744,15 +743,25 @@ int ClusterShapes::getEigenSytemCoordinates(float* xlong, float* xtrans, float*
int ClusterShapes::fit3DProfile(float& chi2, float& E0, float& A, float& B, float& D,
float& xl0, float* xStart, int& index_xStart,
- float X0, float Rm) {
+ float* X0, float* Rm) {
const int npar = 4;
-
- if (_ifNotEigensystem == 1) transformToEigensystem(xStart,index_xStart,X0,Rm);
-
+ if (_ifNotEigensystem == 1){
+ transformToEigensystem(xStart,index_xStart,X0,Rm);
+ }
+
float* E = new float[_nHits];
+ //doesn't fit when _nHits==0
+ //std::cout << "_nhits " << _nHits << std::endl;
+ if(_nHits-1 < npar){ //can't fit because number of degrees of freedom is small
+ E0=0.0;
+
+ delete[] E;
+ int result = 0; // no error handling at the moment
+ return result;
+ }
double par_init[npar];
for (int i = 0; i < npar; ++i) par_init[i] = 0.0; // initialise
@@ -760,16 +769,16 @@ int ClusterShapes::fit3DProfile(float& chi2, float& E0, float& A, float& B, floa
float E0_init = 0.0;
float A_init = 0.0;
float B_init = 0.5; // empirically
- float D_init = 0.5; // empirically
- float t0_init = -10.0/X0; // shift xl0_ini is assumed to be -10.0 without a reason ????
+ float D_init = 0.06; //0.5; //if Rm includes 90% of shower energy, absorption length has 63% of shower energy
+ float t0_init = -10.0/X0[0]; // shift xl0_ini is assumed to be -10.0 without a reason ????
float E0_tmp = 0.0;
- int i_max = 0;
- float t_max = 0.0;
+ // int i_max = 0;
+ //float t_max = 0.0;
for (int i = 0; i < _nHits; ++i) {
if (E0_tmp < _aHit[i]) {
E0_tmp = _aHit[i];
- i_max = i;
+ //i_max = i;
}
E0_init += _aHit[i];
}
@@ -783,11 +792,9 @@ int ClusterShapes::fit3DProfile(float& chi2, float& E0, float& A, float& B, floa
//A_init = t_max*B_init + 1.0;
// third definition
- float Ec = X0 * 0.021/Rm;
+ float Ec = X0[0] * 0.021/Rm[0];
A_init = B_init * log(E0_init/Ec) + 0.5 * B_init + 1.0; // (+0.5 for Photons initiated showers)
-
-
// par_init[0] = E0_init;
E0 = E0_init;
for (int i = 0; i < _nHits; ++i) E[i] = _aHit[i]/E0_init;
@@ -795,28 +802,28 @@ int ClusterShapes::fit3DProfile(float& chi2, float& E0, float& A, float& B, floa
par_init[1] = B_init;
par_init[2] = D_init;
par_init[3] = t0_init;
-
+
// debug
- std::cout << "E0_init : " << E0_init << "\t" << "A_init : " << A_init << "\t"
- << "B_init : " << B_init << "\t" << "D_init : " << D_init << "\t"
- << "xl0_init : " << t0_init*X0 << "\t" << "X0 : " << X0
- << "\t" << "t_max : " << t_max << std::endl << std::endl;
+ // std::cout << "E0_init : " << E0_init << "\t" << "A_init : " << A_init << "\t"
+ // << "B_init : " << B_init << "\t" << "D_init : " << D_init << "\t"
+ // << "xl0_init : " << t0_init*X0 << "\t" << "X0 : " << X0
+ // << "\t" << "t_max : " << t_max << std::endl << std::endl;
- double t0 = xl0/X0;
+ double t0 = xl0/X0[0]; //probably t0 is in ecal
double par[npar];
par[0] = A;
par[1] = B;
par[2] = D;
par[3] = t0;
- fit3DProfileAdvanced(chi2,par_init,par,npar,_t,_s,E,E0);
+ fit3DProfileAdvanced(chi2,par_init,par,npar,&_t[0],&_s[0],E,E0);
A = par[0];
B = par[1];
D = par[2];
- xl0 = par[3] * X0;
+ xl0 = par[3] * X0[0]; //probably shower start is in ecal
delete[] E;
@@ -827,7 +834,7 @@ int ClusterShapes::fit3DProfile(float& chi2, float& E0, float& A, float& B, floa
//=============================================================================
float ClusterShapes::getChi2Fit3DProfileSimple(float a, float b, float c, float d,
- float X0, float Rm) {
+ float* X0, float* Rm) {
float chi2 = 0.0;
@@ -845,7 +852,7 @@ float ClusterShapes::getChi2Fit3DProfileSimple(float a, float b, float c, float
//=============================================================================
float ClusterShapes::getChi2Fit3DProfileAdvanced(float E0, float a, float b, float d,
- float t0, float X0, float Rm) {
+ float t0, float* X0, float* Rm) {
float chi2 = 0.0;
@@ -897,18 +904,18 @@ int ClusterShapes::FitHelix(int max_iter, int status_out, int parametrisation,
int ClusterShapes::FitHelix(int max_iter, int status_out, int parametrisation,
double* parameter, double* dparameter, double& chi2,
- double& distmax, int direction) {
+ double& distmax, int direction) {
// FIXME: version with double typed parameters needed 2006/06/10 OW
if (_nHits < 3) {
- std::cout << "ClusterShapes : helix fit impossible, two few points" ;
- std::cout << std::endl;
- for (int i = 0; i < 5; ++i) {
- parameter[i] = 0.0;
- dparameter[i] = 0.0;
- }
- return 1;
+ std::cout << "ClusterShapes : helix fit impossible, two few points" ;
+ std::cout << std::endl;
+ for (int i = 0; i < 5; ++i) {
+ parameter[i] = 0.0;
+ dparameter[i] = 0.0;
+ }
+ return 1;
}
// find initial parameters
@@ -919,14 +926,14 @@ int ClusterShapes::FitHelix(int max_iter, int status_out, int parametrisation,
// 1st loop
for (int i(0); i < _nHits; ++i) {
- double Rz = sqrt(_xHit[i]*_xHit[i] + _yHit[i]*_yHit[i]);
- if (Rz < Rmin) {
- Rmin = Rz;
- i1 = i;
- }
- if (Rz > Rmax) {
- Rmax = Rz;
- }
+ double Rz = sqrt(_xHit[i]*_xHit[i] + _yHit[i]*_yHit[i]);
+ if (Rz < Rmin) {
+ Rmin = Rz;
+ i1 = i;
+ }
+ if (Rz > Rmax) {
+ Rmax = Rz;
+ }
}
@@ -934,8 +941,8 @@ int ClusterShapes::FitHelix(int max_iter, int status_out, int parametrisation,
// debug
/*
- for (int i(0); i < _nHits; ++i) std::cout << i << " " << _xHit[i] << " " << _yHit[i] << " " << _zHit[i] << std::endl;
- std::cout << std::endl << Rmin << " " << Rmax << " " << i1 << std::endl;
+ for (int i(0); i < _nHits; ++i) std::cout << i << " " << _xHit[i] << " " << _yHit[i] << " " << _zHit[i] << std::endl;
+ std::cout << std::endl << Rmin << " " << Rmax << " " << i1 << std::endl;
*/
@@ -954,8 +961,8 @@ int ClusterShapes::FitHelix(int max_iter, int status_out, int parametrisation,
if ((Rz > Lower) && (Rz < Upper)) {
double dZ = fabs(_zHit[i]-_zHit[i1]);
if (dZ < dZmin) {
- dZmin = dZ;
- i3 = i;
+ dZmin = dZ;
+ i3 = i;
}
}
}
@@ -977,35 +984,35 @@ int ClusterShapes::FitHelix(int max_iter, int status_out, int parametrisation,
// 3d loop
for (int i(0); i < _nHits; ++i) {
- if (_zHit[i] >= z1 && _zHit[i] <= z3) {
- double Rz = sqrt(_xHit[i]*_xHit[i] + _yHit[i]*_yHit[i]);
- double dRz = fabs(Rz - Rref);
- if (dRz < dRmin) {
- i2 = i;
- dRmin = dRz;
- }
+ if (_zHit[i] >= z1 && _zHit[i] <= z3) {
+ double Rz = sqrt(_xHit[i]*_xHit[i] + _yHit[i]*_yHit[i]);
+ double dRz = fabs(Rz - Rref);
+ if (dRz < dRmin) {
+ i2 = i;
+ dRmin = dRz;
}
+ }
}
- int problematic = 0;
+ //int problematic = 0;
if (i2 < 0 || i2 == i1 || i2 == i3) {
- problematic = 1;
+ //problematic = 1;
// std::cout << "here we are " << std::endl;
for (int i(0); i < _nHits; ++i) {
if (i != i1 && i != i3) {
- i2 = i;
- if (_zHit[i2] < z1) {
- int itemp = i1;
- i1 = i2;
- i2 = itemp;
- }
- else if (_zHit[i2] > z3) {
- int itemp = i3;
- i3 = i2;
- i2 = itemp;
- }
- break;
+ i2 = i;
+ if (_zHit[i2] < z1) {
+ int itemp = i1;
+ i1 = i2;
+ i2 = itemp;
+ }
+ else if (_zHit[i2] > z3) {
+ int itemp = i3;
+ i3 = i2;
+ i2 = itemp;
+ }
+ break;
}
}
// std::cout << i1 << " " << i2 << " " << i3 << std::endl;
@@ -1024,23 +1031,23 @@ int ClusterShapes::FitHelix(int max_iter, int status_out, int parametrisation,
double time;
if (det == 0.) {
- time = 500.;
+ time = 500.;
}
else {
- gsl_matrix* A = gsl_matrix_alloc(2,2);
- gsl_vector* B = gsl_vector_alloc(2);
- gsl_vector* T = gsl_vector_alloc(2);
- gsl_matrix_set(A,0,0,ax);
- gsl_matrix_set(A,0,1,-axp);
- gsl_matrix_set(A,1,0,ay);
- gsl_matrix_set(A,1,1,-ayp);
- gsl_vector_set(B,0,x0p-x0);
- gsl_vector_set(B,1,y0p-y0);
- gsl_linalg_HH_solve(A,B,T);
- time = gsl_vector_get(T,0);
- gsl_matrix_free(A);
- gsl_vector_free(B);
- gsl_vector_free(T);
+ gsl_matrix* A = gsl_matrix_alloc(2,2);
+ gsl_vector* B = gsl_vector_alloc(2);
+ gsl_vector* T = gsl_vector_alloc(2);
+ gsl_matrix_set(A,0,0,ax);
+ gsl_matrix_set(A,0,1,-axp);
+ gsl_matrix_set(A,1,0,ay);
+ gsl_matrix_set(A,1,1,-ayp);
+ gsl_vector_set(B,0,x0p-x0);
+ gsl_vector_set(B,1,y0p-y0);
+ gsl_linalg_HH_solve(A,B,T);
+ time = gsl_vector_get(T,0);
+ gsl_matrix_free(A);
+ gsl_vector_free(B);
+ gsl_vector_free(T);
}
double X0 = x0 + ax*time;
@@ -1052,44 +1059,44 @@ int ClusterShapes::FitHelix(int max_iter, int status_out, int parametrisation,
double R0 = sqrt(dX*dX + dY*dY);
/*
- if (problematic == 1) {
- std::cout << i1 << " " << i2 << " " << i3 << std::endl;
- std::cout << _xHit[i1] << " " << _yHit[i1] << " " << _zHit[i1] << std::endl;
- std::cout << _xHit[i2] << " " << _yHit[i2] << " " << _zHit[i2] << std::endl;
- std::cout << _xHit[i3] << " " << _yHit[i3] << " " << _zHit[i3] << std::endl;
- std::cout << "R0 = " << R0 << std::endl;
- }
+ if (problematic == 1) {
+ std::cout << i1 << " " << i2 << " " << i3 << std::endl;
+ std::cout << _xHit[i1] << " " << _yHit[i1] << " " << _zHit[i1] << std::endl;
+ std::cout << _xHit[i2] << " " << _yHit[i2] << " " << _zHit[i2] << std::endl;
+ std::cout << _xHit[i3] << " " << _yHit[i3] << " " << _zHit[i3] << std::endl;
+ std::cout << "R0 = " << R0 << std::endl;
+ }
*/
double phi1 = (double)atan2(_yHit[i1]-Y0,_xHit[i1]-X0);
double phi2 = (double)atan2(_yHit[i2]-Y0,_xHit[i2]-X0);
double phi3 = (double)atan2(_yHit[i3]-Y0,_xHit[i3]-X0);
-// testing bz > 0 hypothesis
+ // testing bz > 0 hypothesis
if ( phi1 > phi2 )
- phi2 = phi2 + 2.0*M_PI;
+ phi2 = phi2 + 2.0*M_PI;
if ( phi1 > phi3 )
- phi3 = phi3 + 2.0*M_PI;
+ phi3 = phi3 + 2.0*M_PI;
if ( phi2 > phi3 )
- phi3 = phi3 + 2.0*M_PI;
+ phi3 = phi3 + 2.0*M_PI;
double bz_plus = (phi3 - phi1) / (_zHit[i3]-_zHit[i1]);
double phi0_plus = phi1 - bz_plus * _zHit[i1];
double dphi_plus = fabs( bz_plus * _zHit[i2] + phi0_plus - phi2 );
-// testing bz < 0 hypothesis
+ // testing bz < 0 hypothesis
phi1 = (double)atan2(_yHit[i1]-Y0,_xHit[i1]-X0);
phi2 = (double)atan2(_yHit[i2]-Y0,_xHit[i2]-X0);
phi3 = (double)atan2(_yHit[i3]-Y0,_xHit[i3]-X0);
if ( phi1 < phi2 )
- phi2 = phi2 - 2.0*M_PI;
+ phi2 = phi2 - 2.0*M_PI;
if ( phi1 < phi3 )
- phi3 = phi3 - 2.0*M_PI;
+ phi3 = phi3 - 2.0*M_PI;
if ( phi2 < phi3 )
- phi3 = phi3 - 2.0*M_PI;
+ phi3 = phi3 - 2.0*M_PI;
double bz_minus = (phi3 - phi1) / (_zHit[i3]-_zHit[i1]);
double phi0_minus = phi1 - bz_minus * _zHit[i1];
@@ -1099,12 +1106,12 @@ int ClusterShapes::FitHelix(int max_iter, int status_out, int parametrisation,
double phi0;
if (dphi_plus < dphi_minus) {
- bz = bz_plus;
- phi0 = phi0_plus;
+ bz = bz_plus;
+ phi0 = phi0_plus;
}
else {
- bz = bz_minus;
- phi0 = phi0_minus;
+ bz = bz_minus;
+ phi0 = phi0_minus;
}
@@ -1133,11 +1140,11 @@ int ClusterShapes::FitHelix(int max_iter, int status_out, int parametrisation,
// debug
/*
- X0 = -1205.28;
- Y0 = 175.317;
- R0 = 1217.97;
- bz = 0.00326074;
- phi0 = -0.144444;
+ X0 = -1205.28;
+ Y0 = 175.317;
+ R0 = 1217.97;
+ bz = 0.00326074;
+ phi0 = -0.144444;
*/
@@ -1176,17 +1183,17 @@ int ClusterShapes::FitHelix(int max_iter, int status_out, int parametrisation,
// debug
/*
- std::cout << std::setprecision(6) << "InitFitCalculated (d0,z0,phi0,omega,tanL) = " << "(" << d0 << "," << z0 << "," << Phi0 << "," << omega << "," << tanL << ")"
- << " " << "sign(omega) = " << direction << std::endl;
+ std::cout << std::setprecision(6) << "InitFitCalculated (d0,z0,phi0,omega,tanL) = " << "(" << d0 << "," << z0 << "," << Phi0 << "," << omega << "," << tanL << ")"
+ << " " << "sign(omega) = " << direction << std::endl;
*/
// debug
/*
- d0 = 0.00016512;
- z0 = 0.000853511;
- Phi0 = 1.11974;
- omega = -4.22171e-05;
- tanL = -0.33436;
+ d0 = 0.00016512;
+ z0 = 0.000853511;
+ Phi0 = 1.11974;
+ omega = -4.22171e-05;
+ tanL = -0.33436;
*/
@@ -1223,7 +1230,7 @@ int ClusterShapes::FitHelix(int max_iter, int status_out, int parametrisation,
for (int ipoint(0); ipoint < _nHits; ipoint++) {
double distRPZ[2];
double Dist = DistanceHelix(_xHit[ipoint],_yHit[ipoint],_zHit[ipoint],
- X0,Y0,R0,bz,phi0,distRPZ);
+ X0,Y0,R0,bz,phi0,distRPZ);
double chi2rphi = distRPZ[0]/_exHit[ipoint];
chi2rphi = chi2rphi*chi2rphi;
double chi2z = distRPZ[1]/_ezHit[ipoint];
@@ -1306,7 +1313,11 @@ int ClusterShapes::FitHelix(int max_iter, int status_out, int parametrisation,
} while ( status==GSL_CONTINUE && iter < max_iter);
- gsl_multifit_covar (s->J, rel_error, covar);
+ //fg: jacobian has been dropped from gsl_multifit_fdfsolver in gsl 2:
+ gsl_matrix * J = gsl_matrix_alloc(s->fdf->n, s->fdf->p);
+ gsl_multifit_fdfsolver_jac( s, J);
+ gsl_multifit_covar( J, rel_error, covar );
+ // gsl_multifit_covar (s->J, rel_error, covar);
@@ -1347,7 +1358,7 @@ int ClusterShapes::FitHelix(int max_iter, int status_out, int parametrisation,
for (int ipoint(0); ipoint < _nHits; ipoint++) {
double distRPZ[2];
double Dist = DistanceHelix(_xHit[ipoint],_yHit[ipoint],_zHit[ipoint],
- X0,Y0,R0,bz,phi0,distRPZ);
+ X0,Y0,R0,bz,phi0,distRPZ);
double chi2rphi = distRPZ[0]/_exHit[ipoint];
chi2rphi = chi2rphi*chi2rphi;
double chi2z = distRPZ[1]/_ezHit[ipoint];
@@ -1446,23 +1457,23 @@ void ClusterShapes::findElipsoid() {
void ClusterShapes::findGravity() {
- _totAmpl = 0. ;
- for (int i(0); i < 3; ++i) {
- _analogGravity[i] = 0.0 ;
- }
- for (int i(0); i < _nHits; ++i) {
- _totAmpl+=_aHit[i] ;
- _analogGravity[0]+=_aHit[i]*_xHit[i] ;
- _analogGravity[1]+=_aHit[i]*_yHit[i] ;
- _analogGravity[2]+=_aHit[i]*_zHit[i] ;
- }
- for (int i(0); i < 3; ++i) {
- _analogGravity[i]/=_totAmpl ;
- }
- _xgr = _analogGravity[0];
- _ygr = _analogGravity[1];
- _zgr = _analogGravity[2];
- _ifNotGravity = 0;
+ _totAmpl = 0. ;
+ for (int i(0); i < 3; ++i) {
+ _analogGravity[i] = 0.0 ;
+ }
+ for (int i(0); i < _nHits; ++i) {
+ _totAmpl+=_aHit[i] ;
+ _analogGravity[0]+=_aHit[i]*_xHit[i] ;
+ _analogGravity[1]+=_aHit[i]*_yHit[i] ;
+ _analogGravity[2]+=_aHit[i]*_zHit[i] ;
+ }
+ for (int i(0); i < 3; ++i) {
+ _analogGravity[i]/=_totAmpl ;
+ }
+ _xgr = _analogGravity[0];
+ _ygr = _analogGravity[1];
+ _zgr = _analogGravity[2];
+ _ifNotGravity = 0;
}
//=============================================================================
@@ -1476,27 +1487,27 @@ void ClusterShapes::findInertia() {
findGravity();
for (int i(0); i < 3; ++i) {
- for (int j(0); j < 3; ++j) {
- aIne[i][j] = 0.0;
- }
+ for (int j(0); j < 3; ++j) {
+ aIne[i][j] = 0.0;
+ }
}
for (int i(0); i < _nHits; ++i) {
- float dX = _xHit[i] - _analogGravity[0];
- float dY = _yHit[i] - _analogGravity[1];
- float dZ = _zHit[i] - _analogGravity[2];
- aIne[0][0] += _aHit[i]*(dY*dY+dZ*dZ);
- aIne[1][1] += _aHit[i]*(dX*dX+dZ*dZ);
- aIne[2][2] += _aHit[i]*(dX*dX+dY*dY);
- aIne[0][1] -= _aHit[i]*dX*dY;
- aIne[0][2] -= _aHit[i]*dX*dZ;
- aIne[1][2] -= _aHit[i]*dY*dZ;
+ float dX = _xHit[i] - _analogGravity[0];
+ float dY = _yHit[i] - _analogGravity[1];
+ float dZ = _zHit[i] - _analogGravity[2];
+ aIne[0][0] += _aHit[i]*(dY*dY+dZ*dZ);
+ aIne[1][1] += _aHit[i]*(dX*dX+dZ*dZ);
+ aIne[2][2] += _aHit[i]*(dX*dX+dY*dY);
+ aIne[0][1] -= _aHit[i]*dX*dY;
+ aIne[0][2] -= _aHit[i]*dX*dZ;
+ aIne[1][2] -= _aHit[i]*dY*dZ;
}
for (int i(0); i < 2; ++i) {
- for (int j = i+1; j < 3; ++j) {
- aIne[j][i] = aIne[i][j];
- }
+ for (int j = i+1; j < 3; ++j) {
+ aIne[j][i] = aIne[i][j];
+ }
}
//****************************************
// analog Inertia
@@ -1523,13 +1534,13 @@ void ClusterShapes::findInertia() {
radius2 = 0.;
for (int i(0); i < 3; ++i) {
- _radius += _analogGravity[i]*_analogGravity[i];
- radius2 += (_analogGravity[i]+_VecAnalogInertia[i])*(_analogGravity[i]+_VecAnalogInertia[i]);
+ _radius += _analogGravity[i]*_analogGravity[i];
+ radius2 += (_analogGravity[i]+_VecAnalogInertia[i])*(_analogGravity[i]+_VecAnalogInertia[i]);
}
if ( radius2 < _radius) {
- for (int i(0); i < 3; ++i)
- _VecAnalogInertia[i] = - _VecAnalogInertia[i];
+ for (int i(0); i < 3; ++i)
+ _VecAnalogInertia[i] = - _VecAnalogInertia[i];
}
_radius = sqrt(_radius);
@@ -1563,99 +1574,99 @@ void ClusterShapes::findWidth() {
float ClusterShapes::findDistance(int i) {
- float cx = 0.0;
- float cy = 0.0;
- float cz = 0.0;
- float dx = 0.0;
- float dy = 0.0;
- float dz = 0.0;
- cx = _VecAnalogInertia[0] ;
- cy = _VecAnalogInertia[1] ;
- cz = _VecAnalogInertia[2] ;
- dx = _analogGravity[0] - _xHit[i] ;
- dy = _analogGravity[1] - _yHit[i] ;
- dz = _analogGravity[2] - _zHit[i] ;
- float tx = cy*dz - cz*dy ;
- float ty = cz*dx - cx*dz ;
- float tz = cx*dy - cy*dx ;
- float tt = sqrt(tx*tx+ty*ty+tz*tz) ;
- float ti = sqrt(cx*cx+cy*cy+cz*cz) ;
- float f = tt / ti ;
- return f ;
+ float cx = 0.0;
+ float cy = 0.0;
+ float cz = 0.0;
+ float dx = 0.0;
+ float dy = 0.0;
+ float dz = 0.0;
+ cx = _VecAnalogInertia[0] ;
+ cy = _VecAnalogInertia[1] ;
+ cz = _VecAnalogInertia[2] ;
+ dx = _analogGravity[0] - _xHit[i] ;
+ dy = _analogGravity[1] - _yHit[i] ;
+ dz = _analogGravity[2] - _zHit[i] ;
+ float tx = cy*dz - cz*dy ;
+ float ty = cz*dx - cx*dz ;
+ float tz = cx*dy - cy*dx ;
+ float tt = sqrt(tx*tx+ty*ty+tz*tz) ;
+ float ti = sqrt(cx*cx+cy*cy+cz*cz) ;
+ float f = tt / ti ;
+ return f ;
}
/**
- Function sdist(xp,yp,zp,cx,cy,cz,xv,yv,zv)
-c----------------------------------------------------------------------
-c Distance from line to point
-c xp, yp, zp -- point is at the line
-c xv, yv, zv -- point is out of line
-********************************************************************
-* Last update V.L.Morgunov 08-Apr-2002 *
-********************************************************************
- real xp,yp,zp,cx,cy,cz,xv,yv,zv,t1,t2,t3,tt,sdist
-
- t1 = cy*(zp-zv)-cz*(yp-yv)
- t2 = cz*(xp-xv)-cx*(zp-zv)
- t3 = cx*(yp-yv)-cy*(xp-xv)
- tt = sqrt(cx**2+cy**2+cz**2)
- sdist = sqrt(t1**2+t2**2+t3**2)/tt
-
- return
- end
+ Function sdist(xp,yp,zp,cx,cy,cz,xv,yv,zv)
+ c----------------------------------------------------------------------
+ c Distance from line to point
+ c xp, yp, zp -- point is at the line
+ c xv, yv, zv -- point is out of line
+ ********************************************************************
+ * Last update V.L.Morgunov 08-Apr-2002 *
+ ********************************************************************
+ real xp,yp,zp,cx,cy,cz,xv,yv,zv,t1,t2,t3,tt,sdist
+
+ t1 = cy*(zp-zv)-cz*(yp-yv)
+ t2 = cz*(xp-xv)-cx*(zp-zv)
+ t3 = cx*(yp-yv)-cy*(xp-xv)
+ tt = sqrt(cx**2+cy**2+cz**2)
+ sdist = sqrt(t1**2+t2**2+t3**2)/tt
+
+ return
+ end
*/
//=============================================================================
float ClusterShapes::vecProduct(float * x1, float * x2) {
- float x1abs(0.);
- float x2abs(0.);
- float prod(0.);
+ float x1abs(0.);
+ float x2abs(0.);
+ float prod(0.);
- for (int i(0); i < 3; ++i) {
- x1abs += x1[i]*x1[i];
- x2abs += x2[i]*x2[i];
- prod += x1[i]*x2[i];
- }
+ for (int i(0); i < 3; ++i) {
+ x1abs += x1[i]*x1[i];
+ x2abs += x2[i]*x2[i];
+ prod += x1[i]*x2[i];
+ }
- x1abs = sqrt(x1abs);
- x2abs = sqrt(x2abs);
+ x1abs = sqrt(x1abs);
+ x2abs = sqrt(x2abs);
- if (x1abs > 0.0 && x2abs > 0.0) {
- prod = prod/(x1abs*x2abs);
- }
- else {
- prod = 0.;
- }
+ if (x1abs > 0.0 && x2abs > 0.0) {
+ prod = prod/(x1abs*x2abs);
+ }
+ else {
+ prod = 0.;
+ }
- return prod;
+ return prod;
}
//=============================================================================
float ClusterShapes::vecProject(float * x, float * axis) {
- float axisabs(0.);
- float prod(0.);
- for (int i(0); i < 3; ++i) {
- axisabs += axis[i]*axis[i];
- prod += x[i]*axis[i];
- }
- axisabs = sqrt(axisabs);
- if (axisabs > 0.0 ) {
- prod = prod/axisabs;
- }
- else {
- prod = 0.;
- }
- return prod;
+ float axisabs(0.);
+ float prod(0.);
+ for (int i(0); i < 3; ++i) {
+ axisabs += axis[i]*axis[i];
+ prod += x[i]*axis[i];
+ }
+ axisabs = sqrt(axisabs);
+ if (axisabs > 0.0 ) {
+ prod = prod/axisabs;
+ }
+ else {
+ prod = 0.;
+ }
+ return prod;
}
//=============================================================================
double ClusterShapes::DistanceHelix(double x, double y, double z, double X0, double Y0,
- double R0, double bz, double phi0, double * distRPZ) {
+ double R0, double bz, double phi0, double * distRPZ) {
double phi = atan2(y-Y0,x-X0);
double R = sqrt( (y-Y0)*(y-Y0) + (x-X0)*(x-X0) );
@@ -1695,9 +1706,8 @@ double ClusterShapes::DistanceHelix(double x, double y, double z, double X0, dou
//=============================================================================
-int ClusterShapes::transformToEigensystem(float* xStart, int& index_xStart, float X0,
- float Rm) {
-
+int ClusterShapes::transformToEigensystem(float* xStart, int& index_xStart, float* X0, float* Rm) {
+
if (_ifNotInertia == 1) findInertia();
float MainAxis[3];
@@ -1711,9 +1721,16 @@ int ClusterShapes::transformToEigensystem(float* xStart, int& index_xStart, floa
MainCentre[1] = _analogGravity[1];
MainCentre[2] = _analogGravity[2];
+ //analoginertia looks something wrong!
+ //change it to the direction of CoG
+ //float ll=sqrt(MainCentre[0]*MainCentre[0]+MainCentre[1]*MainCentre[1]+MainCentre[2]*MainCentre[2]);
+ //MainAxis[0]=MainCentre[0]/ll;
+ //MainAxis[1]=MainCentre[1]/ll;
+ //MainAxis[2]=MainCentre[2]/ll;
+
int ifirst = 0;
float xx[3];
- float prodmin = 0.0;
+ float prodmin = 1.0e+100;
int index = 0;
for (int i(0); i < _nHits; ++i) {
@@ -1731,27 +1748,73 @@ int ClusterShapes::transformToEigensystem(float* xStart, int& index_xStart, floa
xStart[1] = MainCentre[1] + prodmin*MainAxis[1];
xStart[2] = MainCentre[2] + prodmin*MainAxis[2];
index_xStart = index;
+
+ float l=sqrt(MainAxis[0]*MainAxis[0]+MainAxis[1]*MainAxis[1]+MainAxis[2]*MainAxis[2]);
+ //float ll=sqrt(xStart[0]*xStart[0]+xStart[1]*xStart[1]+xStart[2]*xStart[2]);
+ //std::cout << "xstart: " << index_xStart << " " << prodmin << " " << xStart[0] << std::endl;
+ //calculate the surface of hcal
+ float ecalrad=2.058000000e+03; //in mm
+ float plugz=2.650000000e+03; //in mm
+
+ float tmpcos=MainAxis[2]/l;
+ float tmpsin=sqrt(MainAxis[0]*MainAxis[0]+MainAxis[1]*MainAxis[1])/l;
+ float detend=0.0;
+ if(fabs(xStart[2])<2.450000000e+03){ //if in the barrel
+ detend=(ecalrad-sqrt(xStart[0]*xStart[0]+xStart[1]*xStart[1]))/tmpsin;
+ }else{ //if in plug
+ detend=(plugz-fabs(xStart[2]))/fabs(tmpcos);
+ }
+ if(detend<0.0) detend=0.0;
+
+ // std::cout << "check length: " << detend << " "
+ // << xStart[0]/ll << " " << xStart[1]/ll << " " << xStart[2]/ll << " "
+ // << MainAxis[0]/l << " " << MainAxis[1]/l << " " << MainAxis[1]/l << std::endl;
+
for (int i(0); i < _nHits; ++i) {
xx[0] = _xHit[i] - xStart[0];
xx[1] = _yHit[i] - xStart[1];
xx[2] = _zHit[i] - xStart[2];
float xx2(0.);
for (int j(0); j < 3; ++j) xx2 += xx[j]*xx[j];
-
+
_xl[i] = 0.001 + vecProject(xx,MainAxis);
- _xt[i] = sqrt(std::max(0.0,xx2 + 0.1 - _xl[i]*_xl[i]));
+ _xt[i] = sqrt(std::max(0.0,xx2 + 0.01 - _xl[i]*_xl[i]));
// std::cout << i << " " << _xl[i] << " " << _xt[i] << " " << _aHit[i] << " "
// << std::endl;
}
-
- for (int i = 0; i < _nHits; ++i) {
- _t[i] = _xl[i]/X0;
- _s[i] = _xt[i]/Rm;
+
+ //first, check ecal and solve wrong behaviour
+ //std::cout << "param: " << X0[0] << " " << X0[1] << " " << Rm[0] << " " << Rm[1] << std::endl;
+ for (int i = 0; i < _nHits; ++i) {
+ //if(_types[i]==0 || _types[i]==3){ //if the hit is in Ecal
+ _t[i] = _xl[i]/X0[0];
+ _s[i] = _xt[i]/Rm[0];
+ if(detend<_xl[i]){
+ //std::cout << "something wrong!! " << detend << " " << _xl[i] << " " << _t[i] << std::endl;
+ //detend = _xl[i]; //to avoid wrong behaviour
+ }
+ //}
}
-
+
+ //second, check hcal
+ /*for (int i = 0; i < _nHits; ++i) {
+ if(_types[i]==1 || _types[i]==4){ //if the hit is in Hcal
+ if(_xl[i]>detend){
+ _t[i] = detend/X0[0]+(_xl[i]-detend)/X0[1];
+ _s[i] = _xt[i]/Rm[1];
+ }else{
+ _t[i] = _xl[i]/X0[0];
+ _s[i] = _xt[i]/Rm[0];
+ }
+ }
+ // std::cout << _types[i] << " " << _xl[i] << " "
+ // << _xl[i]+sqrt(xStart[0]*xStart[0]+xStart[1]*xStart[1]+xStart[2]*xStart[2]) << " "
+ // << _t[i] << " " << _s[i] << std::endl;
+ }*/
+
_ifNotEigensystem = 0;
-
+
return 0; // no error messages at the moment
}
@@ -1764,7 +1827,7 @@ float ClusterShapes::calculateChi2Fit3DProfileSimple(float a, float b, float c,
// ClusterShapes::transformToEigensystem needs to be executed before
float chi2 = 0.0;
- float Ampl = 0.0;
+ float Ampl = 0.0;
for (int i(0); i < _nHits; ++i) {
// old definition of Ampl and chi2
@@ -1804,11 +1867,11 @@ float ClusterShapes::calculateChi2Fit3DProfileAdvanced(float E0, float a, float
// debug
/*
- std::cout << "OUT : " << Ampl << " " << E0 << " " << a << " " << b << " "
- << d << " " << t0 << " " << _t[i] << " " << invG(a) << " "
- << pow(b*(_t[i]-t0),a-1) << " " << exp(-b*(_t[i]-t0)) << " "
- << exp(-d*_s[i])
- << std::endl;
+ std::cout << "OUT : " << Ampl << " " << E0 << " " << a << " " << b << " "
+ << d << " " << t0 << " " << _t[i] << " " << invG(a) << " "
+ << pow(b*(_t[i]-t0),a-1) << " " << exp(-b*(_t[i]-t0)) << " "
+ << exp(-d*_s[i])
+ << std::endl;
*/
chi2 += ((Ampl - _aHit[i])*(Ampl - _aHit[i]))/(_aHit[i]*_aHit[i]);
@@ -1932,6 +1995,8 @@ int ClusterShapes::fit3DProfileAdvanced(float& chi2, double* par_init, double* p
const gsl_multifit_fdfsolver_type* T = gsl_multifit_fdfsolver_lmsder;
+ //std::cout << T << " " << _nHits << " " << npar << std::endl;
+
gsl_multifit_fdfsolver* Solver = gsl_multifit_fdfsolver_alloc(T,_nHits,npar);
gsl_matrix* covar = gsl_matrix_alloc(npar,npar); // covariance matrix
@@ -1951,6 +2016,8 @@ int ClusterShapes::fit3DProfileAdvanced(float& chi2, double* par_init, double* p
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();
@@ -1958,9 +2025,9 @@ int ClusterShapes::fit3DProfileAdvanced(float& chi2, double* par_init, double* p
// perform fit
do {
iter++;
- std::cout << "Multidimensional Fit Iteration started ... ... ";
+ //std::cout << "Multidimensional Fit Iteration started ... ... ";
status = gsl_multifit_fdfsolver_iterate(Solver);
- std::cout << "--- DONE ---" << std::endl;
+ //std::cout << "--- DONE ---" << std::endl;
if (status) break;
status = gsl_multifit_test_delta (Solver->dx,Solver->x,abs_error,rel_error);
@@ -1974,15 +2041,19 @@ int ClusterShapes::fit3DProfileAdvanced(float& chi2, double* par_init, double* p
par[3] = (float)gsl_vector_get(Solver->x,3);
std::cout << "Status of multidimensional fit : " << status << " "
- << "Iterations : " << iter << std::endl;
+ << "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;
+ << par[1] << "\t" << "D : " << par[2] << "\t" << "t0 : " << par[3]
+ << std::endl << std::endl;
*/
} while ( status==GSL_CONTINUE && iter < max_iter);
- gsl_multifit_covar (Solver->J,rel_error,covar);
+ //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
@@ -2003,3 +2074,222 @@ int ClusterShapes::fit3DProfileAdvanced(float& chi2, double* par_init, double* p
//=============================================================================
+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;
+}
+
+