File indexing completed on 2025-08-05 08:11:14
0001 #include "RandomConeAna.h"
0002
0003 #include <g4main/PHG4TruthInfoContainer.h>
0004 #include <g4main/PHG4Particle.h>
0005 #include <TMath.h>
0006 #include <calotrigger/TriggerRunInfov1.h>
0007 #include <calotrigger/TriggerAnalyzer.h>
0008 #include "TRandom.h"
0009 #include "TH1D.h"
0010 #include "TProfile.h"
0011
0012 #include <calobase/RawTowerGeom.h>
0013 #include <jetbackground/TowerBackground.h>
0014 #include <calobase/RawCluster.h>
0015 #include <calobase/RawClusterContainer.h>
0016 #include <calobase/RawClusterUtility.h>
0017 #include <calobase/RawTowerGeomContainer.h>
0018 #include <calobase/RawTower.h>
0019 #include <calobase/RawTowerContainer.h>
0020 #include <HepMC/SimpleVector.h>
0021
0022
0023 #include <vector>
0024 #include <fun4all/Fun4AllHistoManager.h>
0025 #include <fun4all/Fun4AllReturnCodes.h>
0026 #include <phool/PHCompositeNode.h>
0027 #include <phool/PHIODataNode.h>
0028 #include <phool/PHNode.h>
0029 #include <phool/PHNodeIterator.h>
0030 #include <phool/PHObject.h>
0031 #include <phool/getClass.h>
0032
0033
0034 #include <iostream>
0035
0036 #include <map>
0037
0038
0039 RandomConeAna::RandomConeAna(const std::string &name, const std::string &outfilename):
0040 SubsysReco(name)
0041
0042 {
0043 isSim = false;
0044 _foutname = outfilename;
0045 m_calo_nodename = "TOWERINFO_CALIB";
0046 m_mb_triggerlist[0] = "MBD N&S >= 1";
0047
0048 m_mb_triggerlist[1] = "MBD N&S >= 1, vtx < 10 cm";
0049 }
0050
0051
0052 RandomConeAna::~RandomConeAna()
0053 {
0054
0055 }
0056
0057
0058 int RandomConeAna::Init(PHCompositeNode *topNode)
0059 {
0060
0061 rdm = new TRandom();
0062 hm = new Fun4AllHistoManager("RANDOMCONEHISTOS");
0063
0064 triggeranalyzer = new TriggerAnalyzer();
0065
0066 h_random_et = new TH1D("h_random_et", ";eT;", 200, -10, 10);
0067 h_random_emcal_et = new TH1D("h_random_emcal_et", ";eT;", 200, -10, 10);
0068 h_random_emcalre_et = new TH1D("h_random_emcalre_et", ";eT;", 200, -10, 10);
0069 h_random_hcalin_et = new TH1D("h_random_hcalin_et", ";eT;", 200, -10, 10);
0070 h_random_hcalout_et = new TH1D("h_random_hcalout_et", ";eT;", 200, -10, 10);
0071
0072
0073 hp_random_emcal_fraction = new TProfile("hp_random_emcal_fraction", ";et;", 200, -10, 10);
0074 hp_random_emcalre_fraction = new TProfile("hp_random_emcalre_fraction", ";et;", 200, -10, 10);
0075 hp_random_hcalin_fraction = new TProfile("hp_random_hcalin_fraction", ";eT;", 200, -10, 10);
0076 hp_random_hcalout_fraction = new TProfile("hp_random_hcalout_fraction", ";eT;", 200, -10, 10);
0077 hm->registerHisto(h_random_et);
0078 hm->registerHisto(h_random_emcal_et);
0079 hm->registerHisto(h_random_emcalre_et);
0080 hm->registerHisto(h_random_hcalin_et);
0081 hm->registerHisto(h_random_hcalout_et);
0082 hm->registerHisto(hp_random_emcal_fraction);
0083 hm->registerHisto(hp_random_emcalre_fraction);
0084 hm->registerHisto(hp_random_hcalin_fraction);
0085 hm->registerHisto(hp_random_hcalout_fraction);
0086
0087
0088 return Fun4AllReturnCodes::EVENT_OK;
0089 }
0090
0091
0092 int RandomConeAna::InitRun(PHCompositeNode *topNode)
0093 {
0094 if (Verbosity()) std::cout << __FUNCTION__ << __LINE__<<std::endl;
0095 return Fun4AllReturnCodes::EVENT_OK;
0096 }
0097
0098
0099
0100 void RandomConeAna::SetVerbosity(int verbo){
0101 _verbosity = verbo;
0102 return;
0103 }
0104
0105 void RandomConeAna::reset_tree_vars()
0106 {
0107 if (Verbosity()) std::cout << __FUNCTION__ << __LINE__<<std::endl;
0108
0109
0110 return;
0111 }
0112
0113 int RandomConeAna::process_event(PHCompositeNode *topNode)
0114 {
0115
0116 if (Verbosity()) std::cout << __FILE__ << " "<< __LINE__<<" "<<std::endl;
0117
0118 _i_event++;
0119
0120 if (!isSim)
0121 {
0122 triggeranalyzer->decodeTriggers(topNode);
0123
0124 if (!triggeranalyzer->didTriggerFire(m_mb_triggerlist[0]) && !triggeranalyzer->didTriggerFire(m_mb_triggerlist[1]))
0125 {
0126 return Fun4AllReturnCodes::EVENT_OK;
0127 }
0128 }
0129 RawTowerGeomContainer *tower_geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0130 RawTowerGeomContainer *tower_geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0131 RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0132
0133
0134 TowerInfoContainer *hcalin_towers = findNode::getClass<TowerInfoContainer>(topNode, m_calo_nodename + "_HCALIN");
0135 if (!hcalin_towers)
0136 {
0137 if (Verbosity()) std::cout << "no hcalin towers "<<std::endl;
0138 return Fun4AllReturnCodes::ABORTRUN;
0139 }
0140
0141 TowerInfoContainer *hcalout_towers = findNode::getClass<TowerInfoContainer>(topNode, m_calo_nodename + "_HCALOUT");
0142 if (!hcalout_towers)
0143 {
0144 std::cout << "no hcalout towers "<<std::endl;
0145 return Fun4AllReturnCodes::ABORTRUN;
0146 }
0147 TowerInfoContainer *emcalre_towers = findNode::getClass<TowerInfoContainer>(topNode, m_calo_nodename + "_CEMC_RETOWER");
0148 TowerInfoContainer *emcal_towers = findNode::getClass<TowerInfoContainer>(topNode, m_calo_nodename + "_CEMC");
0149
0150 float jet_phi = rdm->Uniform(-1*TMath::Pi(), TMath::Pi());
0151 float jet_eta = rdm->Uniform(-0.7, 0.7);
0152
0153 int size;
0154
0155 double emcalre_eT = 0;
0156 size = emcalre_towers->size();
0157 for (int i = 0; i < size; i++)
0158 {
0159 TowerInfo* tower = emcalre_towers->get_tower_at_channel(i);
0160
0161 if (!tower || !tower_geomIH)
0162 {
0163 continue;
0164 }
0165 if(!tower->get_isGood()) continue;
0166
0167 unsigned int calokey = emcalre_towers->encode_key(i);
0168
0169 int ieta = emcalre_towers->getTowerEtaBin(calokey);
0170 int iphi = emcalre_towers->getTowerPhiBin(calokey);
0171 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0172 float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0173 float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0174 float dR = sqrt(TMath::Power(jet_eta - tower_eta, 2) + TMath::Power(jet_phi - tower_phi, 2));
0175 if (dR >= 0.4)
0176 {
0177 continue;
0178 }
0179 double tower_eT = tower->get_energy() / std::cosh(tower_eta);
0180 emcalre_eT += tower_eT;
0181 }
0182 double emcal_eT = 0;
0183 size = emcal_towers->size();
0184 for (int i = 0; i < size; i++)
0185 {
0186 TowerInfo* tower = emcal_towers->get_tower_at_channel(i);
0187
0188 if (!tower || !tower_geomEM)
0189 {
0190 continue;
0191 }
0192 if(!tower->get_isGood()) continue;
0193
0194 unsigned int calokey = emcal_towers->encode_key(i);
0195
0196 int ieta = emcal_towers->getTowerEtaBin(calokey);
0197 int iphi = emcal_towers->getTowerPhiBin(calokey);
0198 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta, iphi);
0199 float tower_phi = tower_geomEM->get_tower_geometry(key)->get_phi();
0200 float tower_eta = tower_geomEM->get_tower_geometry(key)->get_eta();
0201 float dR = sqrt(TMath::Power(jet_eta - tower_eta, 2) + TMath::Power(jet_phi - tower_phi, 2));
0202 if (dR >= 0.4)
0203 {
0204 continue;
0205 }
0206 double tower_eT = tower->get_energy() / std::cosh(tower_eta);
0207 emcal_eT += tower_eT;
0208 }
0209 double hcalin_eT = 0;
0210 size = hcalin_towers->size();
0211 for (int i = 0; i < size; i++)
0212 {
0213 TowerInfo* tower = hcalin_towers->get_tower_at_channel(i);
0214
0215 if (!tower || !tower_geomIH)
0216 {
0217 continue;
0218 }
0219 if(!tower->get_isGood()) continue;
0220
0221 unsigned int calokey = hcalin_towers->encode_key(i);
0222
0223 int ieta = hcalin_towers->getTowerEtaBin(calokey);
0224 int iphi = hcalin_towers->getTowerPhiBin(calokey);
0225 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0226 float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0227 float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0228 float dR = sqrt(TMath::Power(jet_eta - tower_eta, 2) + TMath::Power(jet_phi - tower_phi, 2));
0229 if (dR >= 0.4)
0230 {
0231 continue;
0232 }
0233 double tower_eT = tower->get_energy() / std::cosh(tower_eta);
0234 hcalin_eT += tower_eT;
0235 }
0236 double hcalout_eT = 0;
0237 size = hcalout_towers->size();
0238 for (int i = 0; i < size; i++)
0239 {
0240 TowerInfo* tower = hcalout_towers->get_tower_at_channel(i);
0241
0242 if (!tower || !tower_geomIH)
0243 {
0244 continue;
0245 }
0246 if(!tower->get_isGood()) continue;
0247
0248 unsigned int calokey = hcalout_towers->encode_key(i);
0249
0250 int ieta = hcalout_towers->getTowerEtaBin(calokey);
0251 int iphi = hcalout_towers->getTowerPhiBin(calokey);
0252 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0253 float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0254 float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0255 float dR = sqrt(TMath::Power(jet_eta - tower_eta, 2) + TMath::Power(jet_phi - tower_phi, 2));
0256 if (dR >= 0.4)
0257 {
0258 continue;
0259 }
0260 double tower_eT = tower->get_energy() / std::cosh(tower_eta);
0261 hcalout_eT += tower_eT;
0262 }
0263
0264
0265 double totalre_eT = emcalre_eT + hcalin_eT + hcalout_eT;
0266 double total_eT = emcal_eT + hcalin_eT + hcalout_eT;
0267
0268 h_random_et->Fill(totalre_eT);
0269 h_random_emcalre_et->Fill(emcalre_eT);
0270 h_random_emcal_et->Fill(emcal_eT);
0271 h_random_hcalin_et->Fill(hcalin_eT);
0272 h_random_hcalout_et->Fill(hcalout_eT);
0273
0274 hp_random_emcalre_fraction->Fill(totalre_eT, emcalre_eT/totalre_eT);
0275 hp_random_emcal_fraction->Fill(total_eT,emcal_eT/total_eT);
0276 hp_random_hcalin_fraction->Fill(totalre_eT,hcalin_eT/totalre_eT);
0277 hp_random_hcalout_fraction->Fill(totalre_eT,hcalout_eT/totalre_eT);
0278 return Fun4AllReturnCodes::EVENT_OK;
0279 }
0280
0281
0282
0283 void RandomConeAna::GetNodes(PHCompositeNode* topNode)
0284 {
0285
0286
0287 }
0288
0289 int RandomConeAna::ResetEvent(PHCompositeNode *topNode)
0290 {
0291 if (Verbosity() > 0)
0292 {
0293 std::cout << "RandomConeAna::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0294 }
0295
0296
0297 return Fun4AllReturnCodes::EVENT_OK;
0298 }
0299
0300
0301 int RandomConeAna::EndRun(const int runnumber)
0302 {
0303 if (Verbosity() > 0)
0304 {
0305 std::cout << "RandomConeAna::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0306 }
0307
0308 return Fun4AllReturnCodes::EVENT_OK;
0309 }
0310
0311
0312 int RandomConeAna::End(PHCompositeNode *topNode)
0313 {
0314 if (Verbosity() > 0)
0315 {
0316 std::cout << "RandomConeAna::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0317 }
0318 hm->dumpHistos(_foutname.c_str());
0319
0320 return Fun4AllReturnCodes::EVENT_OK;
0321 }
0322
0323
0324 int RandomConeAna::Reset(PHCompositeNode *topNode)
0325 {
0326 if (Verbosity() > 0)
0327 {
0328 std::cout << "RandomConeAna::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0329 }
0330 return Fun4AllReturnCodes::EVENT_OK;
0331 }
0332
