geometric_acceptance/ResonanceGeometricAcceptance/ResonanceGeometricAcceptance.cc

0001 //____________________________________________________________________________..
0002 //
0003 // This is a template for a Fun4All SubsysReco module with all methods from the
0004 // $OFFLINE_MAIN/include/fun4all/SubsysReco.h baseclass
0005 // You do not have to implement all of them, you can just remove unused methods
0006 // here and in ResonanceGeometricAcceptance.h.
0007 //
0008 // ResonanceGeometricAcceptance(const std::string &name = "ResonanceGeometricAcceptance")
0009 // everything is keyed to ResonanceGeometricAcceptance, duplicate names do work but it makes
0010 // e.g. finding culprits in logs difficult or getting a pointer to the module
0011 // from the command line
0012 //
0013 // ResonanceGeometricAcceptance::~ResonanceGeometricAcceptance()
0014 // this is called when the Fun4AllServer is deleted at the end of running. Be
0015 // mindful what you delete - you do loose ownership of object you put on the node tree
0016 //
0017 // int ResonanceGeometricAcceptance::Init(PHCompositeNode *topNode)
0018 // This method is called when the module is registered with the Fun4AllServer. You
0019 // can create historgrams here or put objects on the node tree but be aware that
0020 // modules which haven't been registered yet did not put antyhing on the node tree
0021 //
0022 // int ResonanceGeometricAcceptance::InitRun(PHCompositeNode *topNode)
0023 // This method is called when the first event is read (or generated). At
0024 // this point the run number is known (which is mainly interesting for raw data
0025 // processing). Also all objects are on the node tree in case your module's action
0026 // depends on what else is around. Last chance to put nodes under the DST Node
0027 // We mix events during readback if branches are added after the first event
0028 //
0029 // int ResonanceGeometricAcceptance::process_event(PHCompositeNode *topNode)
0030 // called for every event. Return codes trigger actions, you find them in
0031 // $OFFLINE_MAIN/include/fun4all/Fun4AllReturnCodes.h
0032 //   everything is good:
0033 //     return Fun4AllReturnCodes::EVENT_OK
0034 //   abort event reconstruction, clear everything and process next event:
0035 //     return Fun4AllReturnCodes::ABORT_EVENT;
0036 //   proceed but do not save this event in output (needs output manager setting):
0037 //     return Fun4AllReturnCodes::DISCARD_EVENT;
0038 //   abort processing:
0039 //     return Fun4AllReturnCodes::ABORT_RUN
0040 // all other integers will lead to an error and abort of processing
0041 //
0042 // int ResonanceGeometricAcceptance::ResetEvent(PHCompositeNode *topNode)
0043 // If you have internal data structures (arrays, stl containers) which needs clearing
0044 // after each event, this is the place to do that. The nodes under the DST node are cleared
0045 // by the framework
0046 //
0047 // int ResonanceGeometricAcceptance::EndRun(const int runnumber)
0048 // This method is called at the end of a run when an event from a new run is
0049 // encountered. Useful when analyzing multiple runs (raw data). Also called at
0050 // the end of processing (before the End() method)
0051 //
0052 // int ResonanceGeometricAcceptance::End(PHCompositeNode *topNode)
0053 // This is called at the end of processing. It needs to be called by the macro
0054 // by Fun4AllServer::End(), so do not forget this in your macro
0055 //
0056 // int ResonanceGeometricAcceptance::Reset(PHCompositeNode *topNode)
0057 // not really used - it is called before the dtor is called
0058 //
0059 // void ResonanceGeometricAcceptance::Print(const std::string &what) const
0060 // Called from the command line - useful to print information when you need it
0061 //
0062 // [[maybe_unused]] suppresses compiler warnings if topNode is not used in this method
0063 //
0064 //____________________________________________________________________________..
0065
0066 #include "ResonanceGeometricAcceptance.h"
0067
0068 #include <fun4all/Fun4AllReturnCodes.h>
0069
0070 #include <phool/PHCompositeNode.h>
0071 #include <phool/getClass.h>
0072
0073 #include <trackbase_historic/SvtxTrackMap.h>
0074 #include <trackbase_historic/SvtxTrack.h>
0075 #include <trackbase_historic/PHG4ParticleSvtxMap.h>
0076 #include <trackbase_historic/SvtxPHG4ParticleMap.h>
0077
0078 #include <g4main/PHG4Particle.h>
0079 #include <g4main/PHG4VtxPoint.h>
0080 #include <g4main/PHG4TruthInfoContainer.h>
0081
0082 #include <KFParticle.h>
0083
0084 #include <TFile.h>
0085 #include <Math/Vector4D.h>
0086 #include <TDatabasePDG.h>
0087 #include <TMatrixD.h>
0088
0089 //____________________________________________________________________________..
0090 ResonanceGeometricAcceptance::ResonanceGeometricAcceptance(const std::string &name):
0091  SubsysReco(name)
0092 {
0093   //std::cout << "ResonanceGeometricAcceptance::ResonanceGeometricAcceptance(const std::string &name) Calling ctor" << std::endl;
0094 }
0095
0096 //____________________________________________________________________________..
0097 ResonanceGeometricAcceptance::~ResonanceGeometricAcceptance()
0098 {
0099   //std::cout << "ResonanceGeometricAcceptance::~ResonanceGeometricAcceptance() Calling dtor" << std::endl;
0100 }
0101
0102 //____________________________________________________________________________..
0103 int ResonanceGeometricAcceptance::Init([[maybe_unused]] PHCompositeNode *topNode)
0104 {
0105   //std::cout << "ResonanceGeometricAcceptance::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0106   return Fun4AllReturnCodes::EVENT_OK;
0107 }
0108
0109 //____________________________________________________________________________..
0110 int ResonanceGeometricAcceptance::InitRun([[maybe_unused]] PHCompositeNode *topNode)
0111 {
0112   //std::cout << "ResonanceGeometricAcceptance::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0113
0114   m_trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0115   m_vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0116   m_truth = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0117   m_truth_reco_map = findNode::getClass<PHG4ParticleSvtxMap>(topNode, "PHG4ParticleSvtxMap");
0118   m_reco_truth_map = findNode::getClass<SvtxPHG4ParticleMap>(topNode, "SvtxPHG4ParticleMap");
0119
0120   if(!m_trackmap)
0121   {
0122     std::cout << PHWHERE << " cannot find SvtxTrackMap named " << m_trackmap_name << std::endl;
0123     return Fun4AllReturnCodes::ABORTEVENT;
0124   }
0125
0126   if(!m_vertexmap)
0127   {
0128     std::cout << PHWHERE << " cannot find SvtxVertexMap named " << m_trackmap_name << std::endl;
0129     return Fun4AllReturnCodes::ABORTEVENT;
0130   }
0131
0132   if(!m_truth)
0133   {
0134     std::cout << PHWHERE << " cannot find PHG4TruthInfoContainer named G4TruthInfo" << std::endl;
0135     return Fun4AllReturnCodes::ABORTEVENT;
0136   }
0137
0138   if(!m_truth_reco_map)
0139   {
0140     std::cout << PHWHERE << " cannot find PHG4ParticleSvtxMap named PHG4ParticleSvtxMap" << std::endl;
0141     return Fun4AllReturnCodes::ABORTEVENT;
0142   }
0143
0144   if(!m_reco_truth_map)
0145   {
0146     std::cout << PHWHERE << " cannot find SvtxPHG4ParticleMap named SvtxPHG4ParticleMap" << std::endl;
0147     return Fun4AllReturnCodes::ABORTEVENT;
0148   }
0149
0150   outfile = new TFile(m_outfile_name.c_str(),"RECREATE");
0151
0152   std::string all_basename = m_mother_name+"_all_candidates";
0153   std::string all_basetitle = "All "+m_mother_name+" candidates in fiducial region with selected decay";
0154   std::string passing_basename = m_mother_name+"_passing_candidates";
0155   std::string passing_basetitle = m_mother_name+" candidates passing geometric acceptance cuts";
0156
0157   all_candidates_vspt = makeHistogram(all_basename,all_basetitle,pt_bins);
0158   all_candidates_vseta = makeHistogram(all_basename,all_basetitle,eta_bins);
0159   all_candidates_vsphi = makeHistogram(all_basename,all_basetitle,phi_bins);
0160   all_candidates_vsrapidity = makeHistogram(all_basename,all_basetitle,rapidity_bins);
0161
0162   passing_candidates_vspt = makeHistogram(passing_basename,passing_basetitle,pt_bins);
0163   passing_candidates_vseta = makeHistogram(passing_basename,passing_basetitle,eta_bins);
0164   passing_candidates_vsphi = makeHistogram(passing_basename,passing_basetitle,phi_bins);
0165   passing_candidates_vsrapidity = makeHistogram(passing_basename,passing_basetitle,rapidity_bins);
0166
0167   return Fun4AllReturnCodes::EVENT_OK;
0168 }
0169
0170 //____________________________________________________________________________..
0171 int ResonanceGeometricAcceptance::process_event([[maybe_unused]] PHCompositeNode *topNode)
0172 {
0173   //std::cout << "ResonanceGeometricAcceptance::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0174
0175   std::map<int,std::vector<PHG4Particle*>> all_daughter_map;
0176   std::map<int,std::vector<PHG4Particle*>> passing_daughter_map;
0177
0178   // find mothers
0179   const auto truth_particle_range = m_truth->GetParticleRange();
0180   for(auto particle_iter = truth_particle_range.first; particle_iter != truth_particle_range.second; ++particle_iter)
0181   {
0182     auto id_and_particle = *particle_iter;
0183     PHG4Particle* particle = id_and_particle.second;
0184
0185     // select mothers with specified PDG ID which are primary
0186     const bool is_mother = (particle->get_pid() == m_mother_pdgid) || (m_include_conjugate && particle->get_pid() == -1*m_mother_pdgid);
0187     bool is_primary;
0188     if(m_mother_pdgid == 310)
0189     {
0190       PHG4Particle* parent = m_truth->GetParticle(particle->get_parent_id());
0191       if(parent)
0192       {
0193         const bool parent_is_K0 = (abs(parent->get_pid()) == 311);
0194         const bool parent_is_primary = (parent->get_track_id() == parent->get_primary_id());
0195         is_primary = (particle->get_track_id() == particle->get_primary_id()) || (parent_is_K0 && parent_is_primary);
0196       }
0197       else
0198       {
0199         is_primary = (particle->get_track_id() == particle->get_primary_id());
0200       }
0201     }
0202     else
0203     {
0204       is_primary = (particle->get_track_id() == particle->get_primary_id());
0205     }
0206
0207     if(is_mother && is_primary)
0208     {
0209       all_daughter_map.insert({particle->get_track_id(),{}});
0210       passing_daughter_map.insert({particle->get_track_id(),{}});
0211     }
0212   }
0213
0214   if(Verbosity()>0)
0215   {
0216     std::cout << "found " << all_daughter_map.size() << " truth mothers" << std::endl;
0217   }
0218
0219   // find daughters and pair with mothers
0220   for(auto particle_iter = truth_particle_range.first; particle_iter != truth_particle_range.second; ++particle_iter)
0221   {
0222     auto id_and_particle = *particle_iter;
0223     PHG4Particle* particle = id_and_particle.second;
0224
0225     if(Verbosity()>5)
0226     {
0227       std::cout << "processing truth track " << particle->get_track_id() << " (PDGID " << particle->get_pid() << ")" << std::endl;
0228     }
0229
0230     PHG4Particle* parent = m_truth->GetParticle(particle->get_parent_id());
0231     if(!parent)
0232     {
0233       if(Verbosity()>5)
0234       {
0235         std::cout << "no parent" << std::endl;
0236       }
0237       continue;
0238     }
0239     if((parent->get_pid() == m_mother_pdgid) || (m_include_conjugate && parent->get_pid() == -1*m_mother_pdgid))
0240     {
0241       if(Verbosity()>5)
0242       {
0243         std::cout << "has mother of PDGID " << m_mother_pdgid << std::endl;
0244       }
0245       const bool is_daughter = std::find(m_daughter_pdgids.begin(),m_daughter_pdgids.end(),particle->get_pid()) != m_daughter_pdgids.end();
0246       const bool is_conjugate_daughter = std::find(m_daughter_pdgids.begin(),m_daughter_pdgids.end(),-1*particle->get_pid()) != m_daughter_pdgids.end();
0247       const bool mother_is_primary = (all_daughter_map.count(parent->get_track_id())>0);
0248       if((is_daughter || (m_include_conjugate && is_conjugate_daughter)) && mother_is_primary)
0249       {
0250         if(Verbosity()>5)
0251         {
0252           std::cout << "is a daughter of a primary mother" << std::endl;
0253         }
0254
0255         all_daughter_map[parent->get_track_id()].push_back(particle);
0256
0257         // get all reco matches for daughter which have silicon hits
0258         std::set<SvtxTrack*> matches_with_si = get_reco_matches_with_silicon_hits(particle);
0259         if(matches_with_si.size()>0)
0260         {
0261           passing_daughter_map[parent->get_track_id()].push_back(particle);
0262         }
0263       }
0264     }
0265 /*
0266     if(parent && get_reco_track(particle))
0267     {
0268       SvtxTrack* reco_track = get_reco_track(particle);
0269       std::cout << "truth track " << particle->get_track_id() << " paired with reco track " << reco_track->get_id() << " has parent PDGID " << parent->get_pid() << std::endl;
0270     }
0271 */
0272   }
0273
0274   if(Verbosity()>0)
0275   {
0276     std::cout << "passing daughter pass complete" << std::endl;
0277   }
0278
0279   // fill denominator with primary mothers which decayed in the correct channel and satisfy fiducial cuts given by HistogramInfo variables
0280   for(const auto& [mother_track_id, daughters] : all_daughter_map)
0281   {
0282     PHG4Particle* mother = m_truth->GetParticle(mother_track_id);
0283     ROOT::Math::PxPyPzEVector mother_lorentzvector(mother->get_px(),mother->get_py(),mother->get_pz(),mother->get_e());
0284
0285     PHG4VtxPoint* gvertex = m_truth->GetVtx(mother->get_vtx_id());
0286
0287     bool correct_decay_channel = has_all_daughters(daughters);
0288     bool within_zvertex_limits = (gvertex->get_z()>=zvertex_limits.first && gvertex->get_z()<=zvertex_limits.second);
0289
0290     bool within_pt_limits = (mother_lorentzvector.Pt()>=pt_limits.first && mother_lorentzvector.Pt()<=pt_limits.second);
0291     bool within_eta_limits = (mother_lorentzvector.Eta()>=eta_limits.first && mother_lorentzvector.Eta()<=eta_limits.second);
0292     bool within_phi_limits = (mother_lorentzvector.Phi()>=phi_limits.first && mother_lorentzvector.Phi()<=phi_limits.second);
0293     bool within_rapidity_limits = (mother_lorentzvector.Rapidity()>=rapidity_limits.first && mother_lorentzvector.Rapidity()<=rapidity_limits.second);
0294
0295     if(correct_decay_channel && within_zvertex_limits && (within_pt_limits && within_eta_limits && within_phi_limits && within_rapidity_limits))
0296     {
0297       all_candidates_vspt->Fill(mother_lorentzvector.Pt());
0298       all_candidates_vseta->Fill(mother_lorentzvector.Eta());
0299       all_candidates_vsphi->Fill(mother_lorentzvector.Phi());
0300       all_candidates_vsrapidity->Fill(mother_lorentzvector.Rapidity());
0301     }
0302   }
0303
0304   if(Verbosity()>0)
0305   {
0306     std::cout << "filled " << m_mother_name << " truth histograms" << std::endl;
0307     std::cout << "current nEntries: " << all_candidates_vspt->GetEntries() << std::endl;
0308   }
0309
0310   // If reco vertex exists, fill numerator with primary mothers whose daughters are reconstructible as KFParticle SVs within the fiducial volume.
0311   // This means that they must satisfy all of the following:
0312   // - each daughter must have a reconstructed track
0313   // - the vertex charge must be correct
0314   // - the PID hypotheses available to KFParticle must match with the reco-charge-modified PDG IDs of the daughters
0315   // - the mother crossing must match all the daughter crossings.
0316   // In addition, the mother momentum produced by KFParticle's internal Kalman filter can significantly differ from that obtained from the raw sum of four-vectors;
0317   // therefore, the fiducial volume cuts must be evaluated relative to the KFParticle mother momentum.
0318   if(!m_vertexmap->empty())
0319   {
0320     for(const auto& [mother_track_id, daughters] : passing_daughter_map)
0321     {
0322       bool correct_decay_channel = has_all_daughters(daughters);
0323       if(!correct_decay_channel)
0324       {
0325         continue;
0326       }
0327       std::vector<std::vector<SvtxTrack*>> reco_combinations = get_all_reco_combinations(daughters);
0328
0329       for(std::vector<SvtxTrack*> combination : reco_combinations)
0330       {
0331         bool charge_swaps_ok = check_charge_swaps(daughters,combination);
0332         bool crossings_match = check_crossings_match(combination);
0333
0334         if(correct_decay_channel && charge_swaps_ok && crossings_match)
0335         {
0336           std::vector<KFParticle> kfp_reco_daughters;
0337           for(size_t itr=0; itr<combination.size(); itr++)
0338           {
0339             kfp_reco_daughters.push_back(make_KFParticle(daughters[itr],combination[itr]));
0340           }
0341
0342           // build mother KFParticle
0343
0344           KFParticle kfp_reco_mother;
0345           kfp_reco_mother.SetConstructMethod(2);
0346           for(KFParticle& kfp_daughter : kfp_reco_daughters)
0347           {
0348             kfp_reco_mother.AddDaughter(kfp_daughter);
0349           }
0350 /*
0351           for(KFParticle& kfp_daughter : kfp_reco_daughters)
0352           {
0353             kfp_daughter.SetProductionVertex(kfp_reco_mother);
0354           }
0355 */
0356           kfp_reco_mother.SetPDG(get_mother_PDGID(daughters));
0357           float mother_mass, mother_mass_err;
0358           kfp_reco_mother.GetMass(mother_mass,mother_mass_err);
0359           float mother_pt, mother_pt_err;
0360           kfp_reco_mother.GetPt(mother_pt,mother_pt_err);
0361
0362           // associate mother with best possible reconstructed vertex
0363
0364           int best_vertex_id = get_best_reco_vertex(kfp_reco_mother,combination[0]->get_crossing());
0365           if(Verbosity()>5)
0366           {
0367             std::cout << "best_vertex_id: " << best_vertex_id << std::endl;
0368           }
0369           if(best_vertex_id == -1)
0370           {
0371             continue;
0372           }
0373
0374           SvtxVertex* pv = m_vertexmap->get(best_vertex_id);
0375           KFParticle kfpv = make_KFParticle_Vertex(pv);
0376           //kfp_reco_mother.SetProductionVertex(kfpv);
0377           //kfp_reco_mother.TransportToDecayVertex();
0378
0379           if(Verbosity()>5)
0380           {
0381             std::cout << "reco vertex z: " << kfpv.GetZ() << std::endl;
0382           }
0383
0384           // get mother kinematic variables, compare with fiducial cuts
0385
0386           float reco_mother_pt, reco_mother_pterr;
0387           kfp_reco_mother.GetPt(reco_mother_pt,reco_mother_pterr);
0388           float reco_mother_eta, reco_mother_etaerr;
0389           kfp_reco_mother.GetEta(reco_mother_eta,reco_mother_etaerr);
0390           float reco_mother_phi, reco_mother_phierr;
0391           kfp_reco_mother.GetPhi(reco_mother_phi,reco_mother_phierr);
0392           float reco_mother_rapidity = kfp_reco_mother.GetRapidity();
0393
0394           bool daughter_pt_above_min = std::all_of(combination.begin(),combination.end(),
0395             [&](SvtxTrack* track){ return track->get_pt() >= m_daughter_min_pT; });
0396
0397           bool within_zvertex_limits = (kfpv.GetZ()>=zvertex_limits.first && kfpv.GetZ()<=zvertex_limits.second);
0398
0399           bool within_pt_limits = (reco_mother_pt>=pt_limits.first && reco_mother_pt<=pt_limits.second);
0400           bool within_eta_limits = (reco_mother_eta>=eta_limits.first && reco_mother_eta<=eta_limits.second);
0401           bool within_phi_limits = (reco_mother_phi>=phi_limits.first && reco_mother_phi<=phi_limits.second);
0402           bool within_rapidity_limits = (reco_mother_rapidity>=rapidity_limits.first && reco_mother_rapidity<=rapidity_limits.second);
0403
0404           if(daughter_pt_above_min && within_zvertex_limits && within_pt_limits && within_eta_limits && within_phi_limits && within_rapidity_limits)
0405           {
0406             if(Verbosity()>1)
0407             {
0408               identify(kfp_reco_mother);
0409               for(KFParticle& kfpd : kfp_reco_daughters) identify(kfpd);
0410             }
0411
0412             // if there are multiple reco combinations for a single truth candidate, scale their contribution down by the number of combinations
0413             float fill_weight = 1./reco_combinations.size();
0414
0415             passing_candidates_vspt->Fill(reco_mother_pt,fill_weight);
0416             passing_candidates_vseta->Fill(reco_mother_eta,fill_weight);
0417             passing_candidates_vsphi->Fill(reco_mother_phi,fill_weight);
0418             passing_candidates_vsrapidity->Fill(reco_mother_rapidity,fill_weight);
0419           }
0420         }
0421       }
0422     }
0423   }
0424
0425   if(Verbosity()>0)
0426   {
0427     std::cout << "filled " << m_mother_name << " reco histograms" << std::endl;
0428     std::cout << "current nEntries: " << passing_candidates_vspt->GetEntries() << std::endl;
0429   }
0430
0431   all_daughter_map.clear();
0432   passing_daughter_map.clear();
0433
0434   return Fun4AllReturnCodes::EVENT_OK;
0435 }
0436
0437 int ResonanceGeometricAcceptance::get_nmaps(TrackSeed* si_seed) const
0438 {
0439   int nmaps = 0;
0440   for(auto c_it = si_seed->begin_cluster_keys(); c_it != si_seed->end_cluster_keys(); ++c_it)
0441   {
0442     TrkrDefs::cluskey ckey = *c_it;
0443     if(TrkrDefs::getTrkrId(ckey) == TrkrDefs::TrkrId::mvtxId)
0444     {
0445       nmaps++;
0446     }
0447   }
0448   return nmaps;
0449 }
0450
0451 std::vector<std::vector<SvtxTrack*>> ResonanceGeometricAcceptance::get_all_reco_combinations(const std::vector<PHG4Particle*>& daughters) const
0452 {
0453   std::vector<std::vector<SvtxTrack*>> reco_tracks;
0454   for(PHG4Particle* daughter : daughters)
0455   {
0456     std::set<SvtxTrack*> reco_track_set = get_reco_matches_with_silicon_hits(daughter);
0457     if(reco_track_set.empty()) return {};
0458     std::vector<SvtxTrack*> reco_track_vector;
0459     reco_track_vector.insert(reco_track_vector.end(),reco_track_set.begin(),reco_track_set.end());
0460     reco_tracks.push_back(reco_track_vector);
0461   }
0462
0463   if(Verbosity()>5)
0464   {
0465     std::cout << "reco tracks for each daughter:" << std::endl;
0466     std::cout << "{";
0467     for(size_t i=0;i<reco_tracks.size();i++)
0468     {
0469       std::cout << "{";
0470       for(size_t j=0;j<reco_tracks[i].size();j++)
0471       {
0472         std::cout << reco_tracks[i][j]->get_id();
0473         if(j<reco_tracks[i].size()-1) std::cout << ", ";
0474       }
0475       std::cout << "}";
0476       if(i<reco_tracks.size()-1) std::cout << ", ";
0477     }
0478     std::cout << "}" << std::endl;
0479   }
0480
0481   std::vector<std::vector<SvtxTrack*>> reco_combinations;
0482   if(reco_tracks.size()==0) return {};
0483   for(SvtxTrack* track : reco_tracks[0])
0484   {
0485     reco_combinations.push_back({track});
0486   }
0487   if(reco_tracks.size()>1)
0488   {
0489     for(size_t iv = 1; iv<reco_tracks.size(); iv++)
0490     {
0491       std::vector<std::vector<SvtxTrack*>> reco_iteration;
0492       for(size_t icomb = 0; icomb<reco_combinations.size(); icomb++)
0493       {
0494         for(size_t itr = 0; itr<reco_tracks[iv].size(); itr++)
0495         {
0496           std::vector<SvtxTrack*> reco_combination = reco_combinations[icomb];
0497           reco_combination.push_back(reco_tracks[iv][itr]);
0498           reco_iteration.push_back(reco_combination);
0499         }
0500       }
0501       reco_combinations = reco_iteration;
0502     }
0503   }
0504
0505   if(Verbosity()>5)
0506   {
0507     std::cout << "reco combinations:" << std::endl;
0508     std::cout << "{";
0509     for(size_t i=0;i<reco_combinations.size();i++)
0510     {
0511       std::cout << "{";
0512       for(size_t j=0;j<reco_combinations[i].size();j++)
0513       {
0514         std::cout << reco_combinations[i][j]->get_id();
0515         if(j<reco_combinations[i].size()-1) std::cout << ", ";
0516       }
0517       std::cout << "}";
0518       if(i<reco_combinations.size()-1) std::cout << ", ";
0519     }
0520     std::cout << "}" << std::endl;
0521   }
0522   return reco_combinations;
0523 }
0524
0525 int ResonanceGeometricAcceptance::get_mother_PDGID(const std::vector<PHG4Particle*>& daughters) const
0526 {
0527   bool all_conjugate_daughters = std::all_of(m_daughter_pdgids.begin(),m_daughter_pdgids.end(),[&](int pdgid){ return has_daughter(-1*pdgid,daughters); });
0528   if(all_conjugate_daughters)
0529   {
0530     return -1*m_mother_pdgid;
0531   }
0532   else
0533   {
0534     return m_mother_pdgid;
0535   }
0536 }
0537
0538 void ResonanceGeometricAcceptance::identify(KFParticle particle) const
0539 {
0540   std::cout << "-----------------------------------------------------------------------" << std::endl;
0541   std::cout << "Track ID: " << particle.Id() << std::endl;
0542   std::cout << "PDG ID: " << particle.GetPDG() << ", charge: " << (int) particle.GetQ() << ", mass: " << particle.GetMass() << " GeV" << std::endl;
0543   std::cout << "(px,py,pz) = (" << particle.GetPx() << " +/- " << std::sqrt(particle.GetCovariance(3, 3)) << ", ";
0544   std::cout << particle.GetPy() << " +/- " << std::sqrt(particle.GetCovariance(4, 4)) << ", ";
0545   std::cout << particle.GetPz() << " +/- " << std::sqrt(particle.GetCovariance(5, 5)) << ") GeV" << std::endl;
0546   std::cout << "(x,y,z) = (" << particle.GetX() << " +/- " << std::sqrt(particle.GetCovariance(0, 0)) << ", ";
0547   std::cout << particle.GetY() << " +/- " << std::sqrt(particle.GetCovariance(1, 1)) << ", ";
0548   std::cout << particle.GetZ() << " +/- " << std::sqrt(particle.GetCovariance(2, 2)) << ") cm\n"
0549             << std::endl;
0550   std::cout << "-----------------------------------------------------------------------" << std::endl;
0551 }
0552
0553 int ResonanceGeometricAcceptance::get_best_reco_vertex(KFParticle mother, const int crossing) const
0554 {
0555   int best_masserr_vertex_id = -1;
0556   float best_mass_err = FLT_MAX;
0557
0558   int best_dcachi2_vertex_id = -1;
0559   float best_dcachi2 = FLT_MAX;
0560
0561   for(auto vtx_iter = m_vertexmap->begin(); vtx_iter != m_vertexmap->end(); ++vtx_iter)
0562   {
0563     SvtxVertex* vertex = (*vtx_iter).second;
0564     if(vertex->get_beam_crossing() != crossing)
0565     {
0566       continue;
0567     }
0568     KFParticle mother_copy = mother;
0569     KFParticle kfv = make_KFParticle_Vertex(vertex);
0570     mother_copy.SetProductionVertex(kfv);
0571
0572     if(mother_copy.GetErrMass() < best_mass_err)
0573     {
0574       best_masserr_vertex_id = vertex->get_id();
0575       best_mass_err = mother_copy.GetErrMass();
0576     }
0577
0578     if(mother_copy.GetDeviationFromVertex(kfv) < best_dcachi2)
0579     {
0580       best_dcachi2_vertex_id = vertex->get_id();
0581       best_dcachi2 = mother_copy.GetDeviationFromVertex(kfv);
0582     }
0583   }
0584
0585   if(Verbosity()>5)
0586   {
0587     std::cout << "best vertex ID by DCA chi2: " << best_dcachi2_vertex_id << ", by mass err: " << best_masserr_vertex_id << std::endl;
0588   }
0589
0590   return best_dcachi2_vertex_id;
0591 }
0592
0593 KFParticle ResonanceGeometricAcceptance::make_KFParticle_Vertex(SvtxVertex* vertex) const
0594 {
0595   double parameters[6] = {vertex->get_x(), vertex->get_y(), vertex->get_z()};
0596   double covariance[21];
0597   for(int i=0;i<21;i++)
0598   {
0599     covariance[i] = 0.;
0600   }
0601   int index = 0;
0602   for(int i=0;i<3;i++)
0603   {
0604     for(int j=0;j<=i;j++)
0605     {
0606       covariance[index] = vertex->get_error(i,j);
0607       index++;
0608     }
0609   }
0610
0611   KFParticle kfv;
0612   kfv.Create(parameters,covariance,0,-1);
0613   kfv.Chi2() = vertex->get_chisq();
0614   kfv.NDF() = vertex->get_ndof();
0615
0616   return kfv;
0617 }
0618
0619 KFParticle ResonanceGeometricAcceptance::make_KFParticle(PHG4Particle* particle, SvtxTrack* track) const
0620 {
0621   double mass = TDatabasePDG::Instance()->GetParticle(particle->get_pid())->Mass();
0622
0623   double parameters[6] = {track->get_x(), track->get_y(), track->get_z(), track->get_px(), track->get_py(), track->get_pz()};
0624   double covariance[21];
0625   int index = 0;
0626   for(int i=0;i<6;i++)
0627   {
0628     for(int j=0;j<=i;j++)
0629     {
0630       covariance[index] = track->get_error(i,j);
0631       index++;
0632     }
0633   }
0634   double chi2 = track->get_chisq();
0635   double ndf = track->get_ndf();
0636   int charge = track->get_charge();
0637
0638   KFParticle kfp;
0639   kfp.SetConstructMethod(2);
0640   kfp.Create(parameters,covariance,charge,mass);
0641   kfp.Chi2() = chi2;
0642   kfp.NDF() = ndf;
0643   kfp.SetPDG(particle->get_pid());
0644   kfp.SetId(track->get_id());
0645
0646   return kfp;
0647 }
0648 /*
0649 bool ResonanceGeometricAcceptance::truth_track_is_best_match_to_reco_track(const PHG4Particle* particle) const
0650 {
0651   SvtxTrack* reco_track = get_reco_track(particle);
0652   PHG4Particle* truth_track = nullptr;
0653   std::map<float, std::set<int>> truth_set = m_reco_truth_map->get(reco_track->get_id());
0654   std::cout << "getting truth track ID for reco track " << reco_track->get_id() << std::endl;
0655   for(auto [weight, trackset] : truth_set)
0656   {
0657     std::cout << "reco track " << reco_track->get_id() << " pairs with weight " << weight << " on the following truth tracks:" << std::endl;
0658     for(int trackID : trackset)
0659     {
0660       std::cout << trackID << ", ";
0661     }
0662     std::cout << std::endl;
0663   }
0664   if (!truth_set.empty())
0665   {
0666     std::pair<float, std::set<int>> best_weight = *truth_set.rbegin();
0667     int best_truth_id = *best_weight.second.rbegin();
0668     truth_track = m_truth->GetParticle(best_truth_id);
0669   }
0670   std::cout << "reco track " << reco_track->get_id() << " pairs with truth track " << truth_track->get_track_id() << std::endl;
0671   return (truth_track->get_track_id() == particle->get_track_id());
0672 }
0673 */
0674 std::set<SvtxTrack*> ResonanceGeometricAcceptance::get_reco_matches_with_silicon_hits(const PHG4Particle* particle) const
0675 {
0676   std::set<SvtxTrack*> reco_matches = get_best_reco_matches(particle);
0677
0678   std::set<SvtxTrack*> reco_matches_with_silicon_hits;
0679
0680   for(SvtxTrack* reco_match : reco_matches)
0681   {
0682     TrackSeed* si_seed = reco_match->get_silicon_seed();
0683     if(si_seed && get_nmaps(si_seed)>0)
0684     {
0685       if(Verbosity()>5)
0686       {
0687         std::cout << "matched reco track with ID " << reco_match->get_id() << " has " << si_seed->size_cluster_keys() << " silicon hits, with " << get_nmaps(si_seed) << "MVTX hits" << std::endl;
0688       }
0689       reco_matches_with_silicon_hits.insert(reco_match);
0690     }
0691   }
0692
0693   return reco_matches_with_silicon_hits;
0694 }
0695
0696 std::set<SvtxTrack*> ResonanceGeometricAcceptance::get_best_reco_matches(const PHG4Particle* particle) const
0697 {
0698   std::set<SvtxTrack*> reco_matches;
0699   std::map<float, std::set<unsigned int>> reco_set = m_truth_reco_map->get(particle->get_track_id());
0700   if(Verbosity()>5)
0701   {
0702     std::cout << "getting reco track ID for truth track " << particle->get_track_id() << " (PDGID " << particle->get_pid() << ")" << std::endl;
0703   }
0704   for(auto [weight, trackset] : reco_set)
0705   {
0706     if(Verbosity()>5)
0707     {
0708       std::cout << "truth track " << particle->get_track_id() << " pairs with weight " << weight << " on the following reco tracks:" << std::endl;
0709       for(unsigned int trackID : trackset)
0710       {
0711         std::cout << trackID << ", ";
0712       }
0713       std::cout << std::endl;
0714     }
0715   }
0716   if (!reco_set.empty())
0717   {
0718     for(const auto& [weight, track_set] : reco_set)
0719     {
0720       if(weight>m_truth_weight_threshold)
0721       {
0722         for(unsigned int trackID : track_set)
0723         {
0724           reco_matches.insert(m_trackmap->get(trackID));
0725         }
0726       }
0727     }
0728   }
0729   return reco_matches;
0730 }
0731
0732 bool ResonanceGeometricAcceptance::has_all_daughters(const std::vector<PHG4Particle*>& daughters) const
0733 {
0734   bool all_daughters = std::all_of(m_daughter_pdgids.begin(),m_daughter_pdgids.end(),[&](int pdgid){ return has_daughter(pdgid,daughters); });
0735   bool all_conjugate_daughters = std::all_of(m_daughter_pdgids.begin(),m_daughter_pdgids.end(),[&](int pdgid){ return has_daughter(-1*pdgid,daughters); });
0736   return (all_daughters || (m_include_conjugate && all_conjugate_daughters));
0737 }
0738
0739 bool ResonanceGeometricAcceptance::has_daughter(const int pdgid, const std::vector<PHG4Particle*>& daughters) const
0740 {
0741   return (std::find_if(daughters.begin(),daughters.end(),
0742     [&](PHG4Particle* daughter){ return (daughter->get_pid() == pdgid); })
0743     != daughters.end());
0744 }
0745
0746 bool ResonanceGeometricAcceptance::check_charge_swaps(const std::vector<PHG4Particle*>& daughters, const std::vector<SvtxTrack*>& reco_daughters) const
0747 {
0748   if(Verbosity()>5)
0749   {
0750     std::cout << "charge swap check:" << std::endl;
0751   }
0752   // get the set of "effective" PDG IDs after any charge swaps
0753   if(daughters.size() != reco_daughters.size())
0754   {
0755     std::cout << "ERROR: truth and reco daughters have different number of tracks! Something has gone very wrong here..." << std::endl;
0756     return false;
0757   }
0758   std::vector<int> reco_daughter_PDGIDs;
0759   for(size_t i=0; i<daughters.size(); i++)
0760   {
0761     int gflavor = daughters[i]->get_pid();
0762     SvtxTrack* reco_track = reco_daughters[i];
0763     const int reco_charge = reco_track->get_charge();
0764     if(Verbosity()>5)
0765     {
0766       std::cout << "daughter_gflavor: " << gflavor << ", charge: " << reco_charge << std::endl;
0767     }
0768     const int truth_charge = gflavor > 0 ? 1 : -1;
0769
0770     // charges match if their product is 1
0771     if(reco_charge*truth_charge == 1)
0772     {
0773       reco_daughter_PDGIDs.push_back(gflavor);
0774     }
0775     else
0776     {
0777       reco_daughter_PDGIDs.push_back(-1*gflavor);
0778     }
0779   }
0780
0781   // check to see that we still have the right mix of flavors to form a secondary vertex
0782   bool has_all_daughters = std::all_of(m_daughter_pdgids.begin(),m_daughter_pdgids.end(),
0783     [&](int flavor){ return std::find(reco_daughter_PDGIDs.begin(),reco_daughter_PDGIDs.end(),flavor) != reco_daughter_PDGIDs.end(); });
0784   bool has_all_opposite_daughters = std::all_of(m_daughter_pdgids.begin(),m_daughter_pdgids.end(),
0785     [&](int flavor){ return std::find(reco_daughter_PDGIDs.begin(),reco_daughter_PDGIDs.end(),-1*flavor) != reco_daughter_PDGIDs.end(); });
0786
0787   if(has_all_daughters && Verbosity()>5)
0788   {
0789     std::cout << "has all daughters" << std::endl;
0790   }
0791   if(has_all_opposite_daughters && Verbosity()>5)
0792   {
0793     std::cout << "has all opposite daughters" << std::endl;
0794   }
0795
0796   bool charge_swap_ok = (has_all_daughters || (m_include_conjugate && has_all_opposite_daughters));
0797   if(Verbosity()>5)
0798   {
0799     std::cout << "charge_swap_ok: " << charge_swap_ok << std::endl;
0800   }
0801   return charge_swap_ok;
0802 }
0803
0804 bool ResonanceGeometricAcceptance::check_crossings_match(const std::vector<SvtxTrack*>& reco_daughters) const
0805 {
0806   if(Verbosity()>5)
0807   {
0808     std::cout << "crossing check:" << std::endl;
0809     for(SvtxTrack* reco_daughter : reco_daughters)
0810     {
0811       std::cout << "daughter crossing: " << reco_daughter->get_crossing() << std::endl;
0812     }
0813   }
0814   bool crossings_match = std::all_of(reco_daughters.begin(),reco_daughters.end(),[&](SvtxTrack* reco_track){ return reco_track->get_crossing() == reco_daughters[0]->get_crossing(); });
0815   if(Verbosity()>5)
0816   {
0817     std::cout << "crossings_match: " << crossings_match << std::endl;
0818   }
0819   return crossings_match;
0820 }
0821
0822 //____________________________________________________________________________..
0823 int ResonanceGeometricAcceptance::ResetEvent([[maybe_unused]] PHCompositeNode *topNode)
0824 {
0825   //std::cout << "ResonanceGeometricAcceptance::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0826   return Fun4AllReturnCodes::EVENT_OK;
0827 }
0828
0829 //____________________________________________________________________________..
0830 int ResonanceGeometricAcceptance::EndRun(const int runnumber)
0831 {
0832   std::cout << "ResonanceGeometricAcceptance::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0833   return Fun4AllReturnCodes::EVENT_OK;
0834 }
0835
0836 //____________________________________________________________________________..
0837 int ResonanceGeometricAcceptance::End([[maybe_unused]] PHCompositeNode *topNode)
0838 {
0839   //std::cout << "ResonanceGeometricAcceptance::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0840   if(Verbosity()>0)
0841   {
0842     std::cout << m_mother_name << " (passing / all) = " << passing_candidates_vspt->GetEntries() << " / " << all_candidates_vspt->GetEntries() << std::endl;
0843   }
0844
0845   outfile->cd();
0846   all_candidates_vspt->Write();
0847   all_candidates_vseta->Write();
0848   all_candidates_vsphi->Write();
0849   all_candidates_vsrapidity->Write();
0850   passing_candidates_vspt->Write();
0851   passing_candidates_vseta->Write();
0852   passing_candidates_vsphi->Write();
0853   passing_candidates_vsrapidity->Write();
0854   outfile->Close();
0855
0856   return Fun4AllReturnCodes::EVENT_OK;
0857 }
0858
0859 //____________________________________________________________________________..
0860 int ResonanceGeometricAcceptance::Reset([[maybe_unused]] PHCompositeNode *topNode)
0861 {
0862   //std::cout << "ResonanceGeometricAcceptance::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0863   return Fun4AllReturnCodes::EVENT_OK;
0864 }
0865
0866 //____________________________________________________________________________..
0867 void ResonanceGeometricAcceptance::Print(const std::string &what) const
0868 {
0869   std::cout << "ResonanceGeometricAcceptance::Print(const std::string &what) const Printing info for " << what << std::endl;
0870 }