Bug Summary

File:/volatile/halld/gluex/nightly/2024-03-22/Linux_CentOS7.7-x86_64-gcc4.8.5/halld_recon/src/libraries/BCAL/DBCALCluster_factory.cc
Location:line 797, column 164
Description:Value stored to 'connected' is never read

Annotated Source Code

1/*
2 * DBCALCluster_factory.cc
3 *
4 * Created by Matthew Shepherd on 3/12/11.
5 *
6 */
7
8#include <iostream>
9
10using namespace std;
11
12#include "DANA/DApplication.h"
13#include "BCAL/DBCALGeometry.h"
14#include "BCAL/DBCALHit.h"
15#include "BCAL/DBCALUnifiedHit.h"
16
17#include "BCAL/DBCALCluster_factory.h"
18
19#include "units.h"
20#include <TMath.h>
21
22bool PointSort( const DBCALPoint* p1, const DBCALPoint* p2 ){
23
24 return ( p1->E() > p2->E() );
25}
26
27bool ClusterSort( const DBCALCluster* c1, const DBCALCluster* c2 ){
28
29 return ( c1->E() > c2->E() );
30}
31
32DBCALCluster_factory::DBCALCluster_factory() :
33 m_mergeSig( 5 ),
34 m_moliereRadius( 3.7*k_cm ),
35 m_clust_hit_timecut ( 20.0*k_nsec ),
36 m_timeCut( 8.0*k_nsec ){
37
38 // The phi and theta direction inclusion curves are described in:
39 // http://argus.phys.uregina.ca/gluex/DocDB/0029/002998/003/CAL_meeting_may5.pdf.
40 // The theta direction inclusion curve needs to be a function of theta. C1_parm and
41 // C2_parm are parameters [0] and [1] in dtheta_inclusion_curve.
42 }
43
44jerror_t
45DBCALCluster_factory::init(void){
46
47 m_BCALGeom = NULL__null;
48 return NOERROR;
49
50}
51
52jerror_t
53DBCALCluster_factory::fini( void ){
54
55 return NOERROR;
56}
57
58jerror_t DBCALCluster_factory::brun(JEventLoop *loop, int32_t runnumber) {
59 DApplication* app = dynamic_cast<DApplication*>(loop->GetJApplication());
60 DGeometry* geom = app->GetDGeometry(runnumber);
61 geom->GetTargetZ(m_z_target_center);
62
63 // load BCAL Geometry
64 vector<const DBCALGeometry *> BCALGeomVec;
65 loop->Get(BCALGeomVec);
66 if(BCALGeomVec.size() == 0)
67 throw JException("Could not load DBCALGeometry object!");
68 m_BCALGeom = BCALGeomVec[0];
69
70
71 loop->GetCalib("/BCAL/effective_velocities", effective_velocities);
72
73 loop->GetCalib("/BCAL/attenuation_parameters",attenuation_parameters);
74
75 BCALCLUSTERVERBOSE = 0;
76 gPARMS->SetDefaultParameter("BCALCLUSTERVERBOSE", BCALCLUSTERVERBOSE, "VERBOSE level for BCAL Cluster overlap success and conditions");
77 //command line parameter to investigate what points are being added to clusters and what clusters are being merged together. // Track fitterer helper class
78
79
80 vector<const DTrackFitter *> fitters;
81 loop->Get(fitters);
82
83 if(fitters.size()<1){
84 _DBG_std::cerr<<"libraries/BCAL/DBCALCluster_factory.cc"<<
":"<<84<<" "<<"Unable to get a DTrackFinder object!"<<endl;
85 return RESOURCE_UNAVAILABLE;
86 }
87
88 fitter = fitters[0];
89
90 return NOERROR;
91}
92
93jerror_t
94DBCALCluster_factory::evnt( JEventLoop *loop, uint64_t eventnumber ){
95
96 vector< const DBCALPoint* > twoEndPoint;
97 vector< const DBCALPoint* > usedPoints;
98 loop->Get(twoEndPoint);
99
100 // Want to add singled-ended hits to the Clusters.
101
102 // Looking for hits that are single-ended.
103
104 vector< const DBCALUnifiedHit* > hits;
105 loop->Get(hits);
106
107 vector< const DTrackWireBased* > tracks;
108 loop->Get(tracks);
109
110 // first arrange the list of hits so they are grouped by cell
111 map< int, vector< const DBCALUnifiedHit* > > cellHitMap;
112 for( vector< const DBCALUnifiedHit* >::const_iterator hitPtr = hits.begin();
113 hitPtr != hits.end();
114 ++hitPtr ){
115
116 const DBCALUnifiedHit& hit = (**hitPtr);
117
118 int id = m_BCALGeom->cellId( hit.module, hit.layer, hit.sector );
119
120 if( cellHitMap.find( id ) == cellHitMap.end() ){
121
122 cellHitMap[id] = vector< const DBCALUnifiedHit* >();
123 }
124
125 cellHitMap[id].push_back( *hitPtr );
126 }
127
128 // now we should try to add on single-ended hits ...
129 vector< const DBCALUnifiedHit* > single_ended_hits;
130
131 for( map< int, vector< const DBCALUnifiedHit* > >::iterator mapItr = cellHitMap.begin();
132 mapItr != cellHitMap.end();
133 ++mapItr ){
134
135 if( mapItr->second.size() == 1 ){
136 // only one hit in the cell
137
138 const DBCALUnifiedHit* hit = mapItr->second[0];
139
140 single_ended_hits.push_back(hit);
141
142 }
143 }
144
145 vector<DBCALCluster*> clusters = clusterize( twoEndPoint, usedPoints, single_ended_hits, tracks );
146
147 // load our vector of clusters into the factory member data
148 for( vector<DBCALCluster*>::iterator clust = clusters.begin();
149 clust != clusters.end();
150 ++clust ){
151
152 if( isnan((**clust).t()) == 1 || isnan((**clust).phi()) == 1 || isnan((**clust).theta()) == 1 ) continue;
153 // put in an energy threshold for clusters
154 if( (**clust).E() < 5*k_MeV ) {
155 delete *clust;
156 continue;
157 }
158 vector<const DBCALPoint*>points=(**clust).points();
159 for (unsigned int i=0;i<points.size();i++){
160 (**clust).AddAssociatedObject(points[i]);
161 }
162 _data.push_back(*clust);
163 }
164 return NOERROR;
165}
166
167vector<DBCALCluster*>
168DBCALCluster_factory::clusterize( vector< const DBCALPoint* > points , vector< const DBCALPoint* > usedPoints , vector< const DBCALUnifiedHit* > hits, vector< const DTrackWireBased* > tracks ) const {
169
170 // first sort the points by energy
171 sort( points.begin(), points.end(), PointSort );
172
173 vector<DBCALCluster*> clusters(0);
174
175 // ahh.. more hard coded numbers that should probably
176 // come from a database or something
177 float seedThresh = 1.*k_GeV;
178 float minSeed = 10*k_MeV;
179 //We have a big problem with noise in the outer layer of the detector
180 //(the noise is the greatest in the outer layer, since the number of SiPMs
181 //being summed is also the greatest here).
182 //Thus there are a lot of DBCALPoint's in this layer that are pure noise hits.
183 //The simplest way to deal with this is to prevent outer layer points
184 //from seeding clusters. So hits in the outer layer can be associated
185 //with existing clusters, but cannot create their own cluster.
186 //This is okay since since isolated hits in the outer layer
187 //is not really a signature we expect for many physical showers.
188 //However, if a hit is sufficiently energetic, it is unlikely to be a noise
189 //hit. For this reason, we allow 4th layer hits to seed clusters,
190 //but we need a different (higher) minimum seed energy.
191 float layer4_minSeed = 50*k_MeV;
192 float tracked_phi = 0.;
193 float matched_dphi = .175;
194 float matched_dtheta = .087;
195
196 while( seedThresh > minSeed ) {
197
198 bool usedPoint = false;
199
200 for( vector< const DBCALPoint* >::iterator pt = points.begin();
201 pt != points.end();
202 ++pt ){
203
204 // first see if point should be added to an existing
205 // cluster
206
207 int q = 0;
208
209 // Check if a point is matched to a track
210 for( vector< const DTrackWireBased* >::iterator trk = tracks.begin();
211 trk != tracks.end();
212 ++trk ){
213 DVector3 track_pos(0.0, 0.0, 0.0);
214 double point_r = (**pt).r();
215 double point_z = (**pt).z();
216 vector<DTrackFitter::Extrapolation_t>extrapolations=(*trk)->extrapolations.at(SYS_BCAL);
217 if (fitter->ExtrapolateToRadius(point_r,extrapolations,track_pos)){
218 double dPhi=track_pos.Phi()-(**pt).phi();
219 if (dPhi<-M_PI3.14159265358979323846) dPhi+=2.*M_PI3.14159265358979323846;
220 if (dPhi>M_PI3.14159265358979323846) dPhi-=2.*M_PI3.14159265358979323846;
221 double point_theta_global = fabs(atan2(point_r,point_z + m_z_target_center )); // convert point z-position origin to global frame to match tracks origin
222 double dTheta = fabs(point_theta_global - track_pos.Theta());
223 matched_dphi=0.175+0.175*exp(-0.8*extrapolations[0].momentum.Mag());
224 if(fabs(dPhi) < matched_dphi && dTheta < matched_dtheta){
225 q = 1; // if point and track are matched then set q = 1
226 tracked_phi = extrapolations[0].position.Phi();
227 break;
228 }
229 }
230 }
231
232 for( vector<DBCALCluster*>::iterator clust = clusters.begin();
233 clust != clusters.end();
234 ++clust ){
235
236 if((**clust).Q()==1){
237 if(overlap_charged( **clust,*pt, tracked_phi ) ){
238 usedPoints.push_back( *pt );
239 int point_q = 1;
240 (**clust).addPoint( *pt, point_q );
241 points.erase( pt );
242 usedPoint = true;
243 break;
244 }
245 }
246 if( overlap( **clust, *pt ) ){
247 if (q==1 && (**pt).layer()!=1) q=0;
248 // assign point q=1 if it's in layer 1 because track matching tends to be improved in layer 1 than later layers where the cluster seed is. This would allow us to jump into the charged clustering routines on the fly.
249 usedPoints.push_back( *pt );
250 (**clust).addPoint( *pt , q);
251 points.erase( pt );
252 usedPoint = true;
253 break;
254 }
255 // once we erase a point the iterator is no longer useful
256 // and we start the loop over, so that a point doesn't get added to
257 // multiple clusters. We will recycle through points later to
258 // check if a point was added to its closest cluster.
259 }
260
261 if( usedPoint ) break;
262
263 // if the point doesn't overlap with a cluster see if it can become a
264 // new seed
265 if( (**pt).E() > seedThresh && ((**pt).layer() != 4 || (**pt).E() > layer4_minSeed) ){
266 clusters.push_back(new DBCALCluster( *pt, m_z_target_center, q, m_BCALGeom ) );
267 usedPoints.push_back( *pt );
268 points.erase( pt );
269 usedPoint = true;
270 break;
271 }
272 }
273
274 recycle_points( usedPoints, clusters);
275 // recycle through points that were added to a cluster and check if they
276 // were added to their closest cluster. If they weren't then we remove
277 // the point from its original cluster and add it to its closest cluster.
278
279 double point_reatten_E = 0.;
280 merge( clusters, point_reatten_E );
281 // lower the threshold to look for new seeds if none of
282 // the existing points were used as new clusters or assigned
283 // to existing clusters
284 if( !usedPoint ) seedThresh /= 2;
285 }
286
287 // add the single-ended hits that overlap with a cluster that was made from points
288 for( vector< const DBCALUnifiedHit* >::iterator ht = hits.begin();
289 ht != hits.end();
290 ++ht){
291 bool usedHit = false;
292
293 for( vector<DBCALCluster*>::iterator clust = clusters.begin();
294 clust != clusters.end();
295 ++clust ){
296
297 if( overlap( **clust, *ht ) ){
298
299 int channel_calib = 16*((**ht).module-1)+4*((**ht).layer-1)+(**ht).sector-1; // need to use cellID for objects in DBCALGeometry but the CCDB uses a different channel numbering scheme, so use channel_calib when accessing CCDB tables.
300
301 // given the location of the cluster, we need the best guess
302 // for z with respect to target at this radius
303
304 double z = (**clust).rho()*cos((**clust).theta()) + m_z_target_center;
305 double d = ( ((**ht).end == 0) ? (z - m_BCALGeom->GetBCAL_center() + m_BCALGeom->GetBCAL_length()/2.0) : (m_BCALGeom->GetBCAL_center() + m_BCALGeom->GetBCAL_length()/2.0 - z)); // d gives the distance to upstream or downstream end of BCAL depending on where the hit was with respect to the cluster z position.
306 double lambda = attenuation_parameters[channel_calib][0];
307 double hit_E = (**ht).E;
308 double hit_E_unattenuated = hit_E/exp(-d/lambda); // hit energy unattenuated wrt the cluster z position
309
310 (**clust).addHit( *ht, hit_E_unattenuated );
311 usedHit = true;
312 }
313 if( usedHit ) break;
314 }
315 }
316 return clusters;
317}
318
319void
320DBCALCluster_factory::recycle_points( vector<const DBCALPoint*> usedPoints, vector<DBCALCluster*>& clusters) const{
321
322 if ( clusters.size() <= 1 ) return;
323
324 int q = 2;
325
326 sort( clusters.begin(), clusters.end(), ClusterSort );
327
328 for( vector<const DBCALPoint*>::const_iterator usedpt = usedPoints.begin();
329 usedpt != usedPoints.end();
330 ++usedpt ){
331
332 bool got_overlap=false;
333 double min_phi=1e6;
334
335 for( vector<DBCALCluster*>::iterator clust = clusters.begin();
336 clust != clusters.end();
337 ++clust ){
338
339 if( overlap( **clust, *usedpt ) ){
340 got_overlap=true;
341
342 float deltaPhi = (**clust).phi() - (*usedpt)->phi();
343 if (deltaPhi<-M_PI3.14159265358979323846) deltaPhi+=2.*M_PI3.14159265358979323846;
344 if (deltaPhi>M_PI3.14159265358979323846) deltaPhi-=2.*M_PI3.14159265358979323846;
345 if (fabs(deltaPhi)<min_phi){
346 min_phi=fabs(deltaPhi);
347 }
348 }
349 }
350
351 if(got_overlap==false) break;
352
353 // Find the points closest cluster in distance along the sphere and in phi
354 for( vector<DBCALCluster*>::iterator clust = clusters.begin();
355 clust != clusters.end();
356 ++clust ){
357 bool best_clust = false;
358 vector<const DBCALPoint*>associated_points=(**clust).points();
359
360 float deltaPhi = (**clust).phi() - (*usedpt)->phi();
361 if (deltaPhi<-M_PI3.14159265358979323846) deltaPhi+=2.*M_PI3.14159265358979323846;
362 if (deltaPhi>M_PI3.14159265358979323846) deltaPhi-=2.*M_PI3.14159265358979323846;
363 deltaPhi=fabs(deltaPhi);
364
365 for(unsigned int j = 0 ; j < associated_points.size(); j++){
366 // Check to see if the point we are comparing to the cluster
367 // is already in that cluster.
368 if (fabs((*usedpt)->E()-associated_points[j]->E())<1e-4
369 && fabs(deltaPhi-min_phi)<1e-4) best_clust=true;
370 if(BCALCLUSTERVERBOSE>1)cout << " clust E = " << (**clust).E() <<" assoc point E = " << associated_points[j]->E() << " points E = " << (*usedpt)->E() << " clust match = " << best_clust << endl;
371 }
372 if(best_clust==true) break;
373 // if the point was originally placed in its "best" cluster then we don't want to touch it.
374 if(best_clust==0){
375 int added_point = 0;
376 int removed_point = 0;
377 for(unsigned int i = 0 ; i < associated_points.size(); i++){
378 bool point_match = (fabs((*usedpt)->E()-associated_points[i]->E())<1e-4);
379 if( point_match==0 && added_point==0 && fabs(deltaPhi-min_phi)<1e-4){
380 (**clust).addPoint( *usedpt , q );
381 // if the point found a closer cluster then we add it to the closer cluster.
382 // The point is now an associated object of the closer cluster.
383 added_point=1;
384 }
385 if( point_match==1 && removed_point==0 && fabs(deltaPhi-min_phi)>1e-4){
386 (**clust).removePoint( *usedpt );
387 // Now we remove the point from its original cluster since it has been added
388 // to its closest cluster. The point is no longer an associated object of
389 // the original cluster.
390 removed_point=1;
391 }
392 }
393 }
394 }
395 }
396}
397
398void
399DBCALCluster_factory::merge( vector<DBCALCluster*>& clusters, double point_reatten_E ) const {
400
401 if( clusters.size() <= 1 ) return;
402
403 sort( clusters.begin(), clusters.end(), ClusterSort );
404
405 bool stillMerging = true;
406
407 float low_z_lim = -100.;
408 float high_z_lim = 500.;
409
410 while( stillMerging ){
411
412 stillMerging = false;
413 for( vector<DBCALCluster*>::iterator hClust = clusters.begin();
414 hClust != clusters.end() - 1;
415 ++hClust ){
416
417 vector<const DBCALPoint*>hClust_points=(**hClust).points();
418
419 for( vector<DBCALCluster*>::iterator lClust = hClust + 1;
420 lClust != clusters.end();
421 ++lClust ){
422
423 vector<const DBCALPoint*>lClust_points=(**lClust).points();
424 vector<const DBCALPoint*>hClust_points=(**hClust).points();
425
426 if( overlap( **hClust, **lClust ) ){
427
428 point_reatten_E = 0.;
429
430 if (hClust_points.size() == 1) {
431
432 for( unsigned int i = 0 ; i < hClust_points.size() ; i++){
433
434 if (hClust_points[i]->z() > low_z_lim && hClust_points[i]->z() < high_z_lim) point_reatten_E = 0.;
435 else {
436 int channel_calib = 16*(hClust_points[i]->module()-1)+4*(hClust_points[i]->layer()-1)+hClust_points[i]->sector()-1;
437
438 double fibLen = m_BCALGeom->GetBCAL_length();
439
440 double point_z = hClust_points[i]->z();
441 double zLocal = point_z + m_z_target_center - m_BCALGeom->GetBCAL_center();
442
443 double dUp = 0.5 * fibLen + zLocal;
444 double dDown = 0.5 * fibLen - zLocal;
445 if (dUp>fibLen) dUp=fibLen;
446 if (dUp<0) dUp=0;
447 if (dDown>fibLen) dDown=fibLen;
448 if (dDown<0) dDown=0;
449
450 double lambda = attenuation_parameters[channel_calib][0];
451 double attUp = exp( -dUp / lambda );
452 double attDown = exp( -dDown / lambda );
453
454 double US_unatten_E = hClust_points[i]->E_US()*attUp;
455 double DS_unatten_E = hClust_points[i]->E_DS()*attDown;
456
457 double zLocal_clust = m_BCALGeom->GetBCAL_inner_rad()/tan((**lClust).theta()) + m_z_target_center - m_BCALGeom->GetBCAL_center();
458 double dUp_clust = 0.5 * fibLen + zLocal_clust;
459 double dDown_clust = 0.5 * fibLen - zLocal_clust;
460
461 double attUp_clust = exp( -dUp_clust / lambda );
462 double attDown_clust = exp( -dDown_clust / lambda );
463
464 double US_reattn_E = US_unatten_E/attUp_clust;
465 double DS_reattn_E = DS_unatten_E/attDown_clust;
466 point_reatten_E = 0.5 * ( US_reattn_E + DS_reattn_E);
467
468 }
469 }
470 }
471
472 if (lClust_points.size() == 1) {
473
474 for( unsigned int i = 0 ; i < lClust_points.size() ; i++){
475
476 if (lClust_points[i]->z() > low_z_lim && lClust_points[i]->z() < high_z_lim) point_reatten_E = 0.;
477 else{
478 int channel_calib = 16*(lClust_points[i]->module()-1)+4*(lClust_points[i]->layer()-1)+lClust_points[i]->sector()-1;
479
480 double fibLen = m_BCALGeom->GetBCAL_length();
481
482 double point_z = lClust_points[i]->z();
483 double zLocal = point_z + m_z_target_center - m_BCALGeom->GetBCAL_center();
484
485 double dUp = 0.5 * fibLen + zLocal;
486 double dDown = 0.5 * fibLen - zLocal;
487 if (dUp>fibLen) dUp=fibLen;
488 if (dUp<0) dUp=0;
489 if (dDown>fibLen) dDown=fibLen;
490 if (dDown<0) dDown=0;
491
492 double lambda = attenuation_parameters[channel_calib][0];
493 double attUp = exp( -dUp / lambda );
494 double attDown = exp( -dDown / lambda );
495
496 double US_unatten_E = lClust_points[i]->E_US()*attUp;
497 double DS_unatten_E = lClust_points[i]->E_DS()*attDown;
498
499 double zLocal_clust = m_BCALGeom->GetBCAL_inner_rad()/tan((**hClust).theta()) + m_z_target_center - m_BCALGeom->GetBCAL_center();
500 double dUp_clust = 0.5 * fibLen + zLocal_clust;
501 double dDown_clust = 0.5 * fibLen - zLocal_clust;
502
503 double attUp_clust = exp( -dUp_clust / lambda );
504 double attDown_clust = exp( -dDown_clust / lambda );
505
506 double US_reattn_E = US_unatten_E/attUp_clust;
507 double DS_reattn_E = DS_unatten_E/attDown_clust;
508 point_reatten_E = 0.5 * ( US_reattn_E + DS_reattn_E);
509
510 }
511 }
512 }
513
514 if( (**lClust).Q() == 1 && (**hClust).Q() == 0) {
515 (**lClust).mergeClust(**hClust, point_reatten_E);
516 delete *hClust;
517 clusters.erase( hClust );
518 }
519
520 else {
521 (**hClust).mergeClust(**lClust, point_reatten_E);
522 delete *lClust;
523 clusters.erase( lClust );
524 }
525
526 // now iterators are invalid and we need to bail out of loops
527 stillMerging = true;
528 break;
529 }
530 }
531 if( stillMerging ) break;
532 }
533 }
534}
535
536bool
537DBCALCluster_factory::overlap( const DBCALCluster& highEClust,
538 const DBCALCluster& lowEClust ) const {
539
540 float sigTheta = fabs( highEClust.theta() - lowEClust.theta() ) /
541 sqrt( highEClust.sigTheta() * highEClust.sigTheta() +
542 lowEClust.sigTheta() * lowEClust.sigTheta() );
543
544 // difference in phi is tricky due to overlap at 0/2pi
545 // order based on phi and then take the minimum of the difference
546 // and the difference with 2pi added to the smallest
547
548 float deltaPhi = highEClust.phi() - lowEClust.phi();
549 float deltaPhiAlt = ( highEClust.phi() > lowEClust.phi() ?
550 highEClust.phi() - lowEClust.phi() - 2*TMath::Pi() :
551 lowEClust.phi() - highEClust.phi() - 2*TMath::Pi() );
552
553 deltaPhi = min( fabs( deltaPhi ), fabs( deltaPhiAlt ) );
554
555 float sigPhi = deltaPhi /
556 sqrt( highEClust.sigPhi() * highEClust.sigPhi() +
557 lowEClust.sigPhi() * lowEClust.sigPhi() );
558
559 //We can't rely entirely on sigTheta and sigPhi as defined above.
560 //For high-energy clusters, the position uncertainties will be very small,
561 //so sigTheta/sigPhi will be large, and clusters may not merge properly.
562 //To fix this, force clusters to merge if delta_z and delta_phi are both
563 //very small. This is hopefully only a temporary fix.
564
565 //deltaPhi_force_merge and delta_z_force_merge were determined by looking
566 //at the separation of decay photons from pi0's from a pythia sample.
567 //There are no events where the decay photons have separation
568 //(delta_phi < 0.2 && delta_z < 25 cm), so in most cases it should be safe
569 //to merge clusters together if they are so close.
570 const double deltaPhi_force_merge = 0.1; //radians
571 const double delta_z_force_merge = 15.0*k_cm;
572
573 //A major cause of extra clusters are lower energy hits, which have poor
574 //z-resolution and so are not properly merged. Treat low energy
575 //clusters (< 40 MeV) as a special case. Again, hopefully this is only
576 //a temporary fix until we have a more comprehensive solution.
577 const double delta_z_force_merge_low_E = 40.0*k_cm;
578 const double low_E = .04*k_GeV;
579
580 double z1 = m_BCALGeom->GetBCAL_inner_rad()/tan(highEClust.theta());
581 double z2 = m_BCALGeom->GetBCAL_inner_rad()/tan(lowEClust.theta());
582 double delta_z = fabs(z1-z2);
583
584 bool theta_match = (sigTheta < m_mergeSig) || (delta_z < delta_z_force_merge) || (delta_z < delta_z_force_merge_low_E && lowEClust.E() < low_E);
585
586 bool phi_match = (sigPhi < m_mergeSig) || (deltaPhi < deltaPhi_force_merge);
587
588 //very loose cut to check that the two clusters are in time
589 bool time_match = (highEClust.t() - lowEClust.t()) < m_timeCut;
590
591 if(BCALCLUSTERVERBOSE>1) cout << " clust merge: " << " theta match success = " << theta_match << " phi match = " << phi_match << " time match = " << time_match << " high E = " << highEClust.E() << " low E = " << lowEClust.E() << " highE z = " << z1 << " lowE z = " << z2 << " deltaTheta = " << fabs(highEClust.theta()-lowEClust.theta()) << " sigTheta = " << sigTheta << " highE sigTheta = " << highEClust.sigTheta() << " lowE sigTheta = " << lowEClust.sigTheta() << endl;
592
593 vector<const DBCALPoint*> highE_points;
594 highE_points = (highEClust).points();
595
596 vector<const DBCALPoint*> lowE_points;
597 lowE_points = (lowEClust).points();
598
599 double highE_summed_z = 0.;
600 double highE_summed_phi = 0.;
601 double highE_summed_zphi = 0.;
602 double highE_summed_z_sq = 0.;
603 double highE_slope = 0.;
604 double highE_y_intercept = 0.;
605
606 double lowE_summed_z = 0.;
607 double lowE_summed_phi = 0.;
608 double lowE_summed_zphi = 0.;
609 double lowE_summed_z_sq = 0.;
610 double lowE_slope = 0.;
611 double lowE_y_intercept = 0.;
612
613 int connected = 0;
614// double z_match = 50.;
615 double slope_match = 0.01;
616 double intercept_match = 1.8;
617 double deltaPhi_match = 0.2;
618
619 int lowE_global_sector = 0;
620 int highE_global_sector = 0;
621 int lowE_point_layer = 0;
622
623 for(unsigned int i = 0 ; i < lowE_points.size() ; i ++){
624 // adjust the points phi position to be close to the cluster phi position at the 0/2pi phi boundary
625 if(lowEClust.phi() > lowE_points[i]->phi() ){
626 if( fabs( lowEClust.phi() - lowE_points[i]->phi() - 2*TMath::Pi() ) < TMath::Pi() ) lowE_points[i]->add2Pi();
627 }
628 else{
629 if( fabs( lowE_points[i]->phi() - lowEClust.phi() - 2*TMath::Pi() ) < TMath::Pi() ) lowE_points[i]->sub2Pi();
630 }
631
632 // compute quantities to be used to calculate the direction of the lower energy cluster if we need it for merging.
633 lowE_summed_z += lowE_points[i]->z();
634 lowE_summed_phi += lowE_points[i]->phi();;
635 lowE_summed_zphi += lowE_points[i]->z()*lowE_points[i]->phi();
636 lowE_summed_z_sq += lowE_points[i]->z()*lowE_points[i]->z();
637 if(lowE_points.size()==1) {
638 lowE_global_sector = 4*(lowE_points[i]->module()-1) + lowE_points[i]->sector();
639 lowE_point_layer = lowE_points[i]->layer();
640 }
641 }
642
643 for(unsigned int i = 0 ; i < highE_points.size() ; i ++){
644 // adjust the points phi position to be close to the cluster phi position at the 0/2pi phi boundary
645 if(highEClust.phi() > highE_points[i]->phi() ){
646 if( fabs( highEClust.phi() - highE_points[i]->phi() - 2*TMath::Pi() ) < TMath::Pi() ) highE_points[i]->add2Pi();
647 }
648 else{
649 if( fabs( highE_points[i]->phi() - highEClust.phi() - 2*TMath::Pi() ) < TMath::Pi() ) highE_points[i]->sub2Pi();
650 }
651 // compute quantities to be used to calculate the direction of the higher energy cluster if we need it for merging.
652 highE_summed_z += highE_points[i]->z();
653 highE_summed_phi += highE_points[i]->phi();;
654 highE_summed_zphi += highE_points[i]->z()*highE_points[i]->phi();
655 highE_summed_z_sq += highE_points[i]->z()*highE_points[i]->z();
656 highE_global_sector = 4*(highE_points[i]->module()-1) + highE_points[i]->sector();
657 if(lowE_points.size()==1 && lowE_point_layer == highE_points[i]->layer() && ( lowE_global_sector+1 == highE_global_sector || lowE_global_sector-1 == highE_global_sector ) ) connected = 1; // clustesr that contain only a single point won't have any fit parameters and will make it hard for them to merge, this connected int will force a merge if a single point cluster is connected to a cluster without any points adjacent to it.
658 }
659
660 // calculate slopes and intercepts of the 2 clusters direction and if one of the clusters is matched to a track then we will require their fit parameter quantities
661 // to match for their merging criteria. This allows us to relax their phi position proximity merging criteria since split clusters matched to tracks tend to be
662 // further distributed in the azimuthal direction than neutral clusters.
663
664 highE_slope = (highE_summed_z*highE_summed_phi - highE_points.size()*highE_summed_zphi)/(highE_summed_z*highE_summed_z - highE_points.size()*highE_summed_z_sq);
665 highE_y_intercept = (highE_summed_zphi*highE_summed_z - highE_summed_phi*highE_summed_z_sq)/(highE_summed_z*highE_summed_z - highE_points.size()*highE_summed_z_sq);
666
667 lowE_slope = (lowE_summed_z*lowE_summed_phi - lowE_points.size()*lowE_summed_zphi)/(lowE_summed_z*lowE_summed_z - lowE_points.size()*lowE_summed_z_sq);
668 lowE_y_intercept = (lowE_summed_zphi*lowE_summed_z - lowE_summed_phi*lowE_summed_z_sq)/(lowE_summed_z*lowE_summed_z - lowE_points.size()*lowE_summed_z_sq);
669
670 double delta_slope = fabs(highE_slope - lowE_slope) ;
671 double delta_intercept = fabs(highE_y_intercept - lowE_y_intercept) ;
672
673 highE_points.clear();
674 lowE_points.clear();
675
676
677 // If both clusters trying to merge together were NOT matched to a track then use neutral clusterizer merging critera of theta and phi matching.
678 // If EITHER of the 2 clusters trying to merge together were amtched to a track then use the information about the direction of the cluster for merging.
679
680 if (highEClust.Q() == 0 && lowEClust.Q() == 0 ) return theta_match && phi_match && time_match;
681
682 else return ( ( delta_slope < slope_match && delta_intercept < intercept_match && deltaPhi < deltaPhi_match ) || connected == 1 ) ;
683
684
685}
686
687bool
688DBCALCluster_factory::overlap( const DBCALCluster& clust,
689 const DBCALPoint* point ) const {
690
691
692 float deltaTheta = fabs( clust.theta() - point->theta() );
693 /* sigTheta not used
694 float sigTheta = deltaTheta / sqrt( clust.sigTheta() * clust.sigTheta() +
695 point->sigTheta() * point->sigTheta() );
696 */
697
698 // difference in phi is tricky due to overlap at 0/2pi
699 // order based on phi and then take the minimum of the difference
700 // and the difference with 2pi added to the smallest
701
702 float deltaPhi = clust.phi() - point->phi();
703 float deltaPhiAlt = ( clust.phi() > point->phi() ?
704 clust.phi() - point->phi() - 2*TMath::Pi() :
705 point->phi() - clust.phi() - 2*TMath::Pi() );
706
707 deltaPhi = min( fabs( deltaPhi ), fabs( deltaPhiAlt ) );
708
709 /* sigPhi not used
710 float sigPhi = deltaPhi /
711 sqrt( clust.sigPhi() * clust.sigPhi() +
712 point->sigPhi() * point->sigPhi() );
713 */
714
715 float rho = ( clust.rho() + point->rho() ) / 2;
716 float theta = ( clust.theta() + point->theta() ) / 2;
717
718 float sep = sqrt( ( rho * deltaTheta ) * ( rho * deltaTheta ) +
719 ( rho * sin( theta ) * deltaPhi ) * ( rho * sin( theta ) * deltaPhi ) );
720
721 float sep_term1 = rho*deltaTheta;
722 float sep_term2 = rho*sin(theta)*deltaPhi;
723
724 //very loose cuts to make sure the two hits are in time
725 bool time_match = fabs(clust.t() - point->t()) < m_timeCut;
726
727 double clust_z = clust.rho()*cos(clust.theta());
728
729 //double c1 = C1_parm->Eval(clust_z);
730 double c1=23.389+19.093*tanh(-0.0104*(clust_z-201.722));
731
732 //double c2 = C2_parm->Eval(clust_z);
733 double c2=0.151+0.149*tanh(-0.016*(clust_z-275.194));
734
735 //dtheta_inclusion_curve->SetParameter(0,c1);
736 //dtheta_inclusion_curve->SetParameter(1,c2);
737
738 //double inclusion_val = sep_inclusion_curve->Eval(sep);
739 double inclusion_val=exp(-sep/30.)-0.1;
740
741 //double inclusion_val1 = dtheta_inclusion_curve->Eval(sep_term1);
742 double inclusion_val1=exp(-(sep_term1-0.1)/c1)-c2+.15;
743
744 //double inclusion_val2 = dphi_inclusion_curve->Eval(sep_term2);
745 double inclusion_val2=exp(-(sep_term2-2.)/2.5)-sep_term2*0.002+0.07;
746
747 // We consider fractional energy (point.E/Clust.E) as a function of spatial separation between
748 // a point and cluster to determine if the point should be included in the cluster.
749 // These distributions are tighter in the phihat direction than along thetahat. For more details
750 // on how the selection criteria for cluster,point overlap function go to logbook entry 3396018.
751
752 if(BCALCLUSTERVERBOSE>0) cout << "(m,l,s) = (" <<point->module()<<","<<point->layer()<<","<<point->sector()<<")" << " sep = " << sep << "sep1 = " << sep_term1 << " sep2 = " << sep_term2 << " inclusion value = " << inclusion_val << " inclusion val1= " << inclusion_val1 << " inclusion val2= " << inclusion_val2<< " time match = " << time_match << " clust E = " << clust.E() << " point E = " << point->E() << " energy ratio = " << point->E()/(point->E()+clust.E()) << " clust theta = " << clust.theta()*180./3.14159 << " point theta = " << point->theta()*180./3.14159 << " sep rho*deltaTheta = " << ( rho * deltaTheta ) << endl;
753
754 if(sep>m_moliereRadius && sep<7.*m_moliereRadius &&sep_term2>=2.*m_moliereRadius){
755 return ((point->E()/(point->E()+clust.E())) < inclusion_val1 ) && ((point->E()/(point->E()+clust.E())) < inclusion_val2 ) && time_match && deltaPhi*180./3.14159<10.;
756 }
757
758 else{
759 return ((point->E()/(point->E()+clust.E())) < (inclusion_val1+.2)) && sep < 10.*m_moliereRadius && time_match && sep_term2<2.*m_moliereRadius;
760 }
761
762}
763
764
765bool
766DBCALCluster_factory::overlap_charged( const DBCALCluster& clust,
767 const DBCALPoint* point, float tracked_phi) const {
768
769
770 // difference in phi is tricky due to overlap at 0/2pi
771 // order based on phi and then take the minimum of the difference
772 // and the difference with 2pi added to the smallest
773
774 float phiCut = 0.65417;
775
776 vector<const DBCALPoint*> assoc_points;
777 assoc_points = (clust).points();
778
779 double summed_r = 0.;
780 double summed_phi = 0.;
781 double summed_rphi = 0.;
782 double summed_r_sq = 0.;
783
784 double summed_z = 0.;
785 double summed_zphi = 0.;
786 double summed_z_sq = 0.;
787
788 double slope = 0.;
789 double y_intercept = 0.;
790
791 int point_global_sector = 4*(point->module()-1) + point->sector();
792 int point_layer = point->layer();
793 int connected = 0;
794
795 for(unsigned int i = 0 ; i < assoc_points.size() ; i ++){
796 int assoc_point_global_sector = 4*(assoc_points[i]->module() - 1) + assoc_points[i]->sector();
797 if( point_layer == assoc_points[i]->layer() && ( point_global_sector + 1 == assoc_point_global_sector || point_global_sector - 1 == assoc_point_global_sector) ) connected = 1;
Value stored to 'connected' is never read
798 summed_r += assoc_points[i]->r();
799 summed_z += assoc_points[i]->z();
800 if( tracked_phi > assoc_points[i]->phi() ){
801 if( fabs( tracked_phi - assoc_points[i]->phi() - 2*TMath::Pi() ) < TMath::Pi() ) assoc_points[i]->add2Pi();
802 }
803 else{
804 if( fabs( assoc_points[i]->phi() - tracked_phi - 2*TMath::Pi() ) < TMath::Pi() ) assoc_points[i]->sub2Pi();
805 }
806
807 summed_phi += assoc_points[i]->phi();
808 summed_rphi += assoc_points[i]->r()*assoc_points[i]->phi();
809 summed_r_sq += assoc_points[i]->r()*assoc_points[i]->r();
810 summed_zphi += assoc_points[i]->z()*assoc_points[i]->phi();
811 summed_z_sq += assoc_points[i]->z()*assoc_points[i]->z();
812
813 }
814
815 if(assoc_points.size()<2){
816 slope = (tracked_phi - summed_phi)/(m_BCALGeom->GetBCAL_inner_rad() - summed_r);
817 y_intercept = tracked_phi - slope*m_BCALGeom->GetBCAL_inner_rad();
818 }
819
820 else{
821 slope = (summed_z*summed_phi - assoc_points.size()*summed_zphi)/(summed_z*summed_z - assoc_points.size()*summed_z_sq);
822 y_intercept = (summed_zphi*summed_z - summed_phi*summed_z_sq)/(summed_z*summed_z - assoc_points.size()*summed_z_sq);
823 }
824
825 float fit_phi = 0.;
826
827 if(assoc_points.size() < 2) fit_phi = slope*point->r() + y_intercept;
828 else fit_phi = slope*point->z() + y_intercept;
829
830 assoc_points.clear();
831
832 float deltaPhi = fit_phi-point->phi();
833 float deltaPhiAlt = ( fit_phi > point->phi() ?
834 fit_phi - point->phi() - 2*TMath::Pi() :
835 point->phi() - fit_phi - 2*TMath::Pi() );
836
837 deltaPhi = min( fabs( deltaPhi ), fabs( deltaPhiAlt ) );
838
839 float rho = point->rho();
840 float theta = point->theta();
841
842 float deltaTheta = fabs( clust.theta() - point->theta() );
843
844 float sep = sqrt( ( rho * deltaTheta ) * ( rho * deltaTheta ) +
845 ( rho * sin( theta ) * deltaPhi ) * ( rho * sin( theta ) * deltaPhi ) );
846
847 float sep_term1 = rho*deltaTheta;
848 float sep_term2 = rho*sin(theta)*deltaPhi;
849
850 //very loose cuts to make sure the two hits are in time
851 bool time_match = fabs(clust.t() - point->t()) < m_timeCut;
852
853 bool phi_match = fabs( clust.phi() - point->phi() ) < phiCut;
854
855 double clust_z = clust.rho()*cos(clust.theta());
856
857 //double c1 = C1_parm->Eval(clust_z);
858 double c1=23.389+19.093*tanh(-0.0104*(clust_z-201.722));
859
860 //double c2 = C2_parm->Eval(clust_z);
861 double c2=0.151+0.149*tanh(-0.016*(clust_z-275.194));
862
863 //dtheta_inclusion_curve->SetParameter(0,c1);
864 //dtheta_inclusion_curve->SetParameter(1,c2);
865
866 //double inclusion_val = sep_inclusion_curve->Eval(sep);
867 double inclusion_val=exp(-sep/30.)-0.1;
868
869 //double inclusion_val1 = dtheta_inclusion_curve->Eval(sep_term1);
870 double inclusion_val1=exp(-(sep_term1-0.1)/c1)-c2+.15;
871
872 double inclusion_val2 = exp(-(sep_term2-2.)/1.5) - sep_term2*.007 + .15;
873
874 // We consider fractional energy (point.E/Clust.E) as a function of spatial separation between
875 // a point and cluster to determine if the point should be included in the cluster.
876 // These distributions are tighter in the phihat direction than along thetahat. For more details
877 // on how the selection criteria for cluster,point overlap function go to logbook entry 3396018.
878
879 if(BCALCLUSTERVERBOSE>1) cout << "(m,l,s) = (" <<point->module()<<","<<point->layer()<<","<<point->sector()<<")" << " sep = " << sep << "sep1 = " << sep_term1 << " sep2 = " << sep_term2 << " inclusion value = " << inclusion_val << " inclusion val1= " << inclusion_val1 << " inclusion val2= " << inclusion_val2<< " time match = " << time_match << " clust E = " << clust.E() << " point E = " << point->E() << " energy ratio = " << point->E()/(point->E()+clust.E()) << " clust theta = " << clust.theta()*180./3.14159 << " point theta = " << point->theta()*180./3.14159 << " sep rho*deltaTheta = " << ( rho * deltaTheta ) << endl;
880
881 if(sep>m_moliereRadius && sep<7.*m_moliereRadius &&sep_term2>=2.*m_moliereRadius){
882 return ((point->E()/(point->E()+clust.E())) < (inclusion_val1) ) && ((point->E()/(point->E()+clust.E())) < (inclusion_val2) ) && time_match && phi_match;
883 }
884
885 else{
886 return ((point->E()/(point->E()+clust.E())) < (inclusion_val1 + .2)) && sep < 10.*m_moliereRadius && time_match && sep_term2<2.*m_moliereRadius;
887 }
888
889 return connected == 1;
890
891}
892
893
894bool
895DBCALCluster_factory::overlap( const DBCALCluster& clust,
896 const DBCALUnifiedHit* hit ) const {
897
898 int cellId = m_BCALGeom->cellId( hit->module, hit->layer, hit->sector );
899
900 float cellPhi = m_BCALGeom->phi( cellId );
901 float cellSigPhi = m_BCALGeom->phiSize( cellId );
902
903 // annoying +- 2pi business to try to find the best delta phi
904
905 float deltaPhi = clust.phi() - cellPhi;
906 float deltaPhiAlt = ( clust.phi() > cellPhi ?
907 clust.phi() - cellPhi - 2*TMath::Pi() :
908 cellPhi - clust.phi() - 2*TMath::Pi() );
909 deltaPhi = min( fabs( deltaPhi ), fabs( deltaPhiAlt ) );
910
911 float sigPhi = deltaPhi /
912 sqrt( clust.sigPhi() * clust.sigPhi() + cellSigPhi * cellSigPhi );
913
914 int channel_calib = 16*(hit->module-1)+4*(hit->layer-1)+hit->sector-1; // need to use cellID for objects in DBCALGeometry but the CCDB uses a different channel numbering scheme, so use channel_calib when accessing CCDB tables.
915 // given the location of the cluster, we need the best guess
916 // for z with respect to target at this radius
917 double z = clust.rho()*cos(clust.theta()) + m_z_target_center;
918 double d = ( (hit->end == 0) ? (z - m_BCALGeom->GetBCAL_center() + m_BCALGeom->GetBCAL_length()/2.0) : (m_BCALGeom->GetBCAL_center() + m_BCALGeom->GetBCAL_length()/2.0 - z)); // d gives the distance to upstream or downstream end of BCAL depending on where the hit was with respect to the cluster z position.
919 double time_corr = hit->t - d/effective_velocities[channel_calib]; // hit time corrected to the interaction point in the bar.
920 double time_diff = TMath::Abs(clust.t() - time_corr); // time cut between cluster time and hit time - 20 ns is a very loose time cut.
921
922 return( sigPhi < m_mergeSig && time_diff < m_clust_hit_timecut );
923
924}