DecayFinder.h

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DecayFinder.h

#ifndef DECAYFINDER_DECAYFINDER_H
#define DECAYFINDER_DECAYFINDER_H

// sPHENIX stuff
#include <fun4all/SubsysReco.h>

#include <g4main/PHG4Particle.h>
#include <g4main/PHG4VtxPoint.h>

#include <cstddef>  // for NULL
#include <string>
#include <utility>  // for pair
#include <vector>

class DecayFinderContainer_v1;
class PHCompositeNode;
class PHG4TruthInfoContainer;
class PHHepMCGenEvent;
class PHHepMCGenEventMap;
namespace HepMC
{
  class GenParticle;
}

class DecayFinder : public SubsysReco
{
 public:
  using Decay = std::vector<std::pair<std::pair<int, int>, int>>;

  DecayFinder();

  explicit DecayFinder(const std::string &name);

  ~DecayFinder() override = default;

  int Init(PHCompositeNode *topNode) override;

  int process_event(PHCompositeNode *topNode) override;

  int End(PHCompositeNode *topNode) override;

  int parseDecayDescriptor();

  bool findDecay(PHCompositeNode *topNode);

  bool findParticle(const std::string &particle);

  void searchHepMCRecord(HepMC::GenParticle *particle, std::vector<int> decayProducts,
                         bool &breakLoop, bool &hasPhoton, bool &hasPi0, bool &failedPT, bool &failedETA,
                         std::vector<int> &correctDecayProducts);

  void searchGeant4Record(int barcode, int pid, std::vector<int> decayProducts,
                          bool &breakLoop, bool &hasPhoton, bool &hasPi0, bool &failedPT, bool &failedETA,
                          std::vector<int> &correctDecayProducts);

  bool checkIfCorrectHepMCParticle(HepMC::GenParticle *particle, bool &trackFailedPT, bool &trackFailedETA);

  bool checkIfCorrectGeant4Particle(PHG4Particle *particle, bool &hasPhoton, bool &hasPi0, bool &trackFailedPT, bool &trackFailedETA);

  bool compareDecays(std::vector<int> required, std::vector<int> actual);

  int deleteElement(int arr[], int n, int x);

  void multiplyVectorByScalarAndSort(std::vector<int> &v, int k);

  int get_pdgcode(const std::string &name);

  int get_charge(const std::string &name);

  bool isInRange(float min, float value, float max);

  int createDecayNode(PHCompositeNode *topNode);

  void fillDecayNode(PHCompositeNode *topNode, Decay &decay);

  void printInfo();

  void printNode(PHCompositeNode *topNode);

  // User configuration
  void setDecayDescriptor(const std::string &decayDescriptor) { m_decayDescriptor = decayDescriptor; }
  void triggerOnDecay(bool trigger) { m_triggerOnDecay = trigger; }
  void allowPhotons(bool allow) { m_allowPhotons = allow; }
  void allowPi0(bool allow) { m_allowPi0 = allow; }
  void saveDST(bool save) { m_save_dst = save; }
  void setNodeName(const std::string &name) { m_container_name = name; }

  void setEtaRange(float min, float max)
  {
    m_eta_low_req = min;
    m_eta_high_req = max;
  }

  void setPTmin(float pt) { m_pt_req = pt; }

  void useDecaySpecificEtaRange(bool use) { m_recalcualteEtaRange = use; }

 private:
  PHHepMCGenEventMap *m_geneventmap = nullptr;
  PHHepMCGenEvent *m_genevt = nullptr;
  PHG4TruthInfoContainer *m_truthinfo = nullptr;

  void recalculateEta(double py, double vertex[3]);
  void calculateEffectiveTPCradius(double vertex[3], double &effective_top_r, double &effective_bottom_r);
  bool m_recalcualteEtaRange = true;
  double m_tpc_r = 78.0;
  double m_tpc_z = 105.0;
  double m_effective_top_tpc_r = m_tpc_r;
  double m_effective_bottom_tpc_r = m_tpc_r;

  double m_eta_high_req = 1.1;
  double m_eta_low_req = -1.1;
  double m_pt_req = 0.2;

  int m_counter = 0;
  int m_intermediate_product_counter = 0;
  int m_nCandFail_pT = 0;
  int m_nCandFail_eta = 0;
  int m_nCandFail_pT_and_eta = 0;
  int m_nCandReconstructable = 0;
  int m_nCandHas_Photon = 0;
  int m_nCandHas_Pi0 = 0;
  int m_nCandHas_Photon_and_Pi0 = 0;
  int m_nCandHas_noPhoton_and_noPi0 = 0;

  bool m_getChargeConjugate = false;

  std::string m_decayDescriptor;
  bool m_triggerOnDecay = false;
  bool m_allowPi0 = false;
  bool m_allowPhotons = false;

  int m_mother_ID = 0;
  std::vector<int> m_intermediates_ID;
  std::vector<int> m_daughters_ID;

  int m_nTracksFromMother = 0;
  std::vector<int> m_nTracksFromIntermediates;

  std::vector<int> m_motherDecayProducts;

  bool m_save_dst;
  DecayFinderContainer_v1 *m_decayMap = nullptr;
  Decay decayChain;
  std::string m_nodeName;
  std::string m_container_name;
};

#endif  // DECAYFINDER_DECAYFINDER_H