ParticleID/RICHAnalysis/RICHParticleID.C

0001 #include "RICHParticleID.h"
0002 #include "PidInfo_RICH_v1.h"
0003 #include "PidInfoContainer.h"
0004 #include "TrackProjectorPid.h"
0005 #include "SetupDualRICHAnalyzer.h"
0006 #include "PIDProbabilities.h"
0007
0008 // Fun4All includes
0009 #include <fun4all/Fun4AllServer.h>
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011
0012 #include <phool/getClass.h>
0013 #include <phool/PHCompositeNode.h>
0014
0015 #include <g4main/PHG4Hit.h>
0016 #include <g4main/PHG4HitContainer.h>
0017 #include <g4main/PHG4Particle.h>
0018 #include <g4main/PHG4VtxPoint.h>
0019 #include <g4main/PHG4TruthInfoContainer.h>
0020
0021 #include <g4hough/SvtxTrackMap.h>
0022 #include <g4hough/SvtxTrack.h>
0023 #include <g4hough/SvtxTrackState.h>
0024
0025 // ROOT includes
0026 #include <TTree.h>
0027 #include <TFile.h>
0028 #include <TString.h>
0029 #include <TMath.h>
0030 #include <TH1.h>
0031
0032 // other C++ includes
0033 #include <cassert>
0034 #include <algorithm>
0035
0036 using namespace std;
0037
0038 RICHParticleID::RICHParticleID(std::string tracksname, std::string richname) :
0039   SubsysReco("RICHParticleID" ),
0040   _ievent(0),
0041   _detector(richname),
0042   _trackmap_name(tracksname),
0043   _refractive_index(1),
0044   _pidinfos(nullptr),
0045   _trackproj(nullptr),
0046   _acquire(nullptr),
0047   _particleid(nullptr),
0048   _radius(220.)
0049 {
0050   _richhits_name = "G4HIT_" + _detector;
0051   _pidinfo_node_name = "PIDINFO_" + _detector;
0052 }
0053
0054
0055 int
0056 RICHParticleID::Init(PHCompositeNode *topNode)
0057 {
0058   _verbose = false;
0059   _ievent = 0;
0060
0061   /* create particleid object */
0062   _particleid = new PIDProbabilities();
0063
0064   /* Throw warning if refractive index is not greater than 1 */
0065   if ( _refractive_index <= 1 )
0066     cout << PHWHERE << " WARNING: Refractive index for radiator volume is " << _refractive_index << endl;
0067
0068   return 0;
0069 }
0070
0071
0072 int
0073 RICHParticleID::InitRun(PHCompositeNode *topNode)
0074 {
0075   /* Create nodes for ParticleID */
0076   try
0077     {
0078       CreateNodes(topNode);
0079     }
0080   catch (std::exception &e)
0081     {
0082       std::cout << PHWHERE << ": " << e.what() << std::endl;
0083       throw;
0084     }
0085
0086   /* create track projector object */
0087   _trackproj = new TrackProjectorPid( topNode );
0088
0089   /* create acquire object */
0090   _acquire = new SetupDualRICHAnalyzer();
0091
0092   return 0;
0093 }
0094
0095
0096 int
0097 RICHParticleID::process_event(PHCompositeNode *topNode)
0098 {
0099   _ievent ++;
0100
0101   /* Get truth information (temporary) */
0102   PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0103
0104   /* Get all photon hits in RICH for this event */
0105   PHG4HitContainer* richhits = findNode::getClass<PHG4HitContainer>(topNode,_richhits_name);
0106
0107   /* Get track collection with all tracks in this event */
0108   SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,_trackmap_name);
0109
0110   /* Check if collections found */
0111   if (!truthinfo) {
0112     cout << PHWHERE << "PHG4TruthInfoContainer not found on node tree" << endl;
0113     return Fun4AllReturnCodes::ABORTEVENT;
0114   }
0115   if (!richhits) {
0116     cout << PHWHERE << "PHG4HitContainer not found on node tree" << endl;
0117     return Fun4AllReturnCodes::ABORTEVENT;
0118   }
0119   if (!trackmap) {
0120     cout << PHWHERE << "SvtxTrackMap node not found on node tree" << endl;
0121     return Fun4AllReturnCodes::ABORTEVENT;
0122   }
0123
0124
0125   /* Loop over tracks */
0126   for (SvtxTrackMap::ConstIter track_itr = trackmap->begin(); track_itr != trackmap->end(); track_itr++)
0127     {
0128
0129       bool use_reconstructed_momentum = true;
0130       bool use_truth_momentum = false;
0131       bool use_emission_momentum = false;
0132
0133       bool use_reconstructed_point = true;
0134       bool use_approximate_point = false;
0135
0136
0137       SvtxTrack* track_j = dynamic_cast<SvtxTrack*>(track_itr->second);
0138
0139       /* Check if track_j is a null pointer. */
0140       if (track_j == NULL)
0141         continue;
0142
0143       /* See if there's already a PidInfo object for this track. If yes, retrieve it- if not, create a new one. */
0144       PidInfo *pidinfo_j = _pidinfos->getPidInfo( track_j->get_id() );
0145       if (!pidinfo_j)
0146         {
0147           pidinfo_j = new PidInfo_RICH_v1( track_j->get_id() );
0148           _pidinfos->AddPidInfo( pidinfo_j );
0149     }
0150
0151       /* Fill momv object which is the normalized momentum vector of the track in the RICH (i.e. its direction) */
0152       double momv[3] = {0.,0.,0.};
0153
0154       if (use_reconstructed_momentum) {
0155     /* 'Continue' with next track if RICH projection not found for this track */
0156     if ( ! _trackproj->get_projected_momentum( track_j, momv, TrackProjectorPid::SPHERE, _radius ) )
0157       {
0158         cout << "RICH track projection momentum NOT FOUND; next iteration" << endl;
0159         continue;
0160       }
0161       }
0162       if (use_truth_momentum) {
0163     /* Fill with truth momentum instead of reco */
0164     if ( ! _acquire->get_true_momentum( truthinfo, track_j, momv) )
0165       {
0166         cout << "No truth momentum found for track; next iteration" << endl;
0167         continue;
0168       }
0169       }
0170       if (use_emission_momentum) {
0171     /* Fill with vector constructed from emission points (from truth) */
0172     if ( ! _acquire->get_emission_momentum( truthinfo, richhits, track_j, momv) )
0173       {
0174         cout << "No truth momentum from emission points found for track; next iteration" << endl;
0175         continue;
0176       }
0177       }
0178
0179       double momv_norm = sqrt( momv[0]*momv[0] + momv[1]*momv[1] + momv[2]*momv[2] );
0180       momv[0] /= momv_norm;
0181       momv[1] /= momv_norm;
0182       momv[2] /= momv_norm;
0183
0184
0185       /* Get mean emission point from track in RICH */
0186       double m_emi[3] = {0.,0.,0.};
0187
0188       if (use_reconstructed_point) {
0189     /* 'Continue' with next track if RICH projection not found for this track */
0190     if ( ! _trackproj->get_projected_position( track_j, m_emi, TrackProjectorPid::SPHERE, _radius ) )
0191       {
0192         cout << "RICH track projection position NOT FOUND; next iteration" << endl;
0193         continue;
0194       }
0195       }
0196       if (use_approximate_point) {
0197     m_emi[0] = ((_radius)/momv[2])*momv[0];
0198     m_emi[1] = ((_radius)/momv[2])*momv[1];
0199     m_emi[2] = ((_radius)/momv[2])*momv[2];
0200       }
0201
0202
0203       /* 'Continue' with next track if track doesn't pass through RICH */
0204       if ( ! _trackproj->is_in_RICH( momv ) )
0205     continue;
0206
0207
0208       //cout << "Emission point: " << m_emi[0] << " " << m_emi[1] << " " << m_emi[2] << endl;
0209       //cout << "Momentum vector: " << momv[0] << " " << momv[1] << " " << momv[2] << endl;
0210
0211       /* Vector of reconstructed angles to pass to PID */
0212       vector<float> angles;
0213
0214
0215       /* Loop over all G4Hits in container (RICH photons in event) */
0216       PHG4HitContainer::ConstRange rich_hits_begin_end = richhits->getHits();
0217       PHG4HitContainer::ConstIterator rich_hits_iter;
0218
0219       for (rich_hits_iter = rich_hits_begin_end.first; rich_hits_iter !=  rich_hits_begin_end.second; ++rich_hits_iter)
0220     {
0221
0222       PHG4Hit *hit_i = rich_hits_iter->second;
0223
0224       /* Calculate reconstructed emission angle for output, fill Cherenkov array */
0225       double _theta_reco = _acquire->calculate_emission_angle( m_emi, momv, hit_i );
0226       angles.push_back(_theta_reco);
0227
0228     } // END loop over photons
0229
0230
0231       /* Calculate particle probabilities */
0232       long double probs[4] = {0.,0.,0.,0.};
0233       double momv_magnitude = 0;
0234
0235       /* Emission point momentum only gives a valid unit vector, not magnitude */
0236       if ( use_reconstructed_momentum )
0237     momv_magnitude = momv_norm;
0238       if ( use_truth_momentum )
0239     momv_magnitude = momv_norm;
0240       if ( use_emission_momentum )
0241     {
0242       PHG4Particle* particle = truthinfo->GetParticle( track_j->get_truth_track_id() );
0243       double px = particle->get_px();
0244       double py = particle->get_py();
0245       double pz = particle->get_pz();
0246       momv_magnitude = sqrt( px*px + py*py + pz*pz );
0247     }
0248
0249       if ( !_particleid->particle_probs( angles, momv_magnitude, _refractive_index, probs ) )
0250     {
0251       cout << "No particle ID: ParticleID::particle_probs gives no output" << endl;
0252       continue;
0253     }
0254
0255       /* Set particle probabilities */
0256       pidinfo_j->set_likelihood(PidInfo::ELECTRON, probs[0]);
0257       pidinfo_j->set_likelihood(PidInfo::CHARGEDPION, probs[1]);
0258       pidinfo_j->set_likelihood(PidInfo::CHARGEDKAON, probs[2]);
0259       pidinfo_j->set_likelihood(PidInfo::PROTON, probs[3]);
0260
0261     } // END loop over tracks
0262
0263   return 0;
0264 }
0265
0266 int
0267 RICHParticleID::End(PHCompositeNode *topNode)
0268 {
0269   return 0;
0270 }
0271
0272 void RICHParticleID::CreateNodes(PHCompositeNode *topNode)
0273 {
0274   PHNodeIterator iter(topNode);
0275
0276   PHCompositeNode *dstNode = static_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0277   if (!dstNode)
0278     {
0279       std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0280       throw std::runtime_error("Failed to find DST node in RICHParticleID::CreateNodes");
0281     }
0282
0283   PHNodeIterator dstiter(dstNode);
0284   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode",_detector ));
0285   if(!DetNode){
0286     DetNode = new PHCompositeNode(_detector);
0287     dstNode->addNode(DetNode);
0288   }
0289
0290   _pidinfos = new PidInfoContainer();
0291   PHIODataNode<PHObject> *pidInfoNode = new PHIODataNode<PHObject>(_pidinfos, _pidinfo_node_name.c_str(), "PHObject");
0292   DetNode->addNode(pidInfoNode);
0293 }