9 # define SQR(x) (x)*(x)
31 fHit =
new int[ nhits ];
32 fHitf =
new double[ nhits ];
60 for (
int i=0; i <
fNhits; i++) {
62 JObject::oid_t
id =
getHitID( hits, i ) ;
64 h.
id = (JObject::oid_t)
id;
74 static uint32_t Nwarns=0;
76 std::cout <<
"Warning: DFCALCluster : corrupted cluster hit "
79 std::cout <<
"Last warning!!! (further warnings supressed)"
112 double t2EWeight = 0, tEWeight = 0;
116 double hitEnergy = hitList->
hit[ih].
E*frac;
120 t2EWeight += hitEnergy * hitList->
hit[ih].
t * hitList->
hit[ih].
t;
121 tEWeight += hitEnergy * hitList->
hit[ih].
t;
127 double eMax=0, timeMax=0;
137 centroid.SetXYZ(0., 0., fcalFaceZ );
140 #ifdef LOG2_WEIGHTING
145 double weightSum = 0.0;
146 double centerWeight = 0.0;
147 double neighborMaxWeight = 0.0;
148 double logFraction = exp(-0.23*(energy));
149 double currentOffset = log(energy)/7.0 + 3.7;
159 weight = currentOffset + log(hitList->
hit[ih].
E*frac/energy);
161 centerWeight = weight;
165 (neighborMaxWeight < weight)? weight:neighborMaxWeight;
168 xc += hitList->
hit[ih].
x*weight;
169 yc += hitList->
hit[ih].
y*weight;
177 if (neighborMaxWeight > 0) {
178 xc += (logFraction-1)*(centerWeight-neighborMaxWeight)
180 yc += (logFraction-1)*(centerWeight-neighborMaxWeight)
182 weightSum += (logFraction-1)*(centerWeight-neighborMaxWeight);
184 centroid.SetX(xc/weightSum);
185 centroid.SetY(yc/weightSum);
190 xc += hitList->
hit[ih].
x*(hitList->
hit[ih].
E*frac);
191 yc += hitList->
hit[ih].
y*(hitList->
hit[ih].
E*frac);
194 centroid.SetX(xc/energy);
195 centroid.SetY(yc/energy);
209 double x = hitList->
hit[ih].
x;
210 double y = hitList->
hit[ih].
y;
212 MOM1x += hitList->
hit[ih].
E*frac*
x;
213 MOM1y += hitList->
hit[ih].
E*frac*
y;
214 MOM2x += hitList->
hit[ih].
E*frac*
SQR(x);
215 MOM2y += hitList->
hit[ih].
E*frac*
SQR(y);
217 double phi = atan2( centroid.y() , centroid.x() );
218 double u = x*cos(phi) + y*
sin(phi);
219 double v =-x*
sin(phi) + y*cos(phi);
220 MOM1u += hitList->
hit[ih].
E*frac*
u;
221 MOM1v += hitList->
hit[ih].
E*frac*v;
222 MOM2u += hitList->
hit[ih].
E*frac*
SQR(u);
223 MOM2v += hitList->
hit[ih].
E*frac*
SQR(v);
226 bool something_changed =
false;
227 if (fabs(energy-
fEnergy) > 0.001) {
229 something_changed =
true;
231 if (fabs(eMax-
fEmax) > 0.001) {
235 something_changed =
true;
237 if (fabs(centroid.x()-
fCentroid.x()) > 0.1 ||
238 fabs(centroid.y()-
fCentroid.y()) > 0.1) {
240 something_changed =
true;
243 if (something_changed) {
247 fRMS_t =
sqrt( t2EWeight - ( tEWeight * tEWeight ) );
255 for (
int ih = 0; ih < hitList->
nhits; ih++) {
261 return something_changed;
266 double& Eallowed,
double& Eexpected,
267 double fcalMidplaneZ)
const
271 Eallowed = Eexpected = 0;
274 double x = hitList->
hit[ihit].
x;
275 double y = hitList->
hit[ihit].
y;
283 double u = x*cos(phi) + y*
sin(phi);
284 double v =-x*
sin(phi) + y*cos(phi);
286 double uVar = vVar+
SQR(
SQR(8*theta));
287 double vTail = 4.5+0.9*log(
fEnergy+0.05);
288 double uTail = vTail+
SQR(10*theta);
289 double core = exp(-0.5*
SQR(
SQR(u-u0)/uVar +
SQR(v-v0)/vVar));
290 double tail = exp(-
sqrt(
SQR((u-u0)/uTail)+
SQR((v-v0)/vTail)));
292 Eallowed = 2*
fEmax*core + (0.2+0.5*log(
fEmax+1.))*tail;
294 if ((dist <= 4.) && (Eallowed <
fEmax) ) {
295 std::cerr <<
"Warning: FCAL cluster Eallowed value out of range!\n";
double getHitT(const userhits_t *const hitList, const int ihit) const
static double blockLength()
DFCALCluster(const int nhits)
void shower_profile(const userhits_t *const hitList, const int ihit, double &Eallowed, double &Eexpected, double fcalMidplaneZ) const
void saveHits(const userhits_t *const hit)
bool update(const userhits_t *const hitList, double fcalFaceZ)
int addHit(const int ihit, const double frac)
double getHitY(const userhits_t *const hitList, const int ihit) const
double getHitX(const userhits_t *const hitList, const int ihit) const
vector< DFCALClusterHit_t > my_hits
double getHitE(const userhits_t *const hitList, const int ihit) const
#define MAX_SHOWER_RADIUS
oid_t getHitID(const userhits_t *const hitList, const int ihit) const
int getHitCh(const userhits_t *const hitList, const int ihit) const
double getHitIntOverPeak(const userhits_t *const hitList, const int ihit) const