EMCal-Hot-Towers/src/CaloHotTower.cc

0001 //____________________________________________________________________________..
0002 //
0003 // This is a template for a Fun4All SubsysReco module with all methods from the
0004 // $OFFLINE_MAIN/include/fun4all/SubsysReco.h baseclass
0005 // You do not have to implement all of them, you can just remove unused methods
0006 // here and in CaloHotTower.h.
0007 //
0008 // CaloHotTower(const std::string &name = "CaloHotTower")
0009 // everything is keyed to CaloHotTower, duplicate names do work but it makes
0010 // e.g. finding culprits in logs difficult or getting a pointer to the module
0011 // from the command line
0012 //
0013 // CaloHotTower::~CaloHotTower()
0014 // this is called when the Fun4AllServer is deleted at the end of running. Be
0015 // mindful what you delete - you do loose ownership of object you put on the node tree
0016 //
0017 // int CaloHotTower::Init(PHCompositeNode *topNode)
0018 // This method is called when the module is registered with the Fun4AllServer. You
0019 // can create historgrams here or put objects on the node tree but be aware that
0020 // modules which haven't been registered yet did not put antyhing on the node tree
0021 //
0022 // int CaloHotTower::InitRun(PHCompositeNode *topNode)
0023 // This method is called when the first event is read (or generated). At
0024 // this point the run number is known (which is mainly interesting for raw data
0025 // processing). Also all objects are on the node tree in case your module's action
0026 // depends on what else is around. Last chance to put nodes under the DST Node
0027 // We mix events during readback if branches are added after the first event
0028 //
0029 // int CaloHotTower::process_event(PHCompositeNode *topNode)
0030 // called for every event. Return codes trigger actions, you find them in
0031 // $OFFLINE_MAIN/include/fun4all/Fun4AllReturnCodes.h
0032 //   everything is good:
0033 //     return Fun4AllReturnCodes::EVENT_OK
0034 //   abort event reconstruction, clear everything and process next event:
0035 //     return Fun4AllReturnCodes::ABORT_EVENT;
0036 //   proceed but do not save this event in output (needs output manager setting):
0037 //     return Fun4AllReturnCodes::DISCARD_EVENT;
0038 //   abort processing:
0039 //     return Fun4AllReturnCodes::ABORT_RUN
0040 // all other integers will lead to an error and abort of processing
0041 //
0042 // int CaloHotTower::ResetEvent(PHCompositeNode *topNode)
0043 // If you have internal data structures (arrays, stl containers) which needs clearing
0044 // after each event, this is the place to do that. The nodes under the DST node are cleared
0045 // by the framework
0046 //
0047 // int CaloHotTower::EndRun(const int runnumber)
0048 // This method is called at the end of a run when an event from a new run is
0049 // encountered. Useful when analyzing multiple runs (raw data). Also called at
0050 // the end of processing (before the End() method)
0051 //
0052 // int CaloHotTower::End(PHCompositeNode *topNode)
0053 // This is called at the end of processing. It needs to be called by the macro
0054 // by Fun4AllServer::End(), so do not forget this in your macro
0055 //
0056 // int CaloHotTower::Reset(PHCompositeNode *topNode)
0057 // not really used - it is called before the dtor is called
0058 //
0059 // void CaloHotTower::Print(const std::string &what) const
0060 // Called from the command line - useful to print information when you need it
0061 //
0062 //____________________________________________________________________________..
0063
0064 #include "CaloHotTower.h"
0065
0066 #include <fun4all/Fun4AllReturnCodes.h>
0067 #include <fun4all/Fun4AllServer.h>
0068
0069 #include <phool/PHCompositeNode.h>
0070 #include <phool/getClass.h>
0071 #include <phool/phool.h>
0072
0073 // Tower stuff
0074 #include <calobase/TowerInfo.h>
0075 #include <calobase/TowerInfoDefs.h>
0076 #include <calobase/TowerInfoContainer.h>
0077
0078 // CDB
0079 #include <cdbobjects/CDBTTree.h>
0080 #include <ffamodules/CDBInterface.h>
0081
0082 // ROOT stuff
0083 #include <TFile.h>
0084 #include <TH1.h>
0085
0086 // c++
0087 #include <iostream>
0088 #include <fstream>
0089
0090 using std::cout;
0091 using std::cerr;
0092 using std::endl;
0093 using std::string;
0094 using std::vector;
0095 using std::pair;
0096 using std::min;
0097 using std::max;
0098 using std::to_string;
0099
0100 //____________________________________________________________________________..
0101 CaloHotTower::CaloHotTower(const string &name):
0102  SubsysReco(name),
0103  iEvent(0),
0104  energy_min(9999),
0105  energy_max(0),
0106  bins_energy(3280),
0107  energy_low(0),
0108  energy_high(16400), // 2^14 is the max value
0109  bins_events(m_triggers.size()),
0110  m_emcTowerNode("TOWERS_CEMC"),
0111  m_outputFile("test.root"),
0112  m_calibName_hotMap("CEMC_BadTowerMap"),
0113  triggeranalyzer(nullptr)
0114 {
0115   cout << "CaloHotTower::CaloHotTower(const string &name) Calling ctor" << endl;
0116 }
0117
0118 //____________________________________________________________________________..
0119 CaloHotTower::~CaloHotTower() {
0120   cout << "CaloHotTower::~CaloHotTower() Calling dtor" << endl;
0121 }
0122
0123 Int_t CaloHotTower::readHotTowerIndexFile() {
0124
0125   cout << "Reading Hot and Ref tower indices" << endl;
0126
0127   std::ifstream file(m_hotTowerIndexFile);
0128
0129   // Check if the file was successfully opened
0130   if (!file.is_open()) {
0131       cerr << "Failed to open file: " << m_hotTowerIndexFile << endl;
0132       return 1;
0133   }
0134
0135   string line;
0136
0137   // skip header
0138   std::getline(file, line);
0139
0140   // loop over each run
0141   while (std::getline(file, line)) {
0142       std::istringstream lineStream(line);
0143
0144       UInt_t towerHotIndex;
0145       UInt_t towerRefIndex;
0146       Char_t comma;
0147
0148       if (lineStream >> towerHotIndex >> comma >> towerRefIndex) {
0149           hotTowerIndex.push_back(std::make_pair(towerHotIndex, towerRefIndex));
0150           cout << "Hot: " << towerHotIndex << ", Ref: " << towerRefIndex << endl;
0151       }
0152       else {
0153           cerr << "Failed to parse line: " << line << endl;
0154           return 1;
0155       }
0156   }
0157
0158   // Close the file
0159   file.close();
0160
0161   return 0;
0162 }
0163
0164 //____________________________________________________________________________..
0165 Int_t CaloHotTower::Init(PHCompositeNode *topNode) {
0166   cout << "CaloHotTower::Init(PHCompositeNode *topNode) Initializing" << endl;
0167
0168   // read the hot tower and reference tower indices
0169   Int_t ret = readHotTowerIndexFile();
0170
0171   if(ret) return Fun4AllReturnCodes::ABORTRUN;
0172
0173   // so that the histos actually get written out
0174   Fun4AllServer *se = Fun4AllServer::instance();
0175   se->Print("NODETREE");
0176
0177   hEvents = new TH1F("hEvents","Events; Status; Counts", bins_events, 0, bins_events);
0178   for(UInt_t i = 1; i <= bins_events; ++i) {
0179     hEvents->GetXaxis()->SetBinLabel(i, m_triggers[i-1].c_str());
0180   }
0181
0182   triggeranalyzer = new TriggerAnalyzer();
0183
0184   UInt_t i = 0;
0185   // initialize histograms
0186   for(auto idx : hotTowerIndex) {
0187
0188     UInt_t key    = TowerInfoDefs::encode_emcal(idx.first);
0189     UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
0190     UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
0191
0192     string name  = "HotTower_"+to_string(phibin)+"_"+to_string(etabin);
0193     string title = "Hot Tower: iphi: " + to_string(phibin) + ", ieta: " + to_string(etabin) + "; ADC; Counts";
0194     auto h = new TH1F(name.c_str(), title.c_str(), bins_energy, energy_low, energy_high);
0195     hHotTowerEnergy.push_back(h);
0196
0197     name  = "HotTowerComplement_"+to_string(phibin)+"_"+to_string(etabin);
0198     title = "Hot Tower Complement: iphi: " + to_string(phibin) + ", ieta: " + to_string(etabin) + "; ADC; Counts";
0199     h = new TH1F(name.c_str(), title.c_str(), bins_energy, energy_low, energy_high);
0200     hHotTowerComplementEnergy.push_back(h);
0201
0202     key    = TowerInfoDefs::encode_emcal(idx.second);
0203     etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
0204     phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
0205
0206     name  = "RefTower_"+to_string(phibin)+"_"+to_string(etabin)+"_"+to_string(i);
0207     title = "Reference Tower: iphi: " + to_string(phibin) + ", ieta: " + to_string(etabin) + "; ADC; Counts";
0208     h = new TH1F(name.c_str(), title.c_str(), bins_energy, energy_low, energy_high);
0209     hRefTowerEnergy.push_back(h);
0210
0211     ++i;
0212   }
0213
0214   string calibdir = CDBInterface::instance()->getUrl(m_calibName_hotMap);
0215
0216   if (calibdir.empty()) {
0217       cout << PHWHERE << "CaloHotTower::Init Could not find " << m_calibName_hotMap << endl;
0218       return Fun4AllReturnCodes::ABORTRUN;
0219   }
0220
0221   m_cdbttree_hotMap = new CDBTTree(calibdir);
0222
0223   // print used DB files
0224   CDBInterface::instance()->Print();
0225
0226   return Fun4AllReturnCodes::EVENT_OK;
0227 }
0228
0229 //____________________________________________________________________________..
0230 Int_t CaloHotTower::process_event(PHCompositeNode *topNode) {
0231
0232   if(iEvent%20 == 0) cout << "Progress: " << iEvent << endl;
0233   ++iEvent;
0234
0235   // Get TowerInfoContainer
0236   TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, m_emcTowerNode.c_str());
0237   if (!towers) {
0238     cout << PHWHERE << "CaloHotTower::process_event Could not find node " << m_emcTowerNode << endl;
0239     return Fun4AllReturnCodes::ABORTEVENT;
0240   }
0241
0242   triggeranalyzer->decodeTriggers(topNode);
0243
0244   Bool_t isMB = false;
0245   for(UInt_t i = 0; i < m_triggers.size(); ++i) {
0246     if(triggeranalyzer->didTriggerFire(m_triggers[i])) {
0247       isMB = true;
0248       hEvents->Fill(i);
0249     }
0250   }
0251
0252   // skip event if MBD trigger isn't fired
0253   if(!isMB) return Fun4AllReturnCodes::ABORTEVENT;
0254
0255   for (UInt_t i = 0; i < hotTowerIndex.size(); ++i) {
0256     UInt_t towerIndex = hotTowerIndex[i].first;
0257     UInt_t key        = towers->encode_key(towerIndex);
0258     Int_t hotMap_val  = m_cdbttree_hotMap->GetIntValue(key, "status");
0259     TowerInfo* tower  = towers->get_tower_at_channel(towerIndex);
0260     Float_t energy    = tower->get_energy();
0261     Bool_t isBadChi2  = tower->get_isBadChi2();
0262
0263     // status 2: hot tower
0264     if(hotMap_val == 2 && !isBadChi2) {
0265       hHotTowerEnergy[i]->Fill(energy);
0266       energy_min = min(energy_min, energy);
0267       energy_max = max(energy_max, energy);
0268     }
0269
0270     if(hotMap_val != 2 && !isBadChi2) {
0271       hHotTowerComplementEnergy[i]->Fill(energy);
0272       energy_min = min(energy_min, energy);
0273       energy_max = max(energy_max, energy);
0274     }
0275
0276     towerIndex = hotTowerIndex[i].second;
0277     key        = towers->encode_key(towerIndex);
0278     hotMap_val = m_cdbttree_hotMap->GetIntValue(key, "status");
0279     tower      = towers->get_tower_at_channel(towerIndex);
0280     energy     = tower->get_energy();
0281     isBadChi2  = tower->get_isBadChi2();
0282
0283     // status 0: good tower
0284     if(!hotMap_val && !isBadChi2) {
0285       hRefTowerEnergy[i]->Fill(energy);
0286       energy_min = min(energy_min, energy);
0287       energy_max = max(energy_max, energy);
0288     }
0289   }
0290
0291   return Fun4AllReturnCodes::EVENT_OK;
0292 }
0293
0294 //____________________________________________________________________________..
0295 Int_t CaloHotTower::End(PHCompositeNode *topNode) {
0296   cout << "CaloHotTower::End(PHCompositeNode *topNode) This is the End..." << endl;
0297   cout << "Min Energy: " << energy_min << ", Max Energy: " << energy_max << endl;
0298   cout << "Trigger Summary" << endl;
0299   for(UInt_t i = 0; i < m_triggers.size(); ++i) {
0300     cout << m_triggers[i] << ": " << hEvents->GetBinContent(i+1) << " Events" << endl;
0301   }
0302
0303   TFile output(m_outputFile.c_str(),"recreate");
0304
0305   output.cd();
0306   hEvents->Write();
0307
0308   output.mkdir("Hot");
0309   output.mkdir("HotComplement");
0310   output.mkdir("Ref");
0311
0312   for(UInt_t i = 0; i < hotTowerIndex.size(); ++i) {
0313     output.cd("Hot");
0314     hHotTowerEnergy[i]->Write();
0315
0316     output.cd("HotComplement");
0317     hHotTowerComplementEnergy[i]->Write();
0318
0319     output.cd("Ref");
0320     hRefTowerEnergy[i]->Write();
0321   }
0322
0323   output.Close();
0324
0325   return Fun4AllReturnCodes::EVENT_OK;
0326 }