ParticleID/RICHAnalysis/SetupDualRICHAnalyzer.cc

0001 // Some documentation, possibly //
0002
0003 #include "TMath.h"
0004
0005 #include "SetupDualRICHAnalyzer.h"
0006 #include "dualrich_analyzer.h"
0007
0008 using namespace std;
0009
0010
0011 SetupDualRICHAnalyzer::SetupDualRICHAnalyzer()
0012 {
0013   _analyzer = new eic_dual_rich();
0014 }
0015
0016 double
0017 SetupDualRICHAnalyzer::calculate_emission_angle( double m_emi[3], double momv[3], PHG4Hit * hit_i )
0018 {
0019   /* Input parameters for indirect ray tracing algorithm */
0020   double Ex = m_emi[0];
0021   double Ey = m_emi[1];
0022   double Ez = m_emi[2];
0023
0024   double Dx = hit_i->get_x(0);
0025   double Dy = hit_i->get_y(0);
0026   double Dz = hit_i->get_z(0);
0027
0028   double vx = momv[0];
0029   double vy = momv[1];
0030   double vz = momv[2];
0031
0032   int sec = hit_i->get_detid();
0033   double cx = -18.5*TMath::Sin(sec*TMath::Pi()/4); // mirror center of each octant
0034   double cy = 18.5*TMath::Cos(sec*TMath::Pi()/4);
0035   double cz = 75;
0036
0037   int select_radiator=0;
0038
0039   /* Set mirror parameters */
0040   double R_mirror = 195; // cm
0041   _analyzer->set_mirror(cx, cy, cz, R_mirror);
0042
0043   /* Call algorithm to determine emission angle of photon i w.r.t. track j */
0044   float theta_c = _analyzer->ind_ray(Ex, Ey, Ez, Dx, Dy, Dz, vx, vy, vz, select_radiator); //Indirect Ray Tracing
0045
0046   return theta_c;
0047 }
0048
0049 bool
0050 SetupDualRICHAnalyzer::get_true_momentum( PHG4TruthInfoContainer* truthinfo, SvtxTrack * track, double arr_mom[3] )
0051 {
0052   arr_mom[0] = 0;
0053   arr_mom[1] = 0;
0054   arr_mom[2] = 0;
0055
0056   /* Get the track's particle, then get truth momentum */
0057   PHG4Particle* particle = truthinfo->GetParticle( track->get_truth_track_id() );
0058   if (particle){
0059     arr_mom[0] = particle->get_px();
0060     arr_mom[1] = particle->get_py();
0061     arr_mom[2] = particle->get_pz();
0062     return true;
0063   }
0064
0065   return false;
0066 }
0067
0068 bool
0069 SetupDualRICHAnalyzer::get_emission_momentum( PHG4TruthInfoContainer* truthinfo, PHG4HitContainer* richhits, SvtxTrack * track, double arr_mom[3] )
0070 {
0071   arr_mom[0] = 0;
0072   arr_mom[1] = 0;
0073   arr_mom[2] = 0;
0074
0075   /* Get photon emission points, then get truth momentum using the vector from first to last emitted photon */
0076   vector<double> emix;
0077   vector<double> emiy;
0078   vector<double> emiz;
0079
0080   /* Loop over all G4Hits in container (i.e. RICH photons in event) */
0081   PHG4HitContainer::ConstRange rich_hits_begin_end = richhits->getHits();
0082   PHG4HitContainer::ConstIterator rich_hits_iter;
0083
0084   for (rich_hits_iter = rich_hits_begin_end.first; rich_hits_iter !=  rich_hits_begin_end.second; ++rich_hits_iter)
0085     {
0086       PHG4Hit *hit_i = rich_hits_iter->second;
0087       PHG4Particle *particle = truthinfo->GetParticle( hit_i->get_trkid() );
0088       PHG4VtxPoint *vertex = truthinfo->GetVtx( particle->get_vtx_id() );
0089
0090       /* Should eventually compare that vertices are along the correct tracks, for now we run single tracks */
0091       emix.push_back( vertex->get_x() );
0092       emiy.push_back( vertex->get_y() );
0093       emiz.push_back( vertex->get_z() );
0094     }
0095
0096
0097   vector<double>::iterator first;
0098   vector<double>::iterator last;
0099   double dx=0;
0100   double dy=0;
0101   double dz=0;
0102
0103   /* Use first-to-last (or first-to-11th) */
0104   first = std::min_element(emiz.begin(),emiz.end());
0105   //last = std::max_element(emiz.begin(),emiz.end());
0106   double p1 = std::distance(emiz.begin(),first);
0107   //double p2 = std::distance(emiz.begin(),last);
0108   if (emiz.size() > p1+11){
0109     dx = emix.at(p1+11) - emix.at(p1);
0110     dy = emiy.at(p1+11) - emiy.at(p1);
0111     dz = emiz.at(p1+11) - emiz.at(p1);
0112   }
0113
0114   /* Fill "momentum" */
0115   arr_mom[0] = dx;
0116   arr_mom[1] = dy;
0117   arr_mom[2] = dz;
0118
0119   return true;
0120 }