HF-Particle/HFTrigger/HFTriggerMVA.h

0001 #ifndef HFTRIGGERMVA_H
0002 #define HFTRIGGERMVA_H
0003
0004 //sPHENIX stuff
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006 #include <fun4all/SubsysReco.h>
0007 #include <phool/getClass.h>
0008 #include <trackbase/TrkrDefs.h>
0009 #include <trackbase/TrkrHitSet.h>
0010 #include <trackbase/TrkrHitSetContainer.h>
0011 #include <trackbase_historic/SvtxTrack.h>
0012 #include <trackbase_historic/SvtxTrackMap.h>
0013 #include <trackbase_historic/SvtxVertex.h>
0014 #include <trackbase_historic/SvtxVertexMap.h>
0015
0016 #include <g4eval/SvtxEvalStack.h>
0017 #include <g4eval/SvtxTrackEval.h>
0018
0019 //ROOT stuff
0020 #include <TMVA/Reader.h>
0021 #include <TMVA/Tools.h>
0022
0023 #include <Eigen/Core>
0024 #include <Eigen/Dense>
0025
0026 #include <algorithm>
0027 #include <cfloat>
0028 #include <cmath>
0029 #include <iostream>
0030 #include <map>
0031 #include <numeric>
0032 #include <string>
0033
0034 typedef Eigen::Matrix<float, 6, 1> Track;
0035 typedef Eigen::Matrix<float, 3, 1> Vertex;
0036 typedef Eigen::Matrix<float, 3, 1> TrackX;
0037 typedef Eigen::Matrix<float, 3, 1> TrackP;
0038 typedef Eigen::Matrix<float, 3, 1> DCA;
0039
0040 class SvtxTrack;
0041 class SvtxTrackMap;
0042 class SvtxVertexMap;
0043 class SvtxVertex;
0044
0045 class HFTriggerMVA : public SubsysReco
0046 {
0047  public:
0048   HFTriggerMVA();
0049
0050   explicit HFTriggerMVA(const std::string &name);
0051
0052   virtual ~HFTriggerMVA(){}
0053
0054   int Init(PHCompositeNode *topNode);
0055
0056   int process_event(PHCompositeNode *topNode);
0057
0058   int End(PHCompositeNode *topNode);
0059
0060   bool runTrigger(PHCompositeNode *topNode);
0061
0062   bool runCutsTrigger(float mean_hits, float emcal_energy, float IP, float IP_xy, float DCA);
0063
0064   bool runMVATrigger(TMVA::Reader* reader, std::string method, std::vector<float> inputValues, float cut);
0065
0066   void calculateMultiplicity(PHCompositeNode *topNode, float& meanMultiplicity, float& asymmetryMultiplicity);
0067
0068   float getMaxEMCalEnergy(PHCompositeNode *topNode);
0069
0070   void getIPVariables(Track track1, Track track2, Vertex vertex, float& minIP, float& minIP_xy, float& maxIP, float& maxIP_xy, float& DCA);
0071
0072   Vertex makeVertex(PHCompositeNode *topNode);
0073
0074   std::vector<Vertex> makeAllPrimaryVertices(PHCompositeNode *topNode);
0075
0076   Track makeTrack(PHCompositeNode *topNode);
0077
0078   std::vector<Track> makeAllTracks(PHCompositeNode *topNode);
0079
0080   int decomposeTrack(Track track, TrackX& trackPosition, TrackP& trackMomentum);
0081
0082   float calcualteTrackVertex2DDCA(Track track, Vertex vertex);
0083
0084   float calcualteTrackVertexDCA(Track track, Vertex vertex);
0085
0086   float calcualteTrackTrackDCA(Track trackOne, Track trackTwo);
0087
0088   std::tuple<TMVA::Reader*, std::vector<Float_t>> initMVA(std::vector<std::string> variableList, std::string method, std::string mvaFile);
0089
0090   void printRunRecord();
0091
0092   void printTrigger();
0093
0094   //User configuration
0095   void requireCutsTrigger(bool useTrigger) { m_useCutsTrigger = useTrigger; }
0096   void requireMVAwCaloTrigger(bool useTrigger) { m_useMVAwCaloTrigger = useTrigger; }
0097
0098   void setMeanMult(float value) { meanMult = value; }
0099   void setAsymmMult(float value) { asymmMult = value; }
0100   void setTrackVertex3DDCA(float value) { trackVertex3DDCA = value; }
0101   void setTrackVertex2DDCA(float value) { trackVertex2DDCA = value; }
0102   void setTrackTrackDCA(float value) { trackTrackDCA = value; }
0103   void setminEMCalEnergy(float value) { minEMCalEnergy = value; }
0104   void setMVA_wCaloResponse(float value) { MVA_wCaloResponse = value; }
0105   void setMVA_woutCaloResponse(float value) { MVA_woutCaloResponse = value; }
0106   void setMVA_woutCaloOrMinTrackResponse(float value) { MVA_woutCaloOrMinTrackResponse = value; }
0107
0108  protected:
0109   SvtxEvalStack *m_svtx_evalstack = nullptr;
0110   SvtxTrackEval *trackeval = nullptr;
0111
0112  private:
0113
0114   float meanMult = 12;
0115   float asymmMult = 0.1;
0116   float trackVertex3DDCA = 0.05;
0117   float trackVertex2DDCA = 0.05;
0118   float trackTrackDCA = 0.05;
0119   float minEMCalEnergy = 0.5;
0120   float MVA_wCaloResponse = -0.05;
0121   float MVA_woutCaloResponse = -0.05;
0122   float MVA_woutCaloOrMinTrackResponse = -0.05;
0123
0124   bool m_useCutsTrigger = false;
0125   bool m_useCutswoutTrigger = false;
0126   bool m_useMVAwCaloTrigger = false;
0127   bool m_useMVAwoutCaloTrigger = false;
0128   bool m_useMVAwoutCaloAndMinTrackTrigger = false;
0129
0130   SvtxVertexMap *m_dst_vertexmap = nullptr;
0131   SvtxTrackMap *m_dst_trackmap = nullptr;
0132   SvtxVertex *m_dst_vertex = nullptr;
0133   SvtxTrack *m_dst_track = nullptr;
0134
0135   //Counters
0136   int m_events = 0;
0137   int m_no_trigger = 0;
0138   int m_counter = 0;
0139   int m_cuts_and_mva_wcalo_counter = 0;
0140   int m_no_cuts_and_mva_wcalo_counter = 0;
0141   int m_cuts_and_no_mva_wcalo_counter = 0;
0142   int m_no_cuts_and_no_mva_wcalo_counter = 0;
0143   int m_cuts_and_mva_woutcalo_counter = 0;
0144   int m_no_cuts_and_mva_woutcalo_counter = 0;
0145   int m_cuts_and_no_mva_woutcalo_counter = 0;
0146   int m_no_cuts_and_no_mva_woutcalo_counter = 0;
0147   int m_cuts_and_mva_woutcalo_and_mintracks_counter = 0;
0148   int m_no_cuts_and_mva_woutcalo_and_mintracks_counter = 0;
0149   int m_cuts_and_no_mva_woutcalo_and_mintracks_counter = 0;
0150   int m_no_cuts_and_no_mva_woutcalo_and_mintracks_counter = 0;
0151
0152   //MVA stuff
0153   std::string path = "/sphenix/u/cdean/software/analysis/HF-Particle/HFTrigger/weights/";
0154
0155   std::string mvaType = "BDTG";
0156
0157   std::vector<std::string> varListwCalo = {"min(track_1_IP,track_2_IP)",
0158                                            "min(abs(track_1_IP_xy),abs(track_2_IP_xy))",
0159                                            "max(track_1_IP,track_2_IP)",
0160                                            "max(abs(track_1_IP_xy),abs(track_2_IP_xy))",
0161                                            "max(track_1_EMCAL_energy_cluster,track_2_EMCAL_energy_cluster)",
0162                                            "track_1_track_2_DCA",
0163                                            "INTT_meanHits"};
0164
0165   std::vector<std::string> varListwoutCalo = {"min(track_1_IP,track_2_IP)",
0166                                               "min(abs(track_1_IP_xy),abs(track_2_IP_xy))",
0167                                               "max(track_1_IP,track_2_IP)",
0168                                               "max(abs(track_1_IP_xy),abs(track_2_IP_xy))",
0169                                               "track_1_track_2_DCA",
0170                                               "INTT_meanHits"};
0171
0172   std::vector<std::string> varListwoutCaloAndMinTrack = {"max(track_1_IP,track_2_IP)",
0173                                                          "max(abs(track_1_IP_xy),abs(track_2_IP_xy))",
0174                                                          "track_1_track_2_DCA",
0175                                                          "INTT_meanHits"};
0176
0177   TMVA::Reader *wCaloReader, *woutCaloReader, *woutCaloAndMinTrackReader;
0178   std::vector<float> wCaloFloats, woutCaloFloats, woutCaloAndMinTrackFloats;
0179 };
0180
0181 #endif  //HFTRIGGERMVA_H