ParticleID/RICHAnalysis/RICHEvaluator.C

0001 #include "RICHEvaluator.h"
0002 #include "dualrich_analyzer.h"
0003 #include "TrackProjectorPid.h"
0004 #include "SetupDualRICHAnalyzer.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/PHG4Hit.h>
0014 #include <g4main/PHG4HitContainer.h>
0015 #include <g4main/PHG4Particle.h>
0016 #include <g4main/PHG4VtxPoint.h>
0017 #include <g4main/PHG4TruthInfoContainer.h>
0018
0019 #include <g4hough/SvtxTrackMap.h>
0020 #include <g4hough/SvtxTrack.h>
0021 #include <g4hough/SvtxTrackState.h>
0022
0023 // ROOT includes
0024 #include <TTree.h>
0025 #include <TFile.h>
0026 #include <TString.h>
0027 #include <TMath.h>
0028 #include <TDatabasePDG.h>
0029 #include <TH1.h>
0030
0031 // other C++ includes
0032 #include <cassert>
0033 #include <algorithm>
0034
0035 using namespace std;
0036
0037 RICHEvaluator::RICHEvaluator(std::string tracksname, std::string richname, std::string filename) :
0038   SubsysReco("RICHEvaluator" ),
0039   _ievent(0),
0040   _detector(richname),
0041   _trackmap_name(tracksname),
0042   _refractive_index(1),
0043   _foutname(filename),
0044   _fout_root(nullptr),
0045   _trackproj(nullptr),
0046   _acquire(nullptr),
0047   _pdg(nullptr),
0048   _radius(220.)
0049 {
0050   _richhits_name = "G4HIT_" + _detector;
0051 }
0052
0053
0054 int
0055 RICHEvaluator::Init(PHCompositeNode *topNode)
0056 {
0057   _verbose = false;
0058   _ievent = 0;
0059
0060   /* create output file */
0061   _fout_root = new TFile(_foutname.c_str(), "RECREATE");
0062
0063   /* create output tree */
0064   reset_tree_vars();
0065   init_tree();
0066   init_tree_small();
0067
0068   /* access to PDG databse information */
0069   _pdg = new TDatabasePDG();
0070
0071   /* Throw warning if refractive index is not greater than 1 */
0072   if ( _refractive_index <= 1 )
0073     cout << PHWHERE << " WARNING: Refractive index for radiator volume is " << _refractive_index << endl;
0074
0075   return 0;
0076 }
0077
0078
0079 int
0080 RICHEvaluator::InitRun(PHCompositeNode *topNode)
0081 {
0082   /* create track projector object */
0083   _trackproj = new TrackProjectorPid( topNode );
0084
0085   /* create acquire object */
0086   _acquire = new SetupDualRICHAnalyzer();
0087
0088   return 0;
0089 }
0090
0091
0092 int
0093 RICHEvaluator::process_event(PHCompositeNode *topNode)
0094 {
0095   _ievent ++;
0096
0097
0098   /* Get truth information */
0099   PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0100
0101   /* Get all photon hits in RICH for this event */
0102   PHG4HitContainer* richhits = findNode::getClass<PHG4HitContainer>(topNode,_richhits_name);
0103
0104   /* Get track collection with all tracks in this event */
0105   SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,_trackmap_name);
0106
0107   /* Check if collections found */
0108   if (!truthinfo) {
0109     cout << PHWHERE << "PHG4TruthInfoContainer not found on node tree" << endl;
0110     return Fun4AllReturnCodes::ABORTEVENT;
0111   }
0112   if (!richhits) {
0113     cout << PHWHERE << "PHG4HitContainer not found on node tree" << endl;
0114     return Fun4AllReturnCodes::ABORTEVENT;
0115   }
0116   if (!trackmap) {
0117     cout << PHWHERE << "SvtxTrackMap node not found on node tree" << endl;
0118     return Fun4AllReturnCodes::ABORTEVENT;
0119   }
0120
0121
0122   /* Loop over tracks */
0123   for (SvtxTrackMap::ConstIter track_itr = trackmap->begin(); track_itr != trackmap->end(); track_itr++)
0124     {
0125
0126       bool use_reconstructed_momentum = false;
0127       bool use_truth_momentum = false;
0128       bool use_emission_momentum = true;
0129
0130       bool use_reconstructed_point = false;
0131       bool use_approximate_point = true;
0132
0133
0134       /* Store angles to get RMS value */
0135       TH1F *ch_ang = new TH1F("","",1000,0.0,1.0);
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
0144       /* Fill momv object which is the normalized momentum vector of the track in the RICH (i.e. its direction) */
0145       double momv[3] = {0.,0.,0.};
0146
0147       if (use_reconstructed_momentum) {
0148     /* 'Continue' with next track if RICH projection not found for this track */
0149     if ( ! _trackproj->get_projected_momentum( track_j, momv, TrackProjectorPid::SPHERE, _radius ) )
0150       {
0151         cout << "RICH track projection momentum NOT FOUND; next iteration" << endl;
0152         continue;
0153       }
0154       }
0155       if (use_truth_momentum) {
0156     /* Fill with truth momentum instead of reco */
0157     if ( ! _acquire->get_true_momentum( truthinfo, track_j, momv) )
0158           {
0159             cout << "No truth momentum found for track; next iteration" << endl;
0160             continue;
0161           }
0162       }
0163       if (use_emission_momentum) {
0164         /* Fill with vector constructed from emission points (from truth) */
0165         if ( ! _acquire->get_emission_momentum( truthinfo, richhits, track_j, momv) )
0166           {
0167             cout << "No truth momentum from emission points found for track; next iteration" << endl;
0168             continue;
0169           }
0170       }
0171
0172       double momv_norm = sqrt( momv[0]*momv[0] + momv[1]*momv[1] + momv[2]*momv[2] );
0173       momv[0] /= momv_norm;
0174       momv[1] /= momv_norm;
0175       momv[2] /= momv_norm;
0176
0177
0178       /* Get mean emission point from track in RICH */
0179       double m_emi[3] = {0.,0.,0.};
0180
0181       if (use_reconstructed_point) {
0182     /* 'Continue' with next track if RICH projection not found for this track */
0183     if ( ! _trackproj->get_projected_position( track_j, m_emi, TrackProjectorPid::SPHERE, _radius  ) )
0184       {
0185         cout << "RICH track projection position NOT FOUND; next iteration" << endl;
0186         continue;
0187       }
0188       }
0189       if (use_approximate_point) {
0190         m_emi[0] = ((_radius)/momv[2])*momv[0];
0191         m_emi[1] = ((_radius)/momv[2])*momv[1];
0192         m_emi[2] = ((_radius)/momv[2])*momv[2];
0193       }
0194
0195
0196       /* 'Continue' with next track if track doesn't pass through RICH */
0197       if ( ! _trackproj->is_in_RICH( momv ) )
0198     continue;
0199
0200
0201       /* Calculate truth emission angle and truth mass */
0202       if (truthinfo)
0203     {
0204       _theta_true = calculate_true_emission_angle( truthinfo , track_j , _refractive_index );
0205     }
0206
0207       /* Loop over all G4Hits in container (i.e. RICH photons in event) */
0208       PHG4HitContainer::ConstRange rich_hits_begin_end = richhits->getHits();
0209       PHG4HitContainer::ConstIterator rich_hits_iter;
0210
0211       for (rich_hits_iter = rich_hits_begin_end.first; rich_hits_iter !=  rich_hits_begin_end.second; ++rich_hits_iter)
0212     {
0213       PHG4Hit *hit_i = rich_hits_iter->second;
0214       PHG4Particle* particle = NULL;
0215       PHG4Particle* parent = NULL;
0216       PHG4VtxPoint* vertex = NULL;
0217
0218       if ( truthinfo )
0219         {
0220           particle = truthinfo->GetParticle( hit_i->get_trkid() );
0221           parent = truthinfo->GetParticle( particle->get_parent_id() );
0222           vertex = truthinfo->GetVtx( particle->get_vtx_id() );
0223         }
0224
0225       _hit_x0 =  hit_i->get_x(0);
0226       _hit_y0 =  hit_i->get_y(0);
0227       _hit_z0 =  hit_i->get_z(0);
0228
0229       _emi_x = vertex->get_x();
0230       _emi_y = vertex->get_y();
0231       _emi_z = vertex->get_z();
0232
0233       _track_px = particle->get_px();
0234       _track_py = particle->get_py();
0235       _track_pz = particle->get_pz();
0236
0237       _mtrack_px = parent->get_px();
0238       _mtrack_py = parent->get_py();
0239       _mtrack_pz = parent->get_pz();
0240       _mtrack_ptot = sqrt( _mtrack_px*_mtrack_px + _mtrack_py*_mtrack_py + _mtrack_pz*_mtrack_pz );
0241
0242       _track_e = particle->get_e();
0243       _mtrack_e = parent->get_e();
0244       _edep = hit_i->get_edep();
0245
0246       _bankid = 0;
0247       _volumeid = hit_i->get_detid();
0248       _hitid = hit_i->get_hit_id();
0249       _pid = particle->get_pid();
0250       _mpid = parent->get_pid();
0251       _trackid = particle->get_track_id();
0252       _mtrackid = parent->get_track_id();
0253       _otrackid = track_j->get_id();
0254
0255       /* Set reconstructed emission angle and reconstructed mass for output tree */
0256       _theta_reco = _acquire->calculate_emission_angle( m_emi, momv, hit_i );
0257       ch_ang->Fill(_theta_reco);
0258
0259       /* Fill tree */
0260       _tree_rich->Fill();
0261
0262
0263     } // END loop over photons
0264
0265
0266       _theta_rms = ch_ang->GetRMS();
0267       _theta_mean = ch_ang->GetMean();
0268
0269       /* Fill condensed tree after every track */
0270       _tree_rich_small->Fill();
0271
0272
0273     } // END loop over tracks
0274
0275
0276   return 0;
0277 }
0278
0279
0280 double RICHEvaluator::calculate_true_emission_angle( PHG4TruthInfoContainer* truthinfo, SvtxTrack * track, double index )
0281 {
0282   /* Get truth particle associated with track */
0283   PHG4Particle* particle = truthinfo->GetParticle( track->get_truth_track_id() );
0284
0285   /* Get particle ID */
0286   int pid = particle->get_pid();
0287
0288   /* Get particle mass */
0289   double mass = _pdg->GetParticle(pid)->Mass();
0290
0291   /* Get particle total momentum */
0292   double ptotal = sqrt( particle->get_px() * particle->get_px() +
0293                         particle->get_py() * particle->get_py() +
0294                         particle->get_pz() * particle->get_pz() );
0295
0296   /* Calculate beta = v/c */
0297   double beta = ptotal / sqrt( mass * mass + ptotal * ptotal );
0298
0299   /* Calculate emission angle for Cerenkov light */
0300   double theta_c = acos( 1 / ( index * beta ) );
0301
0302   return theta_c;
0303 }
0304
0305 int
0306 RICHEvaluator::End(PHCompositeNode *topNode)
0307 {
0308   _fout_root->cd();
0309   _tree_rich->Write();
0310   _tree_rich_small->Write();
0311   _fout_root->Close();
0312
0313   return 0;
0314 }
0315
0316
0317 void
0318 RICHEvaluator::reset_tree_vars()
0319 {
0320   _hit_x0 = -999;
0321   _hit_y0 = -999;
0322   _hit_z0 = -999;
0323
0324   _emi_x = -999;
0325   _emi_y = -999;
0326   _emi_z = -999;
0327
0328   _track_px = -999;
0329   _track_py = -999;
0330   _track_pz = -999;
0331
0332   _mtrack_px = -999;
0333   _mtrack_py = -999;
0334   _mtrack_pz = -999;
0335   _mtrack_ptot = -999;
0336
0337   _track_e = -999;
0338   _mtrack_e = -999;
0339   _edep = -999;
0340
0341   _bankid = -999;
0342   _volumeid = -999;
0343   _hitid = -999;
0344   _pid = -999;
0345   _mpid = -999;
0346   _trackid = -999;
0347   _mtrackid = -999;
0348   _otrackid = -999;
0349
0350   _theta_true = -999;
0351   _theta_reco = -999;
0352   _theta_mean = -999;
0353   _theta_rms = -999;
0354
0355   return;
0356 }
0357
0358
0359 int
0360 RICHEvaluator::init_tree()
0361 {
0362   _tree_rich = new TTree("eval_rich","RICH Evaluator info");
0363
0364   _tree_rich->Branch("event", &_ievent, "Event number /I");
0365   _tree_rich->Branch("hit_x", &_hit_x0, "Global x-hit /D");
0366   _tree_rich->Branch("hit_y", &_hit_y0, "Global y-hit /D");
0367   _tree_rich->Branch("hit_z", &_hit_z0, "Global z-hit /D");
0368   _tree_rich->Branch("emi_x", &_emi_x, "Photon emission x-coord. /D");
0369   _tree_rich->Branch("emi_y", &_emi_y, "Photon emission y-coord. /D");
0370   _tree_rich->Branch("emi_z", &_emi_z, "Photon emission z-coord. /D");
0371   _tree_rich->Branch("px", &_track_px, "Particle x-momentum /D");
0372   _tree_rich->Branch("py", &_track_py, "Particle y-momentum /D");
0373   _tree_rich->Branch("pz", &_track_pz, "Particle z-momentum /D");
0374   _tree_rich->Branch("mpx", &_mtrack_px, "Mother x-momentum /D");
0375   _tree_rich->Branch("mpy", &_mtrack_py, "Mother y-momentum /D");
0376   _tree_rich->Branch("mpz", &_mtrack_pz, "Mother z-momentum /D");
0377   _tree_rich->Branch("mptot", &_mtrack_ptot, "Mother total momentum /D");
0378   _tree_rich->Branch("e", &_track_e, "Track energy /D");
0379   _tree_rich->Branch("me", &_mtrack_e, "Mother track energy /D");
0380   _tree_rich->Branch("edep", &_edep, "Energy deposited in material /D");
0381   _tree_rich->Branch("bankid", &_bankid, "Bank ID /I");
0382   _tree_rich->Branch("volumeid", &_volumeid, "Volume ID /I");
0383   _tree_rich->Branch("hitid", &_hitid, "Hit ID /I");
0384   _tree_rich->Branch("pid", &_pid, "Particle ID /I");
0385   _tree_rich->Branch("mpid", &_mpid, "Mother ID /I");
0386   _tree_rich->Branch("trackid", &_trackid, "Track ID /I");
0387   _tree_rich->Branch("mtrackid", &_mtrackid, "Mother track ID /I");
0388   _tree_rich->Branch("otrackid", &_otrackid, "Ordered track ID /I");
0389
0390   _tree_rich->Branch("theta_true", &_theta_true, "True emission angle /D");
0391   _tree_rich->Branch("theta_reco", &_theta_reco, "Reconstructed emission angle /D");
0392
0393   return 0;
0394 }
0395
0396
0397 int
0398 RICHEvaluator::init_tree_small()
0399 {
0400   /* Condensed ROOT tree, 1 entry per track */
0401   _tree_rich_small = new TTree("eval_rich_small","RICH Evaluator info condensed");
0402
0403   _tree_rich_small->Branch("event", &_ievent, "Event number /I");
0404   _tree_rich_small->Branch("otrackid", &_otrackid, "Ordered track ID /I");
0405   _tree_rich_small->Branch("mptot", &_mtrack_ptot, "Total momentum /D");
0406   _tree_rich_small->Branch("theta_true", &_theta_true, "True emission angle /D");
0407   _tree_rich_small->Branch("theta_mean", &_theta_mean, "Reconstructed angle mean /D");
0408   _tree_rich_small->Branch("theta_rms", &_theta_rms, "Reconstructed angle spread /D");
0409
0410   return 0;
0411 }