geometric_acceptance/ResonanceGeometricAcceptance/ResonanceGeometricAcceptance.h

0001 // Tell emacs that this is a C++ source
0002 //  -*- C++ -*-.
0003 #ifndef RESONANCE_GEOMETRICACCEPTANCE_H
0004 #define RESONANCE_GEOMETRICACCEPTANCE_H
0005
0006 #include <fun4all/SubsysReco.h>
0007
0008 #include "../../config/binning.h"
0009 #include "../../util/HistogramTools.h"
0010
0011 #include <trackbase_historic/SvtxTrackMap.h>
0012 #include <trackbase_historic/PHG4ParticleSvtxMap.h>
0013 #include <trackbase_historic/SvtxPHG4ParticleMap.h>
0014
0015 #include <globalvertex/SvtxVertexMap.h>
0016
0017 #include <g4main/PHG4TruthInfoContainer.h>
0018
0019 #include <KFParticle.h>
0020
0021 #include <TFile.h>
0022 #include <string>
0023
0024 class PHCompositeNode;
0025
0026 class ResonanceGeometricAcceptance : public SubsysReco
0027 {
0028  public:
0029
0030   ResonanceGeometricAcceptance(const std::string &name = "ResonanceGeometricAcceptance");
0031
0032   ~ResonanceGeometricAcceptance() override;
0033
0034   /** Called during initialization.
0035       Typically this is where you can book histograms, and e.g.
0036       register them to Fun4AllServer (so they can be output to file
0037       using Fun4AllServer::dumpHistos() method).
0038    */
0039   int Init(PHCompositeNode *topNode) override;
0040
0041   /** Called for first event when run number is known.
0042       Typically this is where you may want to fetch data from
0043       database, because you know the run number. A place
0044       to book histograms which have to know the run number.
0045    */
0046   int InitRun(PHCompositeNode *topNode) override;
0047
0048   /** Called for each event.
0049       This is where you do the real work.
0050    */
0051   int process_event(PHCompositeNode *topNode) override;
0052
0053   /// Clean up internals after each event.
0054   int ResetEvent(PHCompositeNode *topNode) override;
0055
0056   /// Called at the end of each run.
0057   int EndRun(const int runnumber) override;
0058
0059   /// Called at the end of all processing.
0060   int End(PHCompositeNode *topNode) override;
0061
0062   /// Reset
0063   int Reset(PHCompositeNode * /*topNode*/) override;
0064
0065   void Print(const std::string &what = "ALL") const override;
0066
0067   void setTrackMapName(const std::string& name)
0068   {
0069     m_trackmap_name = name;
0070   }
0071
0072   void setVertexMapName(const std::string& name)
0073   {
0074     m_vertexmap_name = name;
0075   }
0076
0077   void setTruthContainerName(const std::string& name)
0078   {
0079     m_truth_name = name;
0080   }
0081
0082   void setOutputFilename(const std::string& name)
0083   {
0084     m_outfile_name = name;
0085   }
0086
0087   void setMotherName(const std::string& name)
0088   {
0089     m_mother_name = name;
0090   }
0091
0092   void setMotherPDGID(const int pdgid)
0093   {
0094     m_mother_pdgid = pdgid;
0095   }
0096
0097   void setDaughterPDGIDs(const std::vector<int>& daughter_pdgids)
0098   {
0099     m_daughter_pdgids = daughter_pdgids;
0100   }
0101
0102   void includeConjugate(const bool include = true)
0103   {
0104     m_include_conjugate = include;
0105   }
0106
0107   void setTruthWeightThreshold(const int threshold)
0108   {
0109     m_truth_weight_threshold = threshold;
0110   }
0111
0112   void setDaughterMinPt(const float minpt)
0113   {
0114     m_daughter_min_pT = minpt;
0115   }
0116
0117   void SetVertexZLimits(const std::pair<float,float> limits)
0118   {
0119     zvertex_limits = limits;
0120   }
0121
0122   void setPtBins(const HistogramInfo& hinfo)
0123   {
0124     pt_bins = hinfo;
0125   }
0126
0127   void setEtaBins(const HistogramInfo& hinfo)
0128   {
0129     eta_bins = hinfo;
0130   }
0131
0132   void setPhiBins(const HistogramInfo& hinfo)
0133   {
0134     phi_bins = hinfo;
0135   }
0136
0137   void setRapidityBins(const HistogramInfo& hinfo)
0138   {
0139     rapidity_bins = hinfo;
0140   }
0141
0142   void setPtLimits(float vmin, float vmax)
0143   {
0144     pt_limits = {vmin, vmax};
0145   }
0146
0147   void setEtaLimits(float vmin, float vmax)
0148   {
0149     eta_limits = {vmin, vmax};
0150   }
0151
0152   void setPhiLimits(float vmin, float vmax)
0153   {
0154     phi_limits = {vmin, vmax};
0155   }
0156
0157   void setRapidityLimits(float vmin, float vmax)
0158   {
0159     rapidity_limits = {vmin, vmax};
0160   }
0161
0162  private:
0163
0164   int get_best_reco_vertex(KFParticle mother, const int crossing) const;
0165   KFParticle make_KFParticle_Vertex(SvtxVertex* vertex) const;
0166   KFParticle make_KFParticle(PHG4Particle* particle, SvtxTrack* track) const;
0167   std::set<SvtxTrack*> get_best_reco_matches(const PHG4Particle* particle) const;
0168   std::set<SvtxTrack*> get_reco_matches_with_silicon_hits(const PHG4Particle* particle) const;
0169   std::vector<std::vector<SvtxTrack*>> get_all_reco_combinations(const std::vector<PHG4Particle*>& daughters) const;
0170   bool has_all_daughters(const std::vector<PHG4Particle*>& daughters) const;
0171   bool has_daughter(const int pdgid, const std::vector<PHG4Particle*>& daughters) const;
0172   bool check_charge_swaps(const std::vector<PHG4Particle*>& daughters, const std::vector<SvtxTrack*>& reco_daughters) const;
0173   bool check_crossings_match(const std::vector<SvtxTrack*>& reco_daughters) const;
0174   void identify(KFParticle kfp) const;
0175   int get_mother_PDGID(const std::vector<PHG4Particle*>& daughters) const;
0176   int get_nmaps(TrackSeed* si_seed) const;
0177   //bool truth_track_is_best_match_to_reco_track(const PHG4Particle* particle) const;
0178
0179   HistogramInfo pt_bins = BinInfo::final_pt_bins;
0180   HistogramInfo eta_bins = BinInfo::final_eta_bins;
0181   HistogramInfo phi_bins = BinInfo::final_phi_bins;
0182   HistogramInfo rapidity_bins = BinInfo::final_rapidity_bins;
0183
0184   std::pair<float,float> zvertex_limits = {-10.,10.}; // cm
0185
0186   std::pair<float,float> pt_limits = {BinInfo::final_pt_bins.bins.front(),BinInfo::final_pt_bins.bins.back()};
0187   std::pair<float,float> eta_limits = {BinInfo::final_eta_bins.bins.front(),BinInfo::final_eta_bins.bins.back()};
0188   std::pair<float,float> phi_limits = {BinInfo::final_phi_bins.bins.front(),BinInfo::final_phi_bins.bins.back()};
0189   std::pair<float,float> rapidity_limits = {BinInfo::final_rapidity_bins.bins.front(),BinInfo::final_rapidity_bins.bins.back()};
0190
0191   TH1F* all_candidates_vspt;
0192   TH1F* all_candidates_vseta;
0193   TH1F* all_candidates_vsphi;
0194   TH1F* all_candidates_vsrapidity;
0195
0196   TH1F* passing_candidates_vspt;
0197   TH1F* passing_candidates_vseta;
0198   TH1F* passing_candidates_vsphi;
0199   TH1F* passing_candidates_vsrapidity;
0200
0201   TFile* outfile;
0202
0203   int m_mother_pdgid = 0;
0204   std::vector<int> m_daughter_pdgids;
0205   bool m_include_conjugate = false;
0206   int m_truth_weight_threshold = 20;
0207
0208   float m_daughter_min_pT = 0.;
0209
0210   std::string m_mother_name = "Particle";
0211   std::string m_trackmap_name = "SvtxTrackMap";
0212   std::string m_vertexmap_name = "SvtxVertexMap";
0213   std::string m_truth_name = "G4TruthInfo";
0214   std::string m_outfile_name = "geometric_acceptance.root";
0215
0216   SvtxTrackMap* m_trackmap;
0217   SvtxVertexMap* m_vertexmap;
0218   PHG4TruthInfoContainer* m_truth;
0219   PHG4ParticleSvtxMap* m_truth_reco_map;
0220   SvtxPHG4ParticleMap* m_reco_truth_map;
0221 };
0222
0223 #endif // GEOMETRICACCEPTANCE_H