pythiaEMCalAna/src/pythiaEMCalAna.h

0001 // Tell emacs that this is a C++ source
0002 //  -*- C++ -*-.
0003 #ifndef ISOCLUSTER_H
0004 #define ISOCLUSTER_H
0005
0006 #include <fun4all/SubsysReco.h>
0007
0008 #include <string>
0009 #include <set>
0010 #include <vector>
0011
0012 #include <HepMC/GenParticle.h>
0013 #include <g4main/PHG4TruthInfoContainer.h>
0014 #include <g4main/PHG4Particle.h>
0015 #include <g4eval/CaloRawClusterEval.h>
0016
0017 class PHCompositeNode;
0018 class TFile;
0019 class TTree;
0020 namespace HepMC {
0021     class GenEvent;
0022 }
0023 class RawClusterContainer;
0024
0025 class pythiaEMCalAna : public SubsysReco
0026 {
0027  public:
0028
0029   pythiaEMCalAna(const std::string &name = "pythiaEMCalAna", const std::string &oname = "pythiaEMCalAnaTrees.root", bool isMC = true, bool hasPythia = true);
0030
0031   ~pythiaEMCalAna() override;
0032
0033   /** Called during initialization.
0034       Typically this is where you can book histograms, and e.g.
0035       register them to Fun4AllServer (so they can be output to file
0036       using Fun4AllServer::dumpHistos() method).
0037    */
0038   int Init(PHCompositeNode *topNode) override;
0039
0040   /** Called for first event when run number is known.
0041       Typically this is where you may want to fetch data from
0042       database, because you know the run number. A place
0043       to book histograms which have to know the run number.
0044    */
0045   int InitRun(PHCompositeNode *topNode) override;
0046
0047   /** Called for each event.
0048       This is where you do the real work.
0049    */
0050   int process_event(PHCompositeNode *topNode) override;
0051
0052   /// Clean up internals after each event.
0053   int ResetEvent(PHCompositeNode *topNode) override;
0054
0055   /// Called at the end of each run.
0056   int EndRun(const int runnumber) override;
0057
0058   /// Called at the end of all processing.
0059   int End(PHCompositeNode *topNode) override;
0060
0061   /// Reset
0062   int Reset(PHCompositeNode * /*topNode*/) override;
0063
0064   void Print(const std::string &what = "ALL") const override;
0065
0066   void setGenEvent(int eventGet)     {getEvent = eventGet;}
0067
0068  private:
0069
0070   TTree *clusters_Towers;
0071   TTree *truth_particles;
0072   PHG4TruthInfoContainer* m_truthInfo;
0073   HepMC::GenEvent* m_theEvent;
0074   RawClusterContainer* m_clusterContainer;
0075   CaloRawClusterEval* m_clusterEval;
0076
0077   //stuff for towers and clusters
0078   std::vector<float> m_tower_energy;
0079   std::vector<float> m_eta_center;
0080   std::vector<float> m_phi_center;
0081   std::vector<float> m_cluster_ID;
0082   std::vector<float> m_cluster_e;
0083   std::vector<float> m_cluster_eta;
0084   std::vector<float> m_cluster_phi;
0085   std::vector<float> m_cluster_ecore;
0086   std::vector<float> m_cluster_chi2;
0087   std::vector<float> m_cluster_prob;
0088   std::vector<float> m_cluster_nTowers;
0089   std::vector<std::vector<float>> m_cluster_allTowersE;
0090   std::vector<std::vector<float>> m_cluster_allTowersEta;
0091   std::vector<std::vector<float>> m_cluster_allTowersPhi;
0092   std::vector<float> m_cluster_nParticles;
0093   std::vector<float> m_cluster_primaryParticlePid;
0094   std::vector<std::vector<float>> m_cluster_allSecondaryPids;
0095   /* std::vector<float> m_cluster_maxE_trackID; */
0096   /* std::vector<float> m_cluster_maxE_PID; */
0097   /* std::vector<std::vector<float>> m_cluster_all_trackIDs; */
0098
0099   //truth particle information
0100   /* std::vector<float> m_truthIsPrimary; */
0101   /* std::vector<float> m_truthTrackID; */
0102   /* std::vector<float> m_truthPid; */
0103   /* std::vector<float> m_truthE; */
0104   /* std::vector<float> m_truthEta; */
0105   /* std::vector<float> m_truthPhi; */
0106   /* std::vector<float> m_truthPt; */
0107   /* std::vector<float> m_truthMass; */
0108   /* std::vector<float> m_truthEndVtx_x; */
0109   /* std::vector<float> m_truthEndVtx_y; */
0110   /* std::vector<float> m_truthEndVtx_z; */
0111   /* std::vector<float> m_truthEndVtx_t; */
0112   /* std::vector<float> m_truthEndVtx_r; */
0113   /* std::vector<std::vector<float>> m_truth_all_clusterIDs; */
0114
0115   // V3
0116   std::vector<float> m_truth_Parent_Barcode;
0117   std::vector<float> m_truth_Parent_Pid;
0118   std::vector<float> m_truth_Parent_M;
0119   std::vector<float> m_truth_Parent_E;
0120   std::vector<float> m_truth_Parent_Eta;
0121   std::vector<float> m_truth_Parent_Phi;
0122   std::vector<float> m_truth_Parent_Pt;
0123   std::vector<float> m_truth_Parent_xF;
0124   std::vector<float> m_truth_Parent_EndVtx_x;
0125   std::vector<float> m_truth_Parent_EndVtx_y;
0126   std::vector<float> m_truth_Parent_EndVtx_z;
0127   /* std::vector<float> m_truth_Parent_EndVtx_t; */
0128   std::vector<float> m_truth_Parent_EndVtx_r;
0129   /* std::vector<float> m_truth_Parent_TotalNDaughters; */
0130   /* std::vector<std::vector<float>> m_truth_Parent_AllDaughterPids; */
0131   /* std::vector<std::vector<float>> m_truth_Parent_AllDaughterEnergyFractions; */
0132   std::vector<float> m_truth_Parent_TotalNClusters;
0133   std::vector<std::vector<float>> m_truth_Parent_AllClusterIDs;
0134   std::vector<std::vector<float>> m_truth_Parent_AllClusterEnergyFractions;
0135   std::vector<float> m_truth_Decay1_Barcode;
0136   std::vector<float> m_truth_Decay1_Pid;
0137   std::vector<float> m_truth_Decay1_M;
0138   std::vector<float> m_truth_Decay1_E;
0139   std::vector<float> m_truth_Decay1_Eta;
0140   std::vector<float> m_truth_Decay1_Phi;
0141   std::vector<float> m_truth_Decay1_Pt;
0142   std::vector<float> m_truth_Decay1_xF;
0143   std::vector<float> m_truth_Decay1_EndVtx_x;
0144   std::vector<float> m_truth_Decay1_EndVtx_y;
0145   std::vector<float> m_truth_Decay1_EndVtx_z;
0146   /* std::vector<float> m_truth_Decay1_EndVtx_t; */
0147   std::vector<float> m_truth_Decay1_EndVtx_r;
0148   std::vector<float> m_truth_Decay1_TotalNClusters;
0149   std::vector<float> m_truth_Decay1_BestClusterID;
0150   std::vector<float> m_truth_Decay1_BestClusterEfraction;
0151   /* std::vector<std::vector<float>> m_truth_Decay1_AllClusterIDs; */
0152   /* std::vector<std::vector<float>> m_truth_Decay1_AllClusterEnergyFractions; */
0153   std::vector<float> m_truth_Decay2_Barcode;
0154   std::vector<float> m_truth_Decay2_Pid;
0155   std::vector<float> m_truth_Decay2_M;
0156   std::vector<float> m_truth_Decay2_E;
0157   std::vector<float> m_truth_Decay2_Eta;
0158   std::vector<float> m_truth_Decay2_Phi;
0159   std::vector<float> m_truth_Decay2_Pt;
0160   std::vector<float> m_truth_Decay2_xF;
0161   std::vector<float> m_truth_Decay2_EndVtx_x;
0162   std::vector<float> m_truth_Decay2_EndVtx_y;
0163   std::vector<float> m_truth_Decay2_EndVtx_z;
0164   /* std::vector<float> m_truth_Decay2_EndVtx_t; */
0165   std::vector<float> m_truth_Decay2_EndVtx_r;
0166   std::vector<float> m_truth_Decay2_TotalNClusters;
0167   std::vector<float> m_truth_Decay2_BestClusterID;
0168   std::vector<float> m_truth_Decay2_BestClusterEfraction;
0169   /* std::vector<std::vector<float>> m_truth_Decay2_AllClusterIDs; */
0170   /* std::vector<std::vector<float>> m_truth_Decay2_AllClusterEnergyFractions; */
0171   /* std::vector<float> m_truth_Cluster1_Id; */
0172   /* std::vector<float> m_truth_Cluster1_E; */
0173   /* std::vector<float> m_truth_Cluster1_Ecore; */
0174   /* std::vector<float> m_truth_Cluster1_Eta; */
0175   /* std::vector<float> m_truth_Cluster1_Phi; */
0176   /* std::vector<float> m_truth_Cluster1_Pt; */
0177   /* std::vector<float> m_truth_Cluster1_xF; */
0178   /* std::vector<float> m_truth_Cluster1_Chi2; */
0179   /* std::vector<float> m_truth_Cluster1_Decay1EnergyFraction; */
0180   /* std::vector<float> m_truth_Cluster1_Decay2EnergyFraction; */
0181   /* std::vector<float> m_truth_Cluster2_Id; */
0182   /* std::vector<float> m_truth_Cluster2_E; */
0183   /* std::vector<float> m_truth_Cluster2_Ecore; */
0184   /* std::vector<float> m_truth_Cluster2_Eta; */
0185   /* std::vector<float> m_truth_Cluster2_Phi; */
0186   /* std::vector<float> m_truth_Cluster2_Pt; */
0187   /* std::vector<float> m_truth_Cluster2_xF; */
0188   /* std::vector<float> m_truth_Cluster2_Chi2; */
0189   /* std::vector<float> m_truth_Cluster2_Decay1EnergyFraction; */
0190   /* std::vector<float> m_truth_Cluster2_Decay2EnergyFraction; */
0191
0192   TFile *fout;
0193   std::string outname;
0194   int getEvent;
0195   /* bool hasHIJING; // needed to handle HIJING embedded samples correctly */
0196   bool isMonteCarlo; // if input is RD we obviously need to skip the truth info
0197   bool hasPYTHIA; // used to determine if primary photons are truly direct photons
0198   // counters for events
0199   long int n_events_total;
0200   long int n_events_minbias;
0201   long int n_events_with_vertex;
0202   long int n_events_with_good_vertex;
0203   long int n_events_positiveCaloE;
0204   // counters for each primary particle type
0205   long int n_primaries;
0206   long int n_primary_photons;
0207   long int n_direct_photons;
0208   long int n_leptons;
0209   long int n_neutral_hadrons;
0210   long int n_neutral_hadrons_geant;
0211   long int n_neutral_hadrons_pythia;
0212   long int n_charged_hadrons;
0213   long int n_charged_hadrons_geant;
0214   long int n_charged_hadrons_pythia;
0215   long int n_pythia_decayed_hadrons;
0216   std::set<int> allTrackIDs; // keep track of all truth particles that produced a shower
0217   std::set<int> primaryBarcodes; // keep track of all PYTHIA-decayed primaries
0218
0219   /* std::vector<std::pair<int,int>> primaryBarcodes; */
0220   /* std::vector<std::pair<int,int>> secondaryBarcodes; */
0221   HepMC::GenParticle* getGenParticle(int barcode);
0222   PHG4Particle* getG4Particle(int barcode);
0223   bool withinAcceptance(PHG4Particle* part);
0224   bool withinAcceptance(HepMC::GenParticle* part);
0225   PHG4VtxPoint* getG4EndVtx(int id);
0226   bool isDirectPhoton(PHG4Particle* part);
0227   /* void addPrimaryFromGeant(PHG4Particle* part); */
0228   /* void addPrimaryFromPythia(HepMC::GenParticle* part); */
0229   void PrintParent(PHG4Particle* par);
0230   void FillTruthParticle(std::string which, PHG4Particle* par);
0231   void FillAllClustersFromParent(PHG4Particle* par);
0232   void FillParent(PHG4Particle* par);
0233   HepMC::GenParticle* GetHepMCParent(PHG4Particle* par);
0234   void FillTruth(HepMC::GenParticle* par);
0235   std::vector<PHG4Particle*> GetAllDaughters(PHG4Particle* parent);
0236   void GetBestDaughters(PHG4Particle* parent, PHG4Particle* &decay1, PHG4Particle* &decay2);
0237   void FillDecay(std::string which, PHG4Particle* decay, PHG4Particle* parent);
0238   RawCluster* FindBestCluster(PHG4Particle* decay, PHG4Particle* parent);
0239 };
0240
0241 #endif // ISOCLUSTER_H