packages/decayfinder/DecayFinder.h

0001 #ifndef DECAYFINDER_DECAYFINDER_H
0002 #define DECAYFINDER_DECAYFINDER_H
0003
0004 // sPHENIX stuff
0005 #include <fun4all/SubsysReco.h>
0006
0007 #include <g4main/PHG4Particle.h>
0008 #include <g4main/PHG4VtxPoint.h>
0009
0010 #include <cstddef>  // for NULL
0011 #include <string>
0012 #include <utility>  // for pair
0013 #include <vector>
0014
0015 class DecayFinderContainer_v1;
0016 class PHCompositeNode;
0017 class PHG4TruthInfoContainer;
0018 class PHHepMCGenEvent;
0019 class PHHepMCGenEventMap;
0020 namespace HepMC
0021 {
0022   class GenParticle;
0023 }
0024
0025 class DecayFinder : public SubsysReco
0026 {
0027  public:
0028   using Decay = std::vector<std::pair<std::pair<int, int>, int>>;
0029
0030   DecayFinder();
0031
0032   explicit DecayFinder(const std::string &name);
0033
0034   ~DecayFinder() override = default;
0035
0036   int Init(PHCompositeNode *topNode) override;
0037
0038   int process_event(PHCompositeNode *topNode) override;
0039
0040   int End(PHCompositeNode *topNode) override;
0041
0042   int parseDecayDescriptor();
0043
0044   bool findDecay(PHCompositeNode *topNode);
0045
0046   bool findParticle(const std::string &particle);
0047
0048   void searchHepMCRecord(HepMC::GenParticle *particle, std::vector<int> decayProducts,
0049                          bool &breakLoop, bool &hasPhoton, bool &hasPi0, bool &failedPT, bool &failedETA,
0050                          std::vector<int> &correctDecayProducts);
0051
0052   void searchGeant4Record(int barcode, int pid, std::vector<int> decayProducts,
0053                           bool &breakLoop, bool &hasPhoton, bool &hasPi0, bool &failedPT, bool &failedETA,
0054                           std::vector<int> &correctDecayProducts);
0055
0056   bool checkIfCorrectHepMCParticle(HepMC::GenParticle *particle, bool &trackFailedPT, bool &trackFailedETA);
0057
0058   bool checkIfCorrectGeant4Particle(PHG4Particle *particle, bool &hasPhoton, bool &hasPi0, bool &trackFailedPT, bool &trackFailedETA);
0059
0060   bool compareDecays(std::vector<int> required, std::vector<int> actual);
0061
0062   int deleteElement(int arr[], int n, int x);
0063
0064   void multiplyVectorByScalarAndSort(std::vector<int> &v, int k);
0065
0066   int get_pdgcode(const std::string &name);
0067
0068   int get_charge(const std::string &name);
0069
0070   bool isInRange(float min, float value, float max);
0071
0072   int createDecayNode(PHCompositeNode *topNode);
0073
0074   void fillDecayNode(PHCompositeNode *topNode, Decay &decay);
0075
0076   void printInfo();
0077
0078   void printNode(PHCompositeNode *topNode);
0079
0080   // User configuration
0081   /**
0082    * Use this function to define the decay you want to find in the HepMC record
0083    * @param[in] decayDescriptor the description of the decay chain, this is a string
0084    * You define the decay with these rules:
0085    * @brief You define a particle decaying with "->", the mother on the left, the decay products on the right
0086    * @brief Set the charge of final state tracks with "^", the particle name on the left and the charge on the right.
0087    *        Accepted charges are +, - and 0
0088    * @brief Use the same rules as above for any intermediatee decays but contain the entire decay within curled
0089    *        brackets, "{}"
0090    * @brief If you also want to find the charge conjugate decay, contain the entire decay descriptor within "[]cc" for
0091    *        charge-conjugate. The "cc" is NOT case sensitive
0092    * @brief The particle names you use must be kept in the TDatabasePDG class from root
0093    *        (https://root.cern.ch/doc/master/classTDatabasePDG.html). Print this table to see available particles with
0094    *        TDatabasePDG::Instance()->Print()
0095    * @brief An example of a decay would be: "[B+ -> {D0_bar -> kaon^+ pion^-} pion^+]cc"
0096    * @note There is an internal list of resonances which, if they appear in the record, will be further analysed. For
0097    *       example, the f0(980)->pipi decay is too quick to have a flight distance and so we would only see the pion
0098    *       pair in the detector. If you are looking for B_s0 -> J/psi pipi then the decay of the f0 will be studied
0099    *       for a pipi final state, basically inclusive decays are handled automatically. If you wish to study the f0
0100    *       decay, add it to your decay descriptor and it will automatically be removed from the "skip list"
0101    */
0102   void setDecayDescriptor(const std::string &decayDescriptor) { m_decayDescriptor = decayDescriptor; }
0103   /**
0104    * @param[in] trigger Set to true to allow further processing of events in which your decay appears, if your decay
0105    *            does not appear, all further processing of this event is skipped. This defaults to false so every event
0106    *            is proccessed in F4A
0107    */
0108   void triggerOnDecay(bool trigger) { m_triggerOnDecay = trigger; }
0109   /**
0110    * @param[in] allow Set to true to allow photons to be associated to your decay
0111    */
0112   void allowPhotons(bool allow) { m_allowPhotons = allow; }
0113   /**
0114    * @param[in] allow Set to true to allow pi zero to be associated to your decay
0115    */
0116   void allowPi0(bool allow) { m_allowPi0 = allow; }
0117   /**
0118    * @param[in] allow Set to true to save any of your decays that are found back to the node tree in a DecayFinderContainer
0119    *           The default name is "decay" and will automatically have "_DecayMap" added to the end
0120    */
0121   void saveDST(bool save) { m_save_dst = save; }
0122   /**
0123    * @param[in] name Change the default name of the DecayFinderContainer.
0124    * @note This name will still have "_DecayMap" added to the end, this cannot be changed
0125    */
0126   void setNodeName(const std::string &name) { m_container_name = name; }
0127
0128   /**
0129    * @param[in] min The minimum eta threshold for track acceptance
0130    * @param[in] min The maximum eta threshold for track acceptance
0131    * @note Set a pseudorapidity threshold range for tracking
0132    */
0133   void setEtaRange(float min, float max)
0134   {
0135     m_eta_low_req = min;
0136     m_eta_high_req = max;
0137   }
0138
0139   /**
0140    * @param[in] pt The minimum pT threshold for track acceptance
0141    * @note Set a minimum pT threshold for tracking
0142    */
0143   void setPTmin(float pt) { m_pt_req = pt; }
0144
0145   /**
0146    *  Set this to true to recalcualte a tracks acceptible pseudorapidity range in sPHENIX
0147    *  on a decay-by-decay basis using the secondary vertex position. Set to false to use the default
0148    *  or user specified range
0149    *  @param[in] True to recalculate the eta requirements per decay, false to use a fixed value (default = true)
0150    */
0151   void useDecaySpecificEtaRange(bool use) { m_recalcualteEtaRange = use; }
0152
0153   std::string passOrFail(bool condition);
0154
0155  private:
0156   PHHepMCGenEventMap *m_geneventmap = nullptr;
0157   PHHepMCGenEvent *m_genevt = nullptr;
0158   PHG4TruthInfoContainer *m_truthinfo = nullptr;
0159
0160   void recalculateEta(double py, double vertex[3]);
0161   void calculateEffectiveTPCradius(double vertex[3], double &effective_top_r, double &effective_bottom_r);
0162   bool m_recalcualteEtaRange = true;
0163   double m_tpc_r = 78.0;
0164   double m_tpc_z = 105.0;
0165   double m_effective_top_tpc_r = m_tpc_r;
0166   double m_effective_bottom_tpc_r = m_tpc_r;
0167
0168   double m_eta_high_req = 1.1;
0169   double m_eta_low_req = -1.1;
0170   double m_pt_req = 0.2;
0171
0172   int m_counter = 0;
0173   int m_intermediate_product_counter = 0;
0174   int m_nCandFail_pT = 0;
0175   int m_nCandFail_eta = 0;
0176   int m_nCandFail_pT_and_eta = 0;
0177   int m_nCandReconstructable = 0;
0178   int m_nCandHas_Photon = 0;
0179   int m_nCandHas_Pi0 = 0;
0180   int m_nCandHas_Photon_and_Pi0 = 0;
0181   int m_nCandHas_noPhoton_and_noPi0 = 0;
0182
0183   bool m_getChargeConjugate = false;
0184
0185   std::string m_decayDescriptor;
0186   bool m_triggerOnDecay = false;
0187   bool m_allowPi0 = false;
0188   bool m_allowPhotons = false;
0189
0190   int m_mother_ID = 0;
0191   std::vector<int> m_intermediates_ID;
0192   std::vector<int> m_daughters_ID;
0193
0194   int m_nTracksFromMother = 0;
0195   std::vector<int> m_nTracksFromIntermediates;
0196
0197   std::vector<int> m_motherDecayProducts;
0198
0199   bool m_save_dst;
0200   DecayFinderContainer_v1 *m_decayMap = nullptr;
0201   Decay decayChain;
0202   std::string m_nodeName;
0203   std::string m_container_name;
0204 };
0205
0206 #endif  // DECAYFINDER_DECAYFINDER_H