ParticleID/RICHAnalysis/FastPid_RICH.C

0001 #include "FastPid_RICH.h"
0002 #include "PidInfo_RICH_v1.h"
0003 #include "PidInfoContainer.h"
0004 #include "TrackProjectorPid.h"
0005
0006 // Fun4All includes
0007 #include <fun4all/Fun4AllServer.h>
0008 #include <fun4all/Fun4AllReturnCodes.h>
0009
0010 #include <phool/getClass.h>
0011 #include <phool/PHCompositeNode.h>
0012
0013 #include <g4main/PHG4Particle.h>
0014 #include <g4main/PHG4TruthInfoContainer.h>
0015
0016 #include <g4hough/SvtxTrackMap.h>
0017 #include <g4hough/SvtxTrack.h>
0018 #include <g4hough/SvtxTrackState.h>
0019
0020 // ROOT includes
0021
0022 // other C++ includes
0023 #include <cassert>
0024 #include <algorithm>
0025
0026 using namespace std;
0027
0028 FastPid_RICH::FastPid_RICH(std::string tracksname, std::string richname) :
0029   SubsysReco("FastPid_RICH" ),
0030   _ievent(0),
0031   _detector(richname),
0032   _trackmap_name(tracksname),
0033   _pidinfos(nullptr),
0034   _trackproj(nullptr),
0035   _radius(220.)
0036 {
0037   _pidinfo_node_name = "PIDINFO_" + _detector;
0038 }
0039
0040
0041 int
0042 FastPid_RICH::Init(PHCompositeNode *topNode)
0043 {
0044   _verbose = false;
0045   _ievent = 0;
0046
0047   return 0;
0048 }
0049
0050
0051 int
0052 FastPid_RICH::InitRun(PHCompositeNode *topNode)
0053 {
0054   /* Create nodes for ParticleID */
0055   try
0056     {
0057       CreateNodes(topNode);
0058     }
0059   catch (std::exception &e)
0060     {
0061       std::cout << PHWHERE << ": " << e.what() << std::endl;
0062       throw;
0063     }
0064
0065   /* create track projector object */
0066   _trackproj = new TrackProjectorPid( topNode );
0067
0068   return 0;
0069 }
0070
0071
0072 int
0073 FastPid_RICH::process_event(PHCompositeNode *topNode)
0074 {
0075   _ievent ++;
0076
0077   /* Get truth information (temporary) */
0078   PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0079
0080   /* Get track collection with all tracks in this event */
0081   SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,_trackmap_name);
0082
0083   /* Check if collections found */
0084   if (!truthinfo) {
0085     cout << PHWHERE << "PHG4TruthInfoContainer not found on node tree" << endl;
0086     return Fun4AllReturnCodes::ABORTEVENT;
0087   }
0088   if (!trackmap) {
0089     cout << PHWHERE << "SvtxTrackMap node not found on node tree" << endl;
0090     return Fun4AllReturnCodes::ABORTEVENT;
0091   }
0092
0093   /* Loop over tracks */
0094   for (SvtxTrackMap::ConstIter track_itr = trackmap->begin(); track_itr != trackmap->end(); track_itr++)
0095     {
0096       SvtxTrack* track_j = dynamic_cast<SvtxTrack*>(track_itr->second);
0097
0098       /* Check if track_j is a null pointer. */
0099       if (track_j == NULL)
0100         continue;
0101
0102       /* See if there's already a PidInfo object for this track. If yes, retrieve it- if not, create a new one. */
0103       PidInfo *pidinfo_j = _pidinfos->getPidInfo( track_j->get_id() );
0104       if (!pidinfo_j)
0105         {
0106           pidinfo_j = new PidInfo_RICH_v1( track_j->get_id() );
0107           _pidinfos->AddPidInfo( pidinfo_j );
0108     }
0109
0110       /* project track to RICH volume and calculate parametrized particle ID response */
0111       /* Get mean emission point from track in RICH volume */
0112       SvtxTrackState *state_j_at_rich = _trackproj->project_track( track_j, TrackProjectorPid::SPHERE, _radius );
0113       if ( ! state_j_at_rich )
0114         {
0115           cout << "RICH track projection position NOT FOUND; next iteration" << endl;
0116       continue;
0117         }
0118
0119       /* Attach track state to PidInfo object */
0120       pidinfo_j->set_track_state( state_j_at_rich );
0121
0122       /* Get truth particle associated with track */
0123       PHG4Particle* particle = truthinfo->GetParticle( track_j->get_truth_track_id() );
0124
0125       /* Get particle ID */
0126       int pid = particle->get_pid();
0127
0128       /* Use parametrized particle ID based on position and momentum at point of track projection */
0129       /* @TODO: Implement parametrized particle ID */
0130       { // beginning of parametrized RICH response
0131
0132     /* approximate RICH acceptance */
0133     float eta_min = 1.45;
0134     float eta_max = 3.5;
0135
0136     /* approximate momentum range for pion / kaon separation */
0137     float p_min = 3.;
0138     float p_max = 50.;
0139
0140     /* check if track projected state in RICH acceptance */
0141     if ( pidinfo_j->get_track_state()->get_eta() > eta_min && pidinfo_j->get_track_state()->get_eta() < eta_max )
0142       {
0143         if ( pidinfo_j->get_track_state()->get_p() > p_min && pidinfo_j->get_track_state()->get_p() < p_max )
0144           {
0145         /* identified kaon */
0146         if ( abs( pid ) == 211 )
0147           pidinfo_j->set_likelihood(PidInfo::CHARGEDPION, 1);
0148
0149         /* identified kaon */
0150         else if ( abs( pid ) == 321 )
0151           pidinfo_j->set_likelihood(PidInfo::CHARGEDKAON, 1);
0152           }
0153       }
0154       } // end of parametrized RICH response
0155
0156       /* print some information to screen */
0157       cout << "True PID: " << pid << endl;
0158       cout << "Position: " << pidinfo_j->get_track_state()->get_x() << ", " << pidinfo_j->get_track_state()->get_y() << ", " << pidinfo_j->get_track_state()->get_z() << endl;
0159       cout << "Momentum: " << pidinfo_j->get_track_state()->get_px() << ", " << pidinfo_j->get_track_state()->get_py() << ", " << pidinfo_j->get_track_state()->get_pz() << endl;
0160       cout << "Eta, |p|: " << pidinfo_j->get_track_state()->get_eta() << ", " << pidinfo_j->get_track_state()->get_p() << endl;
0161       cout << "Likelihood (electron):     " << pidinfo_j->get_likelihood(PidInfo::ELECTRON) << endl;
0162       cout << "Likelihood (charged pion): " << pidinfo_j->get_likelihood(PidInfo::CHARGEDPION) << endl;
0163       cout << "Likelihood (charged kaon): " << pidinfo_j->get_likelihood(PidInfo::CHARGEDKAON) << endl;
0164       cout << "Likelihood (proton):       " << pidinfo_j->get_likelihood(PidInfo::PROTON) << endl;
0165
0166     } // END loop over tracks
0167
0168   return 0;
0169 }
0170
0171 int
0172 FastPid_RICH::End(PHCompositeNode *topNode)
0173 {
0174   return 0;
0175 }
0176
0177 void FastPid_RICH::CreateNodes(PHCompositeNode *topNode)
0178 {
0179   PHNodeIterator iter(topNode);
0180
0181   PHCompositeNode *dstNode = static_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0182   if (!dstNode)
0183     {
0184       std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0185       throw std::runtime_error("Failed to find DST node in RICHParticleID::CreateNodes");
0186     }
0187
0188   PHNodeIterator dstiter(dstNode);
0189   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode",_detector ));
0190   if(!DetNode){
0191     DetNode = new PHCompositeNode(_detector);
0192     dstNode->addNode(DetNode);
0193   }
0194
0195   _pidinfos = new PidInfoContainer();
0196   PHIODataNode<PHObject> *pidInfoNode = new PHIODataNode<PHObject>(_pidinfos, _pidinfo_node_name.c_str(), "PHObject");
0197   DetNode->addNode(pidInfoNode);
0198 }