QA/Jet/RhosinEvent.cc

0001 #include "RhosinEvent.h"
0002
0003 // jetbackground includes
0004 #include <jetbackground/TowerRho.h>
0005 #include <jetbackground/EventRho.h>
0006
0007 // qautils include
0008 #include <qautils/QAHistManagerDef.h>
0009
0010 // fun4all includes
0011 #include <fun4all/Fun4AllHistoManager.h>
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013
0014 // phool includes
0015 #include <phool/PHCompositeNode.h>
0016 #include <phool/getClass.h>
0017
0018 #include <TH1.h>
0019 #include <TStyle.h>
0020
0021 #include <algorithm>
0022 #include <cassert>
0023 #include <cstdlib>
0024 #include <iostream>
0025 #include <string>
0026 #include <vector>
0027
0028 std::string to_lower(const std::string& s)
0029 {
0030   std::string out = s;
0031   std::transform(out.begin(), out.end(), out.begin(),
0032                  static_cast<int(*)(int)>(std::tolower));
0033   return out;
0034 }
0035
0036 RhosinEvent::RhosinEvent(const std::string& moduleName, const std::string& tag)
0037   : SubsysReco(moduleName)
0038   , m_moduleName(moduleName)
0039   , m_histTag(tag)
0040   // , m_name(outputfilename)
0041 {
0042 }
0043
0044 int RhosinEvent::Init(PHCompositeNode* /*topNode*/)
0045 {
0046   // if (Verbosity() > 0)
0047   // {
0048   // std::cout << "RhosinEvent::Init - Output to " << m_outputFileName << std::endl;
0049   // }
0050
0051   // create output file
0052   // PHTFileServer::get().open(m_outputFileName, "RECREATE");
0053   delete m_analyzer; // make cppcheck happy
0054   delete m_manager; // make cppcheck happy
0055   m_analyzer = new TriggerAnalyzer();
0056
0057   gStyle->SetOptTitle(0);
0058   m_manager = QAHistManagerDef::getHistoManager();
0059   if (!m_manager)
0060   {
0061     std::cout << PHWHERE << ": PANIC: couldn't grab histogram manager!" << std::endl;
0062     assert(m_manager);
0063   }
0064
0065   // Initialize histograms
0066   const int N_rho_mult = 320;
0067   const double rho_max_mult = 0.16;
0068   Double_t N_rho_mult_bins[N_rho_mult + 1];
0069   for (int i = 0; i <= N_rho_mult; i++)
0070   {
0071     N_rho_mult_bins[i] = (rho_max_mult / 320.0) * i;
0072   }
0073
0074   const int N_rho_area = 400;
0075   const double rho_max_area = 200;
0076   Double_t N_rho_area_bins[N_rho_area + 1];
0077   for (int i = 0; i <= N_rho_area; i++)
0078   {
0079     N_rho_area_bins[i] = (rho_max_area / 400.0) * i;
0080   }
0081
0082   // make sure module name is lower case
0083   std::string smallModuleName = m_moduleName;
0084   std::transform(
0085       smallModuleName.begin(),
0086       smallModuleName.end(),
0087       smallModuleName.begin(),
0088       ::tolower);
0089
0090   // construct histogram names
0091   std::vector<std::string> vecHistNames = {
0092       "rhomult",
0093       "sigmamult",
0094       "rhoarea",
0095       "sigmaarea"};
0096   for (auto& vecHistName : vecHistNames)
0097   {
0098     vecHistName.insert(0, "h_" + smallModuleName + "_");
0099     if (!m_histTag.empty())
0100     {
0101       vecHistName.append("_" + m_histTag);
0102     }
0103   }
0104
0105   h1_mult_rho = new TH1D(vecHistNames[0].data(), "", N_rho_mult, N_rho_mult_bins);
0106   h1_mult_rho->GetXaxis()->SetTitle("rho_M");
0107   h1_mult_rho->GetYaxis()->SetTitle("Counts");
0108
0109   h1_mult_rho_sigma = new TH1D(vecHistNames[1].data(), "", N_rho_mult, N_rho_mult_bins);
0110   h1_mult_rho_sigma->GetXaxis()->SetTitle("sigma_M");
0111   h1_mult_rho_sigma->GetYaxis()->SetTitle("Counts");
0112
0113   h1_area_rho = new TH1D(vecHistNames[2].data(), "", N_rho_area, N_rho_area_bins);
0114   h1_area_rho->GetXaxis()->SetTitle("rho_A");
0115   h1_area_rho->GetYaxis()->SetTitle("Counts");
0116
0117   h1_area_rho_sigma = new TH1D(vecHistNames[3].data(), "", N_rho_area, N_rho_area_bins);
0118   h1_area_rho_sigma->GetXaxis()->SetTitle("sigma_A");
0119   h1_area_rho_sigma->GetYaxis()->SetTitle("Counts");
0120
0121   if (Verbosity() > 0)
0122   {
0123     std::cout << "RhosinEvent::Init - Histograms created" << std::endl;
0124   }
0125
0126   return Fun4AllReturnCodes::EVENT_OK;
0127 }
0128
0129 int RhosinEvent::process_event(PHCompositeNode* topNode)
0130 {
0131   if (Verbosity() > 1)
0132   {
0133     std::cout << "RhosinEvent::process_event - Process event..." << std::endl;
0134   }
0135
0136   // if needed, check if selected trigger fired
0137   if (m_doTrgSelect)
0138   {
0139     m_analyzer->decodeTriggers(topNode);
0140     bool hasTrigger = JetQADefs::DidTriggerFire(m_trgToSelect, m_analyzer);
0141     if (!hasTrigger)
0142     {
0143       return Fun4AllReturnCodes::EVENT_OK;
0144     }
0145   }
0146
0147
0148   if (m_do_area_rho)
0149   {
0150
0151     TowerRho* towerrho = findNode::getClass<TowerRho>(topNode, m_area_rho_node);
0152     EventRho* eventrho = findNode::getClass<EventRho>(topNode, m_area_rho_node);
0153
0154     if (to_lower(m_area_rho_node).find("tower") != std::string::npos)
0155     {
0156         if (!towerrho)
0157         {
0158           std::cout << "RhosinEvent::process_event - Warning can not find towerrho " << m_area_rho_node << std::endl;
0159         }
0160         else
0161         {
0162           if (Verbosity() > 0)
0163           {
0164               std::cout << "RhosinEvent::process_event - Using towerrho " << m_area_rho_node << std::endl;
0165           }
0166           h1_area_rho->Fill(towerrho->get_rho());
0167           h1_area_rho_sigma->Fill(towerrho->get_sigma());
0168         }
0169     }
0170     else if (to_lower(m_area_rho_node).find("event") != std::string::npos)
0171     {
0172         if (!eventrho)
0173         {
0174           std::cout << "RhosinEvent::process_event - Warning can not find eventrho " << m_area_rho_node << std::endl;
0175         }
0176         else
0177         {
0178           if (Verbosity() > 0)
0179           {
0180               std::cout << "RhosinEvent::process_event - Using eventrho " << m_area_rho_node << std::endl;
0181           }
0182           h1_area_rho->Fill(eventrho->get_rho());
0183           h1_area_rho_sigma->Fill(eventrho->get_sigma());
0184         }
0185     }
0186     if (!eventrho && !towerrho)
0187     {
0188       std::cout << "RhosinEvent::process_event - Error can not find neither towerrho nor eventrho " << m_area_rho_node << std::endl;
0189       return Fun4AllReturnCodes::EVENT_OK;
0190     }
0191   }
0192
0193   if (m_do_mult_rho)
0194   {
0195     TowerRho* towerrho = findNode::getClass<TowerRho>(topNode, m_mult_rho_node);
0196     EventRho* eventrho = findNode::getClass<EventRho>(topNode, m_mult_rho_node);
0197
0198     if (to_lower(m_area_rho_node).find("tower") != std::string::npos)
0199     {
0200         if (!towerrho)
0201         {
0202           std::cout << "RhosinEvent::process_event - Warning can not find towerrho " << m_mult_rho_node << std::endl;
0203         }
0204         else
0205         {
0206           if (Verbosity() > 0)
0207           {
0208               std::cout << "RhosinEvent::process_event - Using towerrho " << m_mult_rho_node << std::endl;
0209           }
0210           h1_mult_rho->Fill(towerrho->get_rho());
0211           h1_mult_rho_sigma->Fill(towerrho->get_sigma());
0212         }
0213     }
0214     else if (to_lower(m_area_rho_node).find("event") != std::string::npos)
0215     {
0216         if (!eventrho)
0217         {
0218           std::cout << "RhosinEvent::process_event - Warning can not find eventrho " << m_mult_rho_node << std::endl;
0219         }
0220         else
0221         {
0222           if (Verbosity() > 0)
0223           {
0224               std::cout << "RhosinEvent::process_event - Using eventrho " << m_mult_rho_node << std::endl;
0225           }
0226           h1_mult_rho->Fill(eventrho->get_rho());
0227           h1_mult_rho_sigma->Fill(eventrho->get_sigma());
0228         }
0229     }
0230
0231     if (!eventrho && !towerrho)
0232     {
0233       std::cout << "RhosinEvent::process_event - Error can not find neither towerrho nor eventrho  " << m_mult_rho_node << std::endl;
0234       return Fun4AllReturnCodes::EVENT_OK;
0235     }
0236   }
0237
0238   if (Verbosity() > 1)
0239   {
0240     std::cout << "RhosinEvent::process_event - Event processed" << std::endl;
0241   }
0242
0243   return Fun4AllReturnCodes::EVENT_OK;
0244 }
0245
0246 int RhosinEvent::End(PHCompositeNode* /*topNode*/)
0247 {
0248   if (Verbosity() > 0)
0249   {
0250     std::cout << "RhosinEvent::EndRun - End run " << std::endl;
0251     // std::cout << "RhosinEvent::EndRun - Writing to " << m_outputFileName << std::endl;
0252   }
0253
0254   // PHTFileServer::get().cd(m_outputFileName);
0255
0256   // h1_mult_rho->Write();
0257   // h1_mult_rho_sigma->Write();
0258   // h1_area_rho->Write();
0259   // h1_area_rho_sigma->Write();
0260
0261   // register histograms
0262   m_manager->registerHisto(h1_mult_rho);
0263   m_manager->registerHisto(h1_mult_rho_sigma);
0264   m_manager->registerHisto(h1_area_rho);
0265   m_manager->registerHisto(h1_area_rho_sigma);
0266
0267   if (Verbosity() > 0)
0268   {
0269     std::cout << "RhosinEvent::EndRun - Done" << std::endl;
0270   }
0271
0272   return Fun4AllReturnCodes::EVENT_OK;
0273 }