analysis/EICAnalysis/TrackProjectorPlaneECAL.cc

0001
0002 #include "TrackProjectorPlaneECAL.h"
0003
0004 /* Fun4All includes */
0005 #include <trackbase_historic/SvtxTrack.h>
0006 #include <trackbase_historic/SvtxTrackMap.h>
0007 #include <trackbase_historic/SvtxTrack_FastSim.h>
0008
0009 #include <phool/PHCompositeNode.h>
0010 #include <phool/PHIODataNode.h>
0011
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013
0014 #include <phgeom/PHGeomUtility.h>
0015
0016 #include <trackbase_historic/SvtxTrack.h>
0017 #include <trackbase_historic/SvtxTrackState_v1.h>
0018 //#include <g4hough/PHG4HoughTransform.h>
0019 #include <trackbase_historic/SvtxTrackMap_v1.h>
0020
0021 #include <calobase/RawCluster.h>
0022 #include <calobase/RawTowerDefs.h>
0023
0024 #include <phfield/PHFieldUtility.h>
0025
0026 #include <phgenfit/Fitter.h>
0027 #include <phgenfit/PlanarMeasurement.h>
0028 #include <phgenfit/Track.h>
0029 #include <phgenfit/SpacepointMeasurement.h>
0030
0031 #include <GenFit/RKTrackRep.h>
0032 #include <GenFit/FieldManager.h>
0033
0034
0035
0036 using namespace std;
0037
0038 TrackProjectorPlaneECAL::TrackProjectorPlaneECAL( PHCompositeNode *topNode ) :
0039   _fitter(nullptr)
0040 {
0041
0042   TGeoManager* tgeo_manager = PHGeomUtility::GetTGeoManager(topNode);
0043   PHField * field = PHFieldUtility::GetFieldMapNode(nullptr, topNode);
0044
0045   _fitter = PHGenFit::Fitter::getInstance(tgeo_manager,field,
0046                                           "DafRef",
0047                                           "RKTrackRep", false);
0048
0049   _fitter->set_verbosity(10);
0050
0051   if (!_fitter) {
0052     cerr << "No fitter found at: " << endl;
0053     cerr << PHWHERE << endl;
0054   }
0055 }
0056
0057 bool
0058 TrackProjectorPlaneECAL::get_projected_position(  SvtxTrack_FastSim * track, RawCluster* cluster, double arr_pos[3], const PROJECTION_SURFACE surf, const float surface_par )
0059 {
0060   /* set position components to 0 */
0061   arr_pos[0] = 0;
0062   arr_pos[1] = 0;
0063   arr_pos[2] = 0;
0064
0065   /* project track */
0066   auto state = project_track( track, cluster, surf, surface_par );
0067
0068   /* Set position at extrapolate position */
0069   if ( state )
0070     {
0071       arr_pos[0] = state->get_x();
0072       arr_pos[1] = state->get_y();
0073       arr_pos[2] = state->get_z();
0074       return true;
0075     }
0076
0077   return false;
0078 }
0079
0080 bool
0081 TrackProjectorPlaneECAL::get_projected_momentum(  SvtxTrack_FastSim * track, RawCluster *cluster, double arr_mom[3], const PROJECTION_SURFACE surf, const float surface_par )
0082 {
0083   /* set momentum components to 0 */
0084   arr_mom[0] = 0;
0085   arr_mom[1] = 0;
0086   arr_mom[2] = 0;
0087
0088   /* project track */
0089   auto state = project_track( track, cluster, surf, surface_par );
0090
0091   /* Set momentum at extrapolate position */
0092   if ( state )
0093     {
0094       arr_mom[0] = state->get_px();
0095       arr_mom[1] = state->get_py();
0096       arr_mom[2] = state->get_pz();
0097       return true;
0098     }
0099
0100   return false;
0101 }
0102 // Iterate through all tracks in the SvtxTrackMap. For every track, we iterate among all of its TrackStates. Each
0103 // track state corresponds to a point in space where the track object has been extrapolated to. Currently, we
0104 // extrapolate via coresoftware/simulation/g4simulation/g4trackfastsim/PHG4TrackFastSim.cc. There are four
0105 // detectors in which we perform extrapolation to, the CEMC, FEMC, EEMC, and FHCAL. For every track, we find
0106 // the best extrapolation to the inputted 'RawCluster' object. We then compare the best states and return the
0107 // best extrapolated track.
0108 SvtxTrack_FastSim*
0109 TrackProjectorPlaneECAL::get_best_track( SvtxTrackMap* trackmap, RawCluster* cluster)
0110 {
0111   std::vector< float > distance_from_track_to_cluster;
0112   bool there_is_a_track=false;
0113
0114   float deltaR = -1;
0115
0116   // Iterate all tracks in the trackmap
0117   for(SvtxTrackMap::ConstIter track_itr = trackmap->begin(); track_itr != trackmap->end(); track_itr++)
0118     {
0119       SvtxTrack_FastSim* the_track = dynamic_cast<SvtxTrack_FastSim*>(track_itr->second);
0120       /* Check if the_track is null ptr */
0121       if(the_track == NULL)
0122      {
0123        distance_from_track_to_cluster.push_back(NAN);
0124        continue;
0125      }
0126       else
0127     {
0128       there_is_a_track = true;
0129     }
0130       float distance_from_state_to_cluster = 9990;
0131       // Iterate all track states in the track object
0132       for(SvtxTrack::ConstStateIter state_itr = the_track->begin_states(); state_itr != the_track->end_states(); state_itr++)
0133     {
0134       float distance_from_state_to_cluster_temp = 9999;
0135
0136       SvtxTrackState* the_state = dynamic_cast<SvtxTrackState*>(state_itr->second);
0137       /* Check if the_state is a NULL pointer */
0138       if(the_state==NULL)
0139         {
0140           cout << "State is NULL, skipping..." << endl;
0141           continue;
0142         }
0143       //cout<<the_state->get_x()<< " : " <<the_state->get_y() << " : " <<the_state->get_z() << endl;
0144       distance_from_state_to_cluster_temp = ( sqrt( (cluster->get_x()-the_state->get_x())*(cluster->get_x()-the_state->get_x()) + (cluster->get_y()-the_state->get_y())*(cluster->get_y()-the_state->get_y()) + (cluster->get_z()-the_state->get_z())*(cluster->get_z()-the_state->get_z())));
0145
0146       /*cout << "Cluster : " << cluster->get_x() << " " << cluster->get_y() << " " << cluster->get_z() << endl;
0147       cout << "State : " << the_state->get_x() << " " << the_state->get_y() << " " << the_state->get_z() << endl;
0148       cout << "Name : " << the_state->get_name() << endl;
0149       cout << " " << endl;*/
0150       if(distance_from_state_to_cluster_temp < distance_from_state_to_cluster)
0151         {
0152           distance_from_state_to_cluster = distance_from_state_to_cluster_temp;
0153         }
0154     }
0155       if(distance_from_state_to_cluster!=9990)
0156     {
0157       distance_from_track_to_cluster.push_back(distance_from_state_to_cluster);
0158
0159     }
0160
0161       /* Position vector of extrapolated track */
0162       //double posv[3] = {0.,0.,0.};
0163       //if(!get_projected_position( the_track,cluster, posv, TrackProjectorPlaneECAL::PLANE_CYLINDER, -1))
0164       //    {
0165       //std::cout << "Track Projection Position NOT FOUND; next iteration" << std::endl;
0166       //      distance_from_track_to_cluster.push_back(NAN);
0167       //      continue;
0168       //    }
0169       // else
0170     //{
0171       // distance_from_track_to_cluster.push_back(  sqrt( (cluster->get_x()-posv[0])*(cluster->get_x()-posv[0]) +
0172       //                           (cluster->get_y()-posv[1])*(cluster->get_y()-posv[1]) +
0173       //                           (cluster->get_z()-posv[2])*(cluster->get_z()-posv[2]) )  );
0174       //
0175       //}
0176     }
0177
0178   /* If no track was found, return NULL */
0179   if(distance_from_track_to_cluster.size()==0)
0180     return NULL;
0181   else if(there_is_a_track==false)
0182     return NULL;
0183   /* Find the track with minimum distance */
0184   /* TODO: What if two tracks are within deltaR */
0185   float min_distance = 9999;
0186   int min_distance_index = 0;
0187   for(unsigned i = 0; i<distance_from_track_to_cluster.size(); i++)
0188     {
0189       if(distance_from_track_to_cluster.at(i)<min_distance)
0190     {
0191       min_distance = distance_from_track_to_cluster.at(i);
0192       min_distance_index = i;
0193     }
0194     }
0195   /* Test min_distance with deltaR */
0196   if(deltaR<0)
0197     return dynamic_cast<SvtxTrack_FastSim*>(trackmap->get(min_distance_index));
0198   else if(min_distance<deltaR)
0199     return dynamic_cast<SvtxTrack_FastSim*>(trackmap->get(min_distance_index));
0200   else
0201     return NULL;
0202 }
0203
0204
0205 SvtxTrackState*
0206 TrackProjectorPlaneECAL::get_best_state( SvtxTrack_FastSim* the_track, RawCluster* cluster)
0207 {
0208   /* Check if the_track is null ptr */
0209   float distance_from_state_to_cluster = 9990;
0210   int best_state_index = -1;
0211   int count = 0;
0212   // Iterate all track states in the track object
0213   for(SvtxTrack::ConstStateIter state_itr = the_track->begin_states(); state_itr != the_track->end_states(); state_itr++)
0214     {
0215       float distance_from_state_to_cluster_temp = 9999;
0216
0217       SvtxTrackState* the_state = dynamic_cast<SvtxTrackState*>(state_itr->second);
0218       /* Check if the_state is a NULL pointer */
0219       if(the_state==NULL)
0220     {
0221       cout << "State is NULL, skipping..." << endl;
0222       continue;
0223     }
0224       //cout<<the_state->get_x()<< " : " <<the_state->get_y() << " : " <<the_state->get_z() << endl;
0225       distance_from_state_to_cluster_temp = ( sqrt( (cluster->get_x()-the_state->get_x())*(cluster->get_x()-the_state->get_x()) + (cluster->get_y()-the_state->get_y())*(cluster->get_y()-the_state->get_y()) + (cluster->get_z()-the_state->get_z())*(cluster->get_z()-the_state->get_z())));
0226       if(distance_from_state_to_cluster_temp < distance_from_state_to_cluster)
0227     {
0228       best_state_index = count;
0229       distance_from_state_to_cluster = distance_from_state_to_cluster_temp;
0230     }
0231       count++;
0232     }
0233   if(distance_from_state_to_cluster==-1||best_state_index==-1)
0234     {
0235       cout << "State Map issue, returning NULL" << endl;
0236     }
0237   count = 0;
0238   for(SvtxTrack::ConstStateIter state_itr = the_track->begin_states(); state_itr != the_track->end_states(); state_itr++)
0239     {
0240       if(count!=best_state_index)
0241     {
0242       count++;
0243       continue;
0244     }
0245       else
0246     {
0247       SvtxTrackState* the_state = dynamic_cast<SvtxTrackState*>(state_itr->second);
0248       return the_state;
0249     }
0250     }
0251
0252   return NULL;
0253 }
0254
0255
0256
0257
0258
0259
0260 char
0261 TrackProjectorPlaneECAL::get_detector()
0262 {
0263   switch(detector){
0264   case CEMC:
0265     return 'C';
0266     break;
0267   case EEMC:
0268     return 'E';
0269     break;
0270   case FEMC:
0271     return 'F';
0272     break;
0273   default:
0274     cout << "WARNING: get_detector was unable to find a defined detector" << endl;
0275     return ' ';
0276   }
0277 }
0278
0279 SvtxTrackState*
0280 TrackProjectorPlaneECAL::project_track(  SvtxTrack_FastSim * track, RawCluster* cluster, const PROJECTION_SURFACE surf, const float surface_par )
0281 {
0282   auto last_state_iter = --track->end_states();
0283   SvtxTrackState * trackstate = last_state_iter->second;
0284
0285   if(!trackstate) {
0286     cout << "No state found here!" << endl;
0287     return NULL;
0288   }
0289
0290   int _pid_guess = 11;
0291   auto rep = unique_ptr<genfit::AbsTrackRep> (new genfit::RKTrackRep(_pid_guess));
0292
0293   unique_ptr<genfit::MeasuredStateOnPlane> msop80 = nullptr;
0294
0295   {
0296     cout << track->size_states() << endl;
0297     TVector3 pos(trackstate->get_x(), trackstate->get_y(), trackstate->get_z());
0298     cout << pos(0) << " : " << pos(1) << " : " << pos(2) << endl;
0299     //pos(0) = 0.0;
0300     //pos(1) = 0.0;
0301     //pos(2) = 0.0;
0302     TVector3 mom(trackstate->get_px(), trackstate->get_py(), trackstate->get_pz());
0303     TMatrixDSym cov(6);
0304     for (int i = 0; i < 6; ++i) {
0305       for (int j = 0; j < 6; ++j) {
0306         cov[i][j] = trackstate->get_error(i, j);
0307       }
0308     }
0309     msop80 = unique_ptr<genfit::MeasuredStateOnPlane> (new genfit::MeasuredStateOnPlane(rep.get()));
0310     msop80->setPosMomCov(pos, mom, cov);
0311   }
0312
0313   /* This is where the actual extrapolation of the track to a surface (cylinder, plane, cone, sphere) happens. */
0314   try {
0315
0316     /* 'rep 'is of type AbsTrackRep (see documentation or header for extrapolation function options ) */
0317     if ( surf == SPHERE )
0318       rep->extrapolateToSphere(*msop80, surface_par, TVector3(0,0,0), false, false);
0319
0320     else if ( surf == CYLINDER )
0321       {
0322     rep->extrapolateToCylinder(*msop80, surface_par, TVector3(0,0,0),  TVector3(0,0,1));
0323       }
0324     else if ( surf == PLANE_CYLINDER )
0325       {
0326     // Get position of cluster
0327         TVector3 cluster_pos(cluster->get_x(),cluster->get_y(),cluster->get_z());
0328     // Get detector normal vector
0329     double xvect=0; double yvect = 0; double zvect = 0;
0330     switch(detector){
0331     case CEMC:
0332       xvect = cos(cluster->get_phi());
0333       yvect = sin(cluster->get_phi());
0334       break;
0335     case EEMC:
0336       zvect = 1;
0337       break;
0338     case FEMC:
0339       cluster_pos(0)=0;
0340       cluster_pos(1)=0;
0341       cluster_pos(2)=295;
0342       zvect = -1;
0343       break;
0344     default:
0345       cout << "WARNING: detector variable native to TrackProjectorPlaneECAL not defined" << endl;
0346       break;
0347     }
0348     // Get normal vector to detector
0349     TVector3 cluster_norm(xvect,yvect,zvect);
0350     // Case in which norm vector remained undefined
0351     if(xvect==0&&yvect==0&&zvect==0)
0352       {
0353         cout << "WARNING: Cluster normal vector uninitialized. Aborting track" << endl;
0354         return NULL;
0355       }
0356
0357
0358     // Extrapolate to detector
0359     rep->extrapolateToCylinder(*msop80, 95, TVector3(0,0,0),  TVector3(0,0,1));
0360     //rep->extrapolateToPlane(*msop80, genfit::SharedPlanePtr( new genfit::DetPlane( cluster_pos, cluster_norm ) ), false, false);
0361
0362       }
0363     else if ( surf == PLANEXY )
0364       {
0365     rep->extrapolateToPlane(*msop80, genfit::SharedPlanePtr( new genfit::DetPlane( TVector3(0, 0, surface_par), TVector3(0, 0, 1) ) ), false, false);
0366
0367       }
0368   } catch (...) {
0369     cout << "track extrapolateToXX failed" << endl;
0370     return NULL;
0371   }
0372
0373   /* Create SvtxTrackState object as storage version of the track state
0374    * to pass projected track state information to other member functions. */
0375   SvtxTrackState_v1 *svtx_state = new SvtxTrackState_v1();
0376
0377   svtx_state->set_x( msop80->getPos().X() );
0378   svtx_state->set_y( msop80->getPos().Y() );
0379   svtx_state->set_z( msop80->getPos().Z() );
0380
0381   svtx_state->set_px( msop80->getMom().x() );
0382   svtx_state->set_py( msop80->getMom().y() );
0383   svtx_state->set_pz( msop80->getMom().z() );
0384
0385   for (int i = 0; i < 6; i++) {
0386     for (int j = i; j < 6; j++) {
0387       svtx_state->set_error(i, j, msop80->get6DCov()[i][j]);
0388     }
0389   }
0390
0391   return svtx_state;
0392 }
0393
0394 void
0395 TrackProjectorPlaneECAL::set_detector( char c )
0396 {
0397   switch(c){
0398   case 'C':
0399   case 'c':
0400     detector = CEMC;
0401     break;
0402   case 'E':
0403   case 'e':
0404     detector = EEMC;
0405     break;
0406   case 'F':
0407   case 'f':
0408     detector = FEMC;
0409     break;
0410   default:
0411     cout << "WARNING: set_detector call received unrecognized char " << c << endl;
0412     break;
0413   }
0414 }
0415
0416 char TrackProjectorPlaneECAL::get_detector_from_cluster(RawCluster* cluster)
0417 {
0418   RawCluster::TowerConstIterator rtiter = cluster->get_towers().first;
0419   //caloid = RawTowerDefs::decode_caloid( rtiter->first );
0420   const char *caloid_to_name = RawTowerDefs::convert_caloid_to_name( RawTowerDefs::decode_caloid( rtiter->first ) ).c_str();
0421   if(strcmp(caloid_to_name,"EEMC")==0) return 'E';
0422   else if(strcmp(caloid_to_name,"CEMC")==0) return 'C';
0423   else if(strcmp(caloid_to_name,"FEMC")==0) return 'F';
0424   else
0425     {
0426       cout << "Unrecognized calorimeter. Default to CEMC" << endl;
0427       return 'C';
0428     }
0429 }