| File: | /home/sdobbs/work/clang/halld_recon/src/libraries/BCAL/DBCALCluster_factory.cc |
| Warning: | line 797, column 164 Value stored to 'connected' is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 1 | /* |
| 2 | * DBCALCluster_factory.cc |
| 3 | * |
| 4 | * Created by Matthew Shepherd on 3/12/11. |
| 5 | * |
| 6 | */ |
| 7 | |
| 8 | #include <iostream> |
| 9 | |
| 10 | using 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 | |
| 22 | bool PointSort( const DBCALPoint* p1, const DBCALPoint* p2 ){ |
| 23 | |
| 24 | return ( p1->E() > p2->E() ); |
| 25 | } |
| 26 | |
| 27 | bool ClusterSort( const DBCALCluster* c1, const DBCALCluster* c2 ){ |
| 28 | |
| 29 | return ( c1->E() > c2->E() ); |
| 30 | } |
| 31 | |
| 32 | DBCALCluster_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 | |
| 44 | jerror_t |
| 45 | DBCALCluster_factory::init(void){ |
| 46 | |
| 47 | m_BCALGeom = NULL__null; |
| 48 | return NOERROR; |
| 49 | |
| 50 | } |
| 51 | |
| 52 | jerror_t |
| 53 | DBCALCluster_factory::fini( void ){ |
| 54 | |
| 55 | return NOERROR; |
| 56 | } |
| 57 | |
| 58 | jerror_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 | |
| 93 | jerror_t |
| 94 | DBCALCluster_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 | |
| 167 | vector<DBCALCluster*> |
| 168 | DBCALCluster_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 | |
| 319 | void |
| 320 | DBCALCluster_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 | |
| 398 | void |
| 399 | DBCALCluster_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 | |
| 536 | bool |
| 537 | DBCALCluster_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 | |
| 687 | bool |
| 688 | DBCALCluster_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 | |
| 765 | bool |
| 766 | DBCALCluster_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 | |
| 894 | bool |
| 895 | DBCALCluster_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 | } |