calorimeter/calo_emc_noisy_tower/emcNoisyTowerFinder.cc

0001 #include "emcNoisyTowerFinder.h"
0002
0003 // Tower stuff
0004 #include <calobase/TowerInfo.h>
0005 #include <calobase/TowerInfoContainer.h>
0006 #include <calobase/TowerInfoDefs.h>
0007
0008 #include <cdbobjects/CDBTTree.h>
0009
0010 #include <ffamodules/CDBInterface.h>
0011
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <fun4all/SubsysReco.h>
0014
0015 #include <phool/getClass.h>
0016 #include <phool/phool.h>
0017
0018 // ROOT
0019 #include <TFile.h>
0020 #include <TH1.h>
0021 #include <TH2.h>
0022 #include <TProfile2D.h>
0023 #include <TSystem.h>
0024
0025 #include <algorithm>  // for sort
0026 #include <cstdlib>    // for exit, size_t
0027 #include <fstream>
0028 #include <iostream>
0029 #include <limits>  // for numeric_limits
0030
0031 //________________________________
0032 emcNoisyTowerFinder::emcNoisyTowerFinder(const std::string &name, const std::string &outputName)
0033   : SubsysReco(name)
0034   , Outfile(outputName)
0035 {
0036 }
0037
0038 //_____________________________
0039 int emcNoisyTowerFinder::InitRun(PHCompositeNode * /*topNode*/)
0040 {
0041   if (Verbosity() > 0)
0042   {
0043     std::cout << "emcNoisyTowerFinder::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0044   }
0045   foutput = new TFile(Outfile.c_str(), "recreate");
0046
0047   h_hits_eta_phi_adc = new TH2F("h_hits_eta_phi_adc", "", Neta, 0, Neta, Nphi, 0, Nphi);
0048   pr_hits_eta_phi_adc = new TProfile2D("pr_hits_eta_phi_adc", "", Neta, 0, Neta, Nphi, 0, Nphi);
0049
0050   h_hits_eta_phi_gev = new TH2F("h_hits_eta_phi_gev", "", Neta, 0, Neta, Nphi, 0, Nphi);
0051   pr_hits_eta_phi_gev = new TProfile2D("pr_hits_eta_phi_gev", "", Neta, 0, Neta, Nphi, 0, Nphi);
0052
0053   std::string default_time_independent_calib = "cemc_pi0_twrSlope_v1_default";
0054   m_fieldname = "Femc_datadriven_qm1_correction";
0055
0056   std::string calibdir = CDBInterface::instance()->getUrl(default_time_independent_calib);
0057   if (calibdir.empty())
0058   {
0059     std::cout << "CaloTowerCalib::::InitRun No EMCal Calibration NOT even a default" << std::endl;
0060     gSystem->Exit(1);
0061     exit(1);
0062   }
0063   cdbttree = new CDBTTree(calibdir);
0064
0065   return Fun4AllReturnCodes::EVENT_OK;
0066 }
0067
0068 //_____________________________
0069 int emcNoisyTowerFinder::process_event(PHCompositeNode *topNode)
0070 {
0071   // Get TowerInfoContainer
0072   TowerInfoContainer *towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_CEMC");
0073   if (!towers)
0074   {
0075     std::cout << PHWHERE << "emcNoisyTowerFinder::process_event Could not find node TOWERS_CEMC" << std::endl;
0076     return Fun4AllReturnCodes::ABORTEVENT;
0077   }
0078   // iterate through all towers, incrementing their Frequency arrays if they record a hit
0079
0080   int tower_range = towers->size();
0081   for (int j = 0; j < tower_range; j++)
0082   {
0083     TowerInfo *tower = towers->get_tower_at_channel(j);
0084     float energy = tower->get_energy();
0085     unsigned int towerkey = towers->encode_key(j);
0086     int ieta = towers->getTowerEtaBin(towerkey);
0087     int iphi = towers->getTowerPhiBin(towerkey);
0088
0089     float calibconst = cdbttree->GetFloatValue(towerkey, m_fieldname);
0090     float calib_energy = calibconst * energy;
0091
0092     if (energy > energy_threshold_adc)
0093     {
0094       h_hits_eta_phi_adc->Fill(ieta + 1, iphi + 1);
0095       pr_hits_eta_phi_adc->Fill(ieta + 1, iphi + 1, 1);
0096     }
0097     else
0098     {
0099       pr_hits_eta_phi_adc->Fill(ieta + 1, iphi + 1, 0);
0100     }
0101
0102     if (calib_energy > energy_threshold_gev)
0103     {
0104       h_hits_eta_phi_gev->Fill(ieta + 1, iphi + 1);
0105       pr_hits_eta_phi_gev->Fill(ieta + 1, iphi + 1, 1);
0106     }
0107     else
0108     {
0109       pr_hits_eta_phi_gev->Fill(ieta + 1, iphi + 1, 0);
0110     }
0111   }
0112
0113   return Fun4AllReturnCodes::EVENT_OK;
0114 }
0115
0116 //__________________________
0117 void emcNoisyTowerFinder::FindHot(std::string &infilename, std::string &outfilename, const std::string &inHist)
0118 {
0119   // TH2F *h_hits_eta_phi_adc = nullptr;
0120   bool isListFile = (infilename.substr(infilename.rfind('.') + 1) == "txt" || infilename.substr(infilename.rfind('.') + 1) == "list");
0121
0122   if (isListFile)
0123   {
0124     std::ifstream fileList(infilename);
0125     std::string line;
0126     while (std::getline(fileList, line))
0127     {
0128       TFile *fin = new TFile(line.c_str());
0129       if (!fin || fin->IsZombie())
0130       {
0131         std::cout << "emcNoisyTowerFinder::FindHot: input file not found or is corrupted " << line << std::endl;
0132         continue;
0133       }
0134       TH2 *tempHist{nullptr};
0135       fin->GetObject(inHist.c_str(), tempHist);
0136       if (!tempHist)
0137       {
0138         std::cout << "emcNoisyTowerFinder::FindHot: input hist not found in file " << line << std::endl;
0139         delete fin;
0140         continue;
0141       }
0142       if (!h_hits_eta_phi_adc)
0143       {
0144         h_hits_eta_phi_adc = (TH2 *) tempHist->Clone();
0145         h_hits_eta_phi_adc->SetDirectory(nullptr);  // Detach from the file to keep it in memory
0146       }
0147       else
0148       {
0149         h_hits_eta_phi_adc->Add(tempHist);
0150       }
0151       delete fin;
0152     }
0153     fileList.close();
0154   }
0155   else
0156   {
0157     TFile *fin = new TFile(infilename.c_str());
0158     if (!fin || fin->IsZombie())
0159     {
0160       std::cout << "emcNoisyTowerFinder::FindHot: input file not found or is corrupted " << infilename << std::endl;
0161       return;
0162     }
0163     fin->GetObject(inHist.c_str(), h_hits_eta_phi_adc);
0164     if (!h_hits_eta_phi_adc)
0165     {
0166       std::cout << "emcNoisyTowerFinder::FindHot: input hist not found " << inHist << std::endl;
0167       delete fin;
0168       return;
0169     }
0170     h_hits_eta_phi_adc->SetDirectory(nullptr);  // Detach from the file to keep it in memory
0171     delete fin;
0172   }
0173
0174   if (!h_hits_eta_phi_adc)
0175   {
0176     std::cout << "emcNoisyTowerFinder::FindHot: no valid histogram found" << std::endl;
0177     return;
0178   }
0179   TH2 *h_hits = h_hits_eta_phi_adc;
0180
0181   TFile *fout = new TFile(outfilename.c_str(), "recreate");
0182
0183   TH2 *h_hot = new TH2F("h_hot", "", Neta, 0, Neta, Nphi, 0, Nphi);
0184   TH2 *h_hitClean{nullptr};  // = new TH2F("h_hitClean", "", Neta, 0, Neta, Nphi, 0, Nphi);
0185   TH2 *h_heatSigma = new TH2F("h_heatSigma", "", Neta, 0, Neta, Nphi, 0, Nphi);
0186   TH1 *h_perMedian = new TH1F("h_perMedian", "", 500, 0, 5);
0187   std::vector<TH1 *> h1_hits;
0188   h1_hits.resize(Neta);
0189   std::vector<TH1 *> h1_hits2;
0190   h1_hits2.resize(Neta);
0191   h_hits->Write();
0192   h_hitClean = (TH2F *) h_hits->Clone("h_hitClean");
0193
0194   float max = h_hits->GetBinContent(h_hits->GetMaximumBin());
0195   float min = h_hits->GetMinimum();
0196   if (min == 0)
0197   {
0198     min = 1;
0199   }
0200
0201   for (int ie = 0; ie < Neta; ie++)
0202   {
0203     h1_hits[ie] = new TH1F((std::string("h1_hits") + std::to_string(ie)).c_str(), "", 200, min, max);
0204     h1_hits2[ie] = new TH1F((std::string("h1_hits2_") + std::to_string(ie)).c_str(), "", 200, min, max);
0205     std::vector<float> vals;
0206     for (int iphi = 0; iphi < Nphi; iphi++)
0207     {
0208       float val = h_hits->GetBinContent(ie + 1, iphi + 1);
0209       h1_hits[ie]->Fill(val);
0210       vals.push_back(val);
0211     }
0212
0213     float median = findMedian(vals);
0214     float mean = h1_hits[ie]->GetMean();
0215     float std = h1_hits[ie]->GetStdDev();
0216     for (int iphi = 0; iphi < Nphi; iphi++)
0217     {
0218       float val = h_hits->GetBinContent(ie + 1, iphi + 1);
0219       if (val / median > 5 || val / median < 0.1)
0220       {
0221         continue;
0222       }
0223       h1_hits2[ie]->Fill(val);
0224     }
0225     mean = h1_hits2[ie]->GetMean();
0226     std = h1_hits2[ie]->GetStdDev();
0227     for (int iphi = 0; iphi < Nphi; iphi++)
0228     {
0229       h_hot->SetBinContent(ie + 1, iphi + 1, 0);
0230       float val = h_hits->GetBinContent(ie + 1, iphi + 1);
0231       h_perMedian->Fill(val / median);
0232       float sigma = (val - mean) / std;
0233       h_heatSigma->SetBinContent(ie + 1, iphi + 1, sigma);
0234       if (sigma > sigma_bad_thresh)  // hot tower
0235       {
0236         h_hot->SetBinContent(ie + 1, iphi + 1, 2);
0237         h_hitClean->SetBinContent(ie + 1, iphi + 1, 0);
0238       }
0239       if ((sigma < -1 * sigma_bad_thresh || val < mean * percent_cold_thresh) && val != 0)  // cold tower
0240       {
0241         h_hot->SetBinContent(ie + 1, iphi + 1, 3);
0242         h_hitClean->SetBinContent(ie + 1, iphi + 1, 0);
0243       }
0244       if (val == 0)  // dead tower
0245       {
0246         h_hot->SetBinContent(ie + 1, iphi + 1, 1);
0247       }
0248     }
0249   }
0250   h_hitClean->Write();
0251
0252   fout->Write();
0253
0254   ////////////////////////////////////////////
0255   // make cdb tree
0256   size_t pos = outfilename.find_last_of('.');
0257   std::string f_cdbout_name = outfilename;
0258   f_cdbout_name.insert(pos, "cdb");
0259
0260   CDBTTree *cdbttree_out = new CDBTTree(f_cdbout_name);
0261   std::string m_fieldname_out = "status";
0262   for (int i = 0; i < Neta; i++)
0263   {
0264     for (int j = 0; j < Nphi; j++)
0265     {
0266       unsigned int key = TowerInfoDefs::encode_emcal(i, j);
0267       if (Neta == 24)
0268       {
0269         key = TowerInfoDefs::encode_hcal(i, j);
0270       }
0271       int val = h_hot->GetBinContent(i + 1, j + 1);
0272       float sigma = h_heatSigma->GetBinContent(i + 1, j + 1);
0273       cdbttree_out->SetIntValue(key, m_fieldname_out, val);
0274       cdbttree_out->SetFloatValue(key, m_caloName + "_sigma", sigma);
0275     }
0276   }
0277   cdbttree_out->Commit();
0278   cdbttree_out->WriteCDBTTree();
0279   delete cdbttree_out;
0280
0281   fout->Close();
0282 }
0283
0284 //____________________________________________________
0285 int emcNoisyTowerFinder::End(PHCompositeNode * /*topNode*/)
0286 {
0287   foutput->Write();
0288   // fchannels -> cd();
0289   // fchannels->Close();
0290   // delete fchannels;
0291   // fchannels=NULL;
0292
0293   std::cout << "emcNoisyTowerFinder::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0294   return Fun4AllReturnCodes::EVENT_OK;
0295 }
0296 //____________________________________________________
0297
0298 float emcNoisyTowerFinder::findMedian(const std::vector<float> &arr)
0299 {
0300   if (arr.empty())
0301   {
0302     return std::numeric_limits<float>::quiet_NaN();  // Return NaN if the array is empty
0303   }
0304   std::vector<float> sortedArr = arr;
0305   std::sort(sortedArr.begin(), sortedArr.end());
0306
0307   size_t n = sortedArr.size();
0308   if (n % 2 == 0)
0309   {
0310     return (sortedArr[(n / 2) - 1] + sortedArr[n / 2]) / 2.0;
0311   }
0312
0313   return sortedArr[n / 2];
0314 }