ParticleID/RICHAnalysis/TrackProjectorPid.cc

0001 #include "TrackProjectorPid.h"
0002
0003 /* Fun4All includes */
0004 #include <phool/PHCompositeNode.h>
0005 #include <phool/PHIODataNode.h>
0006
0007 #include <fun4all/Fun4AllReturnCodes.h>
0008
0009 #include <phgeom/PHGeomUtility.h>
0010
0011 #include <g4hough/SvtxTrack.h>
0012 #include <g4hough/SvtxTrackState_v1.h>
0013 #include <g4hough/PHG4HoughTransform.h>
0014
0015 #include <phfield/PHFieldUtility.h>
0016
0017 #include <phgenfit/Fitter.h>
0018 #include <phgenfit/PlanarMeasurement.h>
0019 #include <phgenfit/Track.h>
0020 #include <phgenfit/SpacepointMeasurement.h>
0021
0022 #include <GenFit/RKTrackRep.h>
0023 #include <GenFit/FieldManager.h>
0024
0025 using namespace std;
0026
0027 TrackProjectorPid::TrackProjectorPid( PHCompositeNode *topNode ) :
0028   _fitter(nullptr)
0029 {
0030   TGeoManager* tgeo_manager = PHGeomUtility::GetTGeoManager(topNode);
0031
0032   PHField * field = PHFieldUtility::GetFieldMapNode(nullptr, topNode);
0033
0034   _fitter = PHGenFit::Fitter::getInstance(tgeo_manager,field,
0035                                           "DafRef",
0036                                           "RKTrackRep", false);
0037
0038   _fitter->set_verbosity(10);
0039
0040   if (!_fitter) {
0041     cerr << "No fitter found at: " << endl;
0042     cerr << PHWHERE << endl;
0043   }
0044 }
0045
0046 bool
0047 TrackProjectorPid::get_projected_position(  SvtxTrack * track, double arr_pos[3], const PROJECTION_SURFACE surf, const float surface_par )
0048 {
0049   /* set position components to 0 */
0050   arr_pos[0] = 0;
0051   arr_pos[1] = 0;
0052   arr_pos[2] = 0;
0053
0054   /* project track */
0055   auto state = project_track( track, surf, surface_par );
0056
0057   /* Set position at extrapolate position */
0058   if ( state )
0059     {
0060       arr_pos[0] = state->get_x();
0061       arr_pos[1] = state->get_y();
0062       arr_pos[2] = state->get_z();
0063       return true;
0064     }
0065
0066   return false;
0067 }
0068
0069 bool
0070 TrackProjectorPid::get_projected_momentum(  SvtxTrack * track, double arr_mom[3], const PROJECTION_SURFACE surf, const float surface_par )
0071 {
0072   /* set momentum components to 0 */
0073   arr_mom[0] = 0;
0074   arr_mom[1] = 0;
0075   arr_mom[2] = 0;
0076
0077   /* project track */
0078   auto state = project_track( track, surf, surface_par );
0079
0080   /* Set momentum at extrapolate position */
0081   if ( state )
0082     {
0083       arr_mom[0] = state->get_px();
0084       arr_mom[1] = state->get_py();
0085       arr_mom[2] = state->get_pz();
0086       return true;
0087     }
0088
0089   return false;
0090 }
0091
0092 SvtxTrackState*
0093 TrackProjectorPid::project_track(  SvtxTrack * track, const PROJECTION_SURFACE surf, const float surface_par )
0094 {
0095   /* Do projection */
0096   std::vector<double> point;
0097   point.assign(3, -9999.);
0098
0099   auto last_state_iter = --track->end_states();
0100
0101   SvtxTrackState * trackstate = last_state_iter->second;
0102
0103   if(!trackstate) {
0104     cout << "No state found here!" << endl;
0105     return NULL;
0106   }
0107
0108   int _pid_guess = -211;
0109   auto rep = unique_ptr<genfit::AbsTrackRep> (new genfit::RKTrackRep(_pid_guess));
0110
0111   unique_ptr<genfit::MeasuredStateOnPlane> msop80 = nullptr;
0112
0113   {
0114     TVector3 pos(trackstate->get_x(), trackstate->get_y(), trackstate->get_z());
0115     TVector3 mom(trackstate->get_px(), trackstate->get_py(), trackstate->get_pz());
0116
0117     TMatrixDSym cov(6);
0118     for (int i = 0; i < 6; ++i) {
0119       for (int j = 0; j < 6; ++j) {
0120         cov[i][j] = trackstate->get_error(i, j);
0121       }
0122     }
0123
0124     msop80 = unique_ptr<genfit::MeasuredStateOnPlane> (new genfit::MeasuredStateOnPlane(rep.get()));
0125
0126     msop80->setPosMomCov(pos, mom, cov);
0127   }
0128
0129   /* This is where the actual extrapolation of the track to a surface (cylinder, plane, cone, sphere) happens. */
0130   try {
0131
0132     /* 'rep 'is of type AbsTrackRep (see documentation or header for extrapolation function options ) */
0133     if ( surf == SPHERE )
0134       rep->extrapolateToSphere(*msop80, surface_par, TVector3(0,0,0), false, false);
0135
0136     else if ( surf == CYLINDER )
0137       rep->extrapolateToCylinder(*msop80, surface_par, TVector3(0,0,0),  TVector3(0,0,1));
0138
0139
0140     else if ( surf == PLANEXY )
0141       rep->extrapolateToPlane(*msop80, genfit::SharedPlanePtr( new genfit::DetPlane( TVector3(0, 0, surface_par), TVector3(0, 0, 1) ) ), false, false);
0142
0143   } catch (...) {
0144     cout << "track extrapolateToXX failed" << endl;
0145     return NULL;
0146   }
0147
0148   /* Create SvtxTrackState object as storage version of the track state
0149    * to pass projected track state information to other member functions. */
0150   SvtxTrackState_v1 *svtx_state = new SvtxTrackState_v1();
0151
0152   svtx_state->set_x( msop80->getPos().X() );
0153   svtx_state->set_y( msop80->getPos().Y() );
0154   svtx_state->set_z( msop80->getPos().Z() );
0155
0156   svtx_state->set_px( msop80->getMom().x() );
0157   svtx_state->set_py( msop80->getMom().y() );
0158   svtx_state->set_pz( msop80->getMom().z() );
0159
0160   for (int i = 0; i < 6; i++) {
0161     for (int j = i; j < 6; j++) {
0162       svtx_state->set_error(i, j, msop80->get6DCov()[i][j]);
0163     }
0164   }
0165
0166   return svtx_state;
0167 }
0168
0169 bool
0170 TrackProjectorPid::is_in_RICH( double momv[3] )
0171 {
0172   double eta = atanh( momv[2] );
0173
0174   // rough pseudorapidity cut
0175   if (eta > 1.45 && eta < 7.5)
0176     return true;
0177
0178   return false;
0179 }