analysis/EICAnalysis/TrackProjectionTools.C

0001 #include "TrackProjectionTools.h"
0002
0003 /* STL includes */
0004 #include <cassert>
0005
0006 /* Fun4All includes */
0007 #include <trackbase_historic/SvtxTrackMap_v1.h>
0008
0009 #include <phool/getClass.h>
0010
0011 #include <calobase/RawTowerGeomContainer.h>
0012 #include <calobase/RawTowerContainer.h>
0013 #include <calobase/RawTowerGeom.h>
0014 #include <calobase/RawTowerv1.h>
0015
0016 #include <calobase/RawClusterContainer.h>
0017 #include <calobase/RawCluster.h>
0018
0019 #include <g4vertex/GlobalVertexMap.h>
0020 #include <g4vertex/GlobalVertex.h>
0021
0022 #include <g4eval/CaloEvalStack.h>
0023 #include <g4eval/CaloRawTowerEval.h>
0024
0025 using namespace std;
0026
0027 TrackProjectionTools::TrackProjectionTools( PHCompositeNode *topNode )
0028 {
0029   _topNode = topNode;
0030 }
0031
0032
0033 SvtxTrack* TrackProjectionTools::FindClosestTrack( RawCluster* cluster, float& best_track_dr )
0034 {
0035
0036   /* best matching track */
0037   SvtxTrack* best_track = NULL;
0038   best_track_dr = NAN;
0039
0040   /* find name of calorimeter for this cluster */
0041   string caloname = "NONE";
0042
0043   /* C++11 range loop */
0044   for (auto& towit : cluster->get_towermap() )
0045     {
0046       caloname = RawTowerDefs::convert_caloid_to_name( RawTowerDefs::decode_caloid(towit.first) );
0047       break;
0048     }
0049
0050   /* Get track collection with all tracks in this event */
0051   SvtxTrackMap* trackmap =
0052     findNode::getClass<SvtxTrackMap>(_topNode,"SvtxTrackMap_FastSim");
0053   if (!trackmap)
0054     {
0055       cout << PHWHERE << "SvtxTrackMap node not found on node tree"
0056            << endl;
0057     }
0058
0059   /* Loop over all tracks from BARREL tracking and see if one points to the same
0060    * cluster as the reference clusters (i.e. matching ID in the same calorimeter) */
0061   /*
0062   for (SvtxTrackMap::ConstIter track_itr = trackmap->begin();
0063        track_itr != trackmap->end(); track_itr++)
0064     {
0065       SvtxTrack* track =  dynamic_cast<SvtxTrack*>(track_itr->second);
0066
0067       if ( caloname == "CEMC" &&
0068            track->get_cal_cluster_id(SvtxTrack::CEMC) == cluster->get_id() )
0069         {
0070           best_track = track;
0071         }
0072     }
0073   */
0074
0075   /* If track found with barrel tracking, return it here- if not, proceed with forward tracking below. */
0076   if ( best_track )
0077     return best_track;
0078
0079   /* Cluster / track matching for barrel calorimeters and tracking */
0080   float max_dr = 10;
0081
0082   /* cluster position for easy reference */
0083   float cx = cluster->get_x();
0084   float cy = cluster->get_y();
0085
0086   /* If track map found: Loop over all tracks to find best match for cluster (forward calorimeters)*/
0087   if ( trackmap &&
0088        ( caloname == "FEMC" || caloname == "EEMC" ) )
0089     {
0090       for (SvtxTrackMap::ConstIter track_itr = trackmap->begin();
0091            track_itr != trackmap->end(); track_itr++)
0092         {
0093           /* get pointer to track */
0094           SvtxTrack* track =  dynamic_cast<SvtxTrack*>(track_itr->second);
0095
0096           /* distance between track and cluster */
0097           float dr = NAN;
0098
0099           /* loop over track states (projections) sotred for this track */
0100           for (SvtxTrack::ConstStateIter state_itr = track->begin_states();
0101                state_itr != track->end_states(); state_itr++)
0102             {
0103               /* get pointer to current track state */
0104               SvtxTrackState *temp = dynamic_cast<SvtxTrackState*>(state_itr->second);
0105
0106               /* check if track state projection name matches calorimeter where cluster was found */
0107               if( (temp->get_name()==caloname) )
0108                 {
0109                   dr = sqrt( pow( cx - temp->get_x(), 2 ) + pow( cy - temp->get_y(), 2 ) );
0110           break;
0111                 }
0112             }
0113
0114           /* check dr and update best_track and best_track_dr if this track is closest to cluster */
0115           if ( ( best_track_dr != best_track_dr ) ||
0116            ( dr < max_dr &&
0117          dr < best_track_dr )
0118            )
0119             {
0120               best_track = track;
0121               best_track_dr = dr;
0122             }
0123         }
0124     }
0125
0126
0127   /* If track found with barrel tracking, return it here- if not, proceed with alternative barrel cluster-track matching below. */
0128   if ( best_track )
0129     return best_track;
0130
0131
0132   /* Cluster / track matching for barrel calorimeters and tracking */
0133   float max_dr_barrel = 10;
0134
0135   float ctheta = atan2( cluster->get_r() , cluster->get_z() );
0136   float ceta =  -log( tan( ctheta / 2.0 ) );
0137   float cphi = cluster->get_phi();
0138
0139   /* If track map found: Loop over all tracks to find best match for cluster (barrel calorimeters)*/
0140   if ( trackmap &&
0141        ( caloname == "CEMC" ) )
0142     {
0143       for (SvtxTrackMap::ConstIter track_itr = trackmap->begin();
0144            track_itr != trackmap->end(); track_itr++)
0145         {
0146           /* get pointer to track */
0147           SvtxTrack* track =  dynamic_cast<SvtxTrack*>(track_itr->second);
0148
0149           /* distance between track and cluster */
0150           float dr = NAN;
0151
0152           /* loop over track states (projections) sotred for this track */
0153           for (SvtxTrack::ConstStateIter state_itr = track->begin_states();
0154                state_itr != track->end_states(); state_itr++)
0155             {
0156               /* get pointer to current track state */
0157               SvtxTrackState *temp = dynamic_cast<SvtxTrackState*>(state_itr->second);
0158
0159               /* check if track state projection name matches calorimeter where cluster was found */
0160               if( (temp->get_name()==caloname) )
0161                 {
0162                   dr = sqrt( pow( ceta - temp->get_eta(), 2 ) + pow( cphi - temp->get_phi(), 2 ) );
0163                   break;
0164                 }
0165             }
0166
0167           /* check dr and update best_track and best_track_dr if this track is closest to cluster */
0168       if ( ( best_track_dr != best_track_dr ) ||
0169            (dr < max_dr_barrel &&
0170                dr < best_track_dr) )
0171             {
0172               best_track = track;
0173               best_track_dr = dr;
0174             }
0175         }
0176     }
0177
0178   return best_track;
0179 }
0180
0181 RawCluster*
0182 TrackProjectionTools::getCluster( SvtxTrack* track, string detName )
0183 {
0184   /* track direction */
0185   float eta = 0;
0186   float phi = 0;
0187
0188   /* closest cluster */
0189   RawCluster* cluster_min_dr = NULL;
0190
0191   // pull the clusters
0192   string clusternodename = "CLUSTER_" + detName;
0193   RawClusterContainer *clusterList = findNode::getClass<RawClusterContainer>(_topNode,clusternodename.c_str());
0194   if (!clusterList) {
0195     cerr << PHWHERE << " ERROR: Can't find node " << clusternodename << endl;
0196     return NULL;
0197   }
0198
0199   // loop over all clusters and find nearest
0200   double min_r = NAN;
0201
0202   for (unsigned int k = 0; k < clusterList->size(); ++k) {
0203
0204     RawCluster *cluster = clusterList->getCluster(k);
0205
0206     double dphi = atan2(sin(phi-cluster->get_phi()),cos(phi-cluster->get_phi()));
0207
0208     double cluster_theta = atan2( cluster->get_r() , cluster->get_z() );
0209     double cluster_eta =  -log(tan(cluster_theta/2.0));
0210
0211     double deta = eta-cluster_eta;
0212     double r = sqrt(pow(dphi,2)+pow(deta,2));
0213
0214     if (r < min_r) {
0215       min_r = r;
0216       cluster_min_dr = cluster;
0217     }
0218   }
0219
0220 //  if(detName == "FEMC") {
0221 //    if(!secondFlag){
0222 //      femc_cldr = min_r;
0223 //      femc_clE = min_e;
0224 //      femc_clN = min_n;
0225 //    }
0226 //    else{
0227 //      femc_cldr2 = min_r;
0228 //      femc_clE2 = min_e;
0229 //      femc_clN2 = min_n;
0230 //    }
0231 //  }
0232 //  else{
0233 //    if(!secondFlag){
0234 //      fhcal_cldr = min_r;
0235 //      fhcal_clE = min_e;
0236 //      fhcal_clN = min_n;
0237 //    }
0238 //    else{
0239 //      fhcal_cldr2 = min_r;
0240 //      fhcal_clE2 = min_e;
0241 //      fhcal_clN2 = min_n;
0242 //    }
0243 //  }
0244 //
0245   return cluster_min_dr;
0246 }
0247
0248
0249
0250 float
0251 TrackProjectionTools::getE33Barrel( string detName, float phi, float eta ){
0252
0253   string towernodename = "TOWER_CALIB_" + detName;
0254   // Grab the towers
0255   RawTowerContainer* towerList = findNode::getClass<RawTowerContainer>(_topNode, towernodename.c_str());
0256   if (!towerList)
0257     {
0258       std::cout << PHWHERE << ": Could not find node " << towernodename.c_str() << std::endl;
0259       return -1;
0260     }
0261   string towergeomnodename = "TOWERGEOM_" + detName;
0262   RawTowerGeomContainer *towergeo = findNode::getClass<RawTowerGeomContainer>(_topNode, towergeomnodename.c_str());
0263   if (! towergeo)
0264     {
0265       cout << PHWHERE << ": Could not find node " << towergeomnodename.c_str() << endl;
0266       return -1;
0267     }
0268
0269   // calculate 3x3 and 5x5 tower energies
0270   int binphi = towergeo->get_phibin(phi);
0271   int bineta = towergeo->get_etabin(eta);
0272
0273   float energy_3x3 = 0.0;
0274   float energy_5x5 = 0.0;
0275
0276   for (int iphi = binphi-2; iphi <= binphi+2; ++iphi) {
0277     for (int ieta = bineta-2; ieta <= bineta+2; ++ieta) {
0278
0279       // wrap around
0280       int wrapphi = iphi;
0281       if (wrapphi < 0) {
0282         wrapphi = towergeo->get_phibins() + wrapphi;
0283       }
0284       if (wrapphi >= towergeo->get_phibins()) {
0285         wrapphi = wrapphi - towergeo->get_phibins();
0286       }
0287
0288       // edges
0289       if (ieta < 0) continue;
0290       if (ieta >= towergeo->get_etabins()) continue;
0291
0292       RawTower* tower = towerList->getTower(ieta,wrapphi);
0293       if (tower) {
0294         energy_5x5 += tower->get_energy();
0295         if (abs(iphi - binphi)<=1 and abs(ieta - bineta)<=1 )
0296           energy_3x3 += tower->get_energy();
0297       }
0298
0299     }
0300   }
0301
0302   return energy_3x3;
0303 }
0304
0305
0306 float
0307 TrackProjectionTools::getE33Forward( string detName, float tkx, float tky )
0308 {
0309   float twr_sum = 0;
0310
0311   string towernodename = "TOWER_CALIB_" + detName;
0312   // Grab the towers
0313   RawTowerContainer* towers = findNode::getClass<RawTowerContainer>(_topNode, towernodename.c_str());
0314   if (!towers)
0315     {
0316       std::cout << PHWHERE << ": Could not find node " << towernodename.c_str() << std::endl;
0317       return -1;
0318     }
0319   string towergeomnodename = "TOWERGEOM_" + detName;
0320   RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(_topNode, towergeomnodename.c_str());
0321   if (! towergeom)
0322     {
0323       cout << PHWHERE << ": Could not find node " << towergeomnodename.c_str() << endl;
0324       return -1;
0325     }
0326
0327   // Locate the central tower
0328   float r_dist = 9999.0;
0329   int twr_j = -1;
0330   int twr_k = -1;
0331   RawTowerDefs::CalorimeterId calo_id_ = RawTowerDefs::convert_name_to_caloid( detName );
0332
0333   RawTowerContainer::ConstRange begin_end  = towers->getTowers();
0334   RawTowerContainer::ConstIterator itr = begin_end.first;
0335   for (; itr != begin_end.second; ++itr)
0336     {
0337       RawTowerDefs::keytype towerid = itr->first;
0338       RawTowerGeom *tgeo = towergeom->get_tower_geometry(towerid);
0339
0340       float x = tgeo->get_center_x();
0341       float y = tgeo->get_center_y();
0342
0343       float temp_rdist = sqrt(pow(tkx-x,2) + pow(tky-y,2)) ;
0344       if(temp_rdist< r_dist){
0345         r_dist = temp_rdist;
0346         twr_j = RawTowerDefs::decode_index1(towerid);
0347         twr_k = RawTowerDefs::decode_index2(towerid);
0348       }
0349
0350       if( (fabs(tkx-x)<(tgeo->get_size_x()/2.0)) &&
0351           (fabs(tky-y)<(tgeo->get_size_y()/2.0)) ) break;
0352
0353     }
0354
0355   // Use the central tower to sum up the 3x3 energy
0356   if(twr_j>=0 && twr_k>=0){
0357     for(int ij = -1; ij <=1; ij++){
0358       for(int ik = -1; ik <=1; ik++){
0359         RawTowerDefs::keytype temp_towerid = RawTowerDefs::encode_towerid( calo_id_ , twr_j + ij , twr_k + ik );
0360         RawTower *rawtower = towers->getTower(temp_towerid);
0361         if(rawtower) twr_sum += rawtower->get_energy();
0362       }
0363     }
0364   }
0365
0366   return twr_sum;
0367 }