Adding the photon kinematics for the EMCal clusters QA type plots

offline/QA/Jet/Makefile.am

@@ -31,7 +31,8 @@ pkginclude_HEADERS = \
   TrksInJetQAJetManager.h \
   RhosinEvent.h \
   JetSeedCount.h \
-  JetKinematicCheck.h
+  JetKinematicCheck.h \
+  photonjetskinematics.h
 
 libjetqa_la_SOURCES = \
   ConstituentsinJets.cc\
@@ -47,7 +48,8 @@ libjetqa_la_SOURCES = \
   TrksInJetQAJetManager.cc \
   RhosinEvent.cc \
   JetSeedCount.cc \
-  JetKinematicCheck.cc
+  JetKinematicCheck.cc \
+  photonjetskinematics.cc
 
 libjetqa_la_LIBADD = \
   -L$(libdir) \

offline/QA/Jet/photonjetskinematics.cc

@@ -0,0 +1,220 @@
+//____________________________________________________________________________..
+
+#include "photonjetskinematics.h"
+#include <fun4all/PHTFileServer.h>
+
+// Fun4all includes
+
+#include <fun4all/Fun4AllHistoManager.h>
+
+#include <fun4all/Fun4AllReturnCodes.h>
+
+#include <phool/PHCompositeNode.h>
+
+#include <TFile.h>
+#include <phool/getClass.h>
+#include "TH1.h"
+#include "TMath.h"
+//#include <jetbase/JetContainer.h>
+
+
+// Tower includes
+#include <calobase/RawCluster.h>
+#include <calobase/RawClusterContainer.h>
+#include <calobase/RawClusterUtility.h>
+#include <calobase/RawTower.h>
+#include <calobase/RawTowerContainer.h>
+#include <calobase/RawTowerGeom.h>
+#include <calobase/RawTowerGeomContainer.h>
+#include <calobase/TowerInfo.h>
+#include <calobase/TowerInfoContainer.h>
+#include <calobase/TowerInfoDefs.h>
+
+// Cluster includes
+#include <calobase/RawCluster.h>
+#include <calobase/RawClusterContainer.h>
+
+// QA includes
+#include <qautils/QAHistManagerDef.h>
+
+// c++ includes
+#include <algorithm>
+#include <cassert>
+#include <iostream>
+#include <vector>
+//____________________________________________________________________________..
+
+photonjetskinematics::photonjetskinematics(const std::string &m_modulename, const std::string &m_inputnode) :
+  SubsysReco(m_modulename)
+  , modulename(m_modulename)
+  , inputnode(m_inputnode)
+  , histtag("AllTrig")
+  , trgToSelect(JetQADefs::GL1::MBDNSJet1)
+  , doTrgSelect(false)
+{
+  if (Verbosity() > 1) std::cout << "photonjetskinematics::photonjetskinematics(const std::string &name, const std::string&outputfilename) Calling ctor" << std::endl;
+}
+//____________________________________________________________________________..
+photonjetskinematics::~photonjetskinematics()
+{
+  std::cout << "photonjetskinematics::~photonjetskinematics() Calling dtor" << std::endl;
+}
+
+//____________________________________________________________________________..
+int photonjetskinematics::Init(PHCompositeNode* /*topNode*/)
+{
+  if (Verbosity() > 1) std::cout << "photonjetskinematics::Init(PHCompositeNode *topNode) Initializing" << std::endl;
+  manager = QAHistManagerDef::getHistoManager();
+  if (!manager)
+    {
+      std::cerr << PHWHERE << "PANIC: couldn't grab histogram manager!" << std::endl;
+      assert(manager);
+    }
+
+  // make sure module name is lower case
+  std::string smallModuleName = modulename;
+  std::transform(
+		 smallModuleName.begin(),
+		 smallModuleName.end(),
+		 smallModuleName.begin(),
+		 ::tolower);
+  // construct histogram names
+  std::vector<std::string> vecHistNames = {
+    "emcal_cluster_chi2",
+    "emcal_cluster_energy",
+    "emcal_cluster_eta_phi",
+    "emcal_cluster_eta",
+    "emcal_cluster_phi"};
+
+  for (auto& histName : vecHistNames)
+    {
+      histName.insert(0, "h_" + smallModuleName + "_");
+      if (!histtag.empty())
+	{
+	  histName.append("_" + histtag);
+	}
+    }
+
+  //initializing histograms
+
+  h_emcal_cluster_chi2 = new TH1D(vecHistNames[0].data(), "h_emcal_cluster_chi2", 30, 0, 150);
+  h_emcal_cluster_chi2->GetXaxis()->SetTitle("#chi^{2}");
+
+  h_emcal_cluster_energy = new TH1D(vecHistNames[1].data(), "h_emcal_cluster_energy", 100, 0, 10);
+  h_emcal_cluster_energy->GetXaxis()->SetTitle("E_{T} [GeV]");
+
+  h_emcal_cluster_eta_phi = new TH2F(vecHistNames[2].data(), "h_emcal_cluster_eta_phi", 140, -1.2, 1.2, 64, -M_PI, M_PI);
+  h_emcal_cluster_eta_phi->GetXaxis()->SetTitle("#eta");
+  h_emcal_cluster_eta_phi->GetYaxis()->SetTitle("#Phi");
+  h_emcal_cluster_eta_phi->SetOption("COLZ");
+
+  h_emcal_cluster_eta = new TH1D(vecHistNames[3].data(), "Eta Distribution", 140, -1.2, 1.2);
+  h_emcal_cluster_eta->GetXaxis()->SetTitle("#eta");
+
+  h_emcal_cluster_phi = new TH1D(vecHistNames[4].data(), "Phi Distribution", 64, -M_PI, M_PI);
+  h_emcal_cluster_phi->GetXaxis()->SetTitle("#phi");
+
+  return Fun4AllReturnCodes::EVENT_OK;
+
+}
+
+//____________________________________________________________________________..
+int photonjetskinematics::InitRun(PHCompositeNode* /*topNode*/)
+{
+
+  if (Verbosity() > 1) std::cout << "photonjetskinematics::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
+  return Fun4AllReturnCodes::EVENT_OK;
+}
+
+//____________________________________________________________________________..
+int photonjetskinematics::process_event(PHCompositeNode *topNode)
+{  
+
+  RawClusterContainer* clusterContainer = findNode::getClass<RawClusterContainer>(topNode, inputnode);
+  if (!clusterContainer)
+    {
+      std::cout << PHWHERE << "funkyCaloStuff::process_event - Fatal Error - CLUSTER_CEMC node is missing. " << std::endl;
+      return 0;
+    }
+
+  // if needed, check if trigger fired
+  if (doTrgSelect)
+    {
+      bool hasTrigger = JetQADefs::DidTriggerFire(trgToSelect, topNode);
+      if (!hasTrigger)
+	{
+	  return Fun4AllReturnCodes::EVENT_OK;
+	}
+    }
+
+  RawClusterContainer::ConstIterator clusterIter;
+  RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
+
+  for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
+    {
+      RawCluster* recoCluster = clusterIter->second;
+      CLHEP::Hep3Vector vertex(0, 0, 0);
+      CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
+
+      float clusE = E_vec_cluster.mag();
+      float clus_eta = E_vec_cluster.pseudoRapidity();
+      float clus_phi = E_vec_cluster.phi();
+      float clus_chisq = recoCluster->get_chi2();
+
+      //Filling Histograms, only taking into account E vec cluster, not reco cluster
+
+       h_emcal_cluster_chi2->Fill(clus_chisq);
+       h_emcal_cluster_energy->Fill(clusE);
+       h_emcal_cluster_eta_phi->Fill(clus_eta, clus_phi);
+       h_emcal_cluster_eta->Fill(clus_eta);  // 1D eta plot
+       h_emcal_cluster_phi->Fill(clus_phi);  // 1D phi plot
+   }
+
+  // std::cout << "photonjetskinematics::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
+   return Fun4AllReturnCodes::EVENT_OK;
+}
+
+//____________________________________________________________________________..
+int photonjetskinematics::ResetEvent(PHCompositeNode* /*topNode*/)
+{
+  // std::cout << "photonjetskinematics::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
+  return Fun4AllReturnCodes::EVENT_OK;
+}
+
+//____________________________________________________________________________..
+int photonjetskinematics::EndRun(const int runnumber)
+{
+  if (Verbosity() > 1) std::cout << "photonjetskinematics::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
+  return Fun4AllReturnCodes::EVENT_OK;
+}
+
+//____________________________________________________________________________..
+int photonjetskinematics::End(PHCompositeNode* /*topNode*/)
+{
+  // Commenting out the following line because Jenkins keeps failing the build test :( 
+  // if (Verbosity() > 1) std::cout << "photonjetskinematics::End - Output to " << outfilename << std::endl;
+
+  //Outputting the histograms
+  manager->registerHisto(h_emcal_cluster_eta_phi);
+  manager->registerHisto(h_emcal_cluster_energy);
+  manager->registerHisto(h_emcal_cluster_chi2);
+  manager->registerHisto(h_emcal_cluster_eta);  
+  manager->registerHisto(h_emcal_cluster_phi); 
+
+  if (Verbosity() > 1) std::cout << "photonjetskinematics::End(PHCompositeNode *topNode) This is the End..." << std::endl; 
+  return Fun4AllReturnCodes::EVENT_OK;
+
+}
+
+//____________________________________________________________________________..
+int photonjetskinematics::Reset(PHCompositeNode* /*topNode*/)
+{
+  if (Verbosity() > 1) std::cout << "photonjetskinematics::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
+  return Fun4AllReturnCodes::EVENT_OK;
+}
+
+//____________________________________________________________________________..
+void photonjetskinematics::Print(const std::string &what) const
+{
+  if (Verbosity() > 1) std::cout << "photonjetskinematics::Print(const std::string &what) const Printing info for " << what << std::endl;
+}

offline/QA/Jet/photonjetskinematics.h

@@ -0,0 +1,94 @@
+// Tell emacs that this is a C++ source
+//  -*- C++ -*-.
+#ifndef PHOTONJETSKINEMATICS_H
+#define PHOTONJETSKINEMATICS_H
+
+#include <fun4all/SubsysReco.h>
+
+#include <string>
+#include <TFile.h>
+#include <TH1.h>
+#include <TH2D.h>
+#include "JetQADefs.h"
+
+class Fun4AllHistoManager;
+class PHCompositeNode;
+class TH1;
+class TH2;
+class TH3;
+
+class photonjetskinematics : public SubsysReco
+{
+ public:
+
+
+  photonjetskinematics(const std::string &modulename = "photonjetskinematics", const std::string &inputnode = "CLUSTERINFO_CEMC");
+  ~photonjetskinematics() override;
+
+  /** Called during initialization.
+      Typically this is where you can book histograms, and e.g.
+      register them to Fun4AllServer (so they can be output to file
+      using Fun4AllServer::dumpHistos() method).
+   */
+  int Init(PHCompositeNode *topNode) override;
+
+  /** Called for first event when run number is known.
+      Typically this is where you may want to fetch data from
+      database, because you know the run number. A place
+      to book histograms which have to know the run number.
+   */
+  int InitRun(PHCompositeNode *topNode) override;
+
+  /** Called for each event.
+      This is where you do the real work.
+   */
+  int process_event(PHCompositeNode *topNode) override;
+
+  /// Clean up internals after each event.
+  int ResetEvent(PHCompositeNode *topNode) override;
+
+  /// Called at the end of each run.
+  int EndRun(const int runnumber) override;
+
+  /// Called at the end of all processing.
+  int End(PHCompositeNode *topNode) override;
+
+  /// Reset
+  int Reset(PHCompositeNode * /*topNode*/) override;
+
+  void Print(const std::string &what = "ALL") const override;
+
+  /// specifies a trigger to select
+  void SetTrgToSelect(const uint32_t trig = JetQADefs::GL1::MBDNSPhoton1)
+  {
+    doTrgSelect = true;
+    trgToSelect = trig;
+  }
+
+  /// set histogram tag
+  void SetHistTag(const std::string& tag)
+  {
+    histtag = tag;
+  }
+
+ private:
+  // std::string outfilename;
+  // TFile *outfile;
+ // hist manager
+
+ Fun4AllHistoManager* manager;
+ std::string modulename;
+ std::string inputnode;
+ std::string histtag;
+ uint32_t trgToSelect;
+ bool doTrgSelect;
+
+ ///Output histograms
+ TH1 *h_emcal_cluster_chi2 = nullptr;
+ TH1 *h_emcal_cluster_energy = nullptr;
+ TH2 *h_emcal_cluster_eta_phi = nullptr;
+ TH1 *h_emcal_cluster_eta = nullptr;    // <-- Declare eta histogram
+ TH1 *h_emcal_cluster_phi = nullptr;    // <-- Declare phi histogram
+};
+
+#endif // PHOTONJETSKINEMATICS_H