HCalHotTowerFinder/src/hcal_towerid.cc

0001 #include "hcal_towerid.h"
0002 #include <fun4all/Fun4AllReturnCodes.h>
0003 #include <phool/PHCompositeNode.h>
0004
0005 //Fun4All
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007 #include <fun4all/Fun4AllServer.h>
0008 #include <fun4all/Fun4AllHistoManager.h>
0009 #include <phool/PHCompositeNode.h>
0010 #include <phool/getClass.h>
0011 #include <phool/phool.h>
0012 #include <ffaobjects/EventHeader.h>
0013
0014 //ROOT
0015 #include <TFile.h>
0016 #include <TTree.h>
0017 #include <TH1F.h>
0018 #include <TH2F.h>
0019
0020 //Tower stuff
0021 #include <calobase/TowerInfoContainer.h>
0022 #include <calobase/TowerInfo.h>
0023 #include <calobase/TowerInfoDefs.h>
0024
0025 #include <cdbobjects/CDBTTree.h>
0026
0027 //________________________________
0028 hcal_towerid::hcal_towerid(const std::string &name, const std::string &cdbtreename_i, const std::string &cdbtreename_o, float adccut_i, float adccut_o, float sigmas_lo, float sigmas_hi, float inner_f,float outer_f):
0029 SubsysReco(name)
0030   ,T(nullptr)
0031   ,Outfile(name)
0032   ,cdbtreename_i(cdbtreename_i)
0033   ,cdbtreename_o(cdbtreename_o)
0034   ,adccut_i(adccut_i)
0035   ,adccut_o(adccut_o)       // Minimum ADC for Kurary fibers to register a hit
0036   ,sigmas_lo(sigmas_lo)     // # of standard deviations from the mode for a cold tower
0037   ,sigmas_hi(sigmas_hi)     // # of standard deviations from the mode for a hot tower
0038   ,inner_f(inner_f)
0039   ,outer_f(outer_f)     // Fiducial cut for Sectors and IBs
0040 {
0041   std::cout << "towerid::towerid(const std::string &name) Calling ctor" << std::endl;
0042 }
0043 //__________________________________
0044 hcal_towerid::~hcal_towerid()
0045 {
0046   std::cout << "towerid::~towerid() Calling dtor" << std::endl;
0047 }
0048 //_____________________________
0049 int hcal_towerid::Init(PHCompositeNode *topNode)
0050 {
0051     //initialize output file with hot channels
0052   out = new TFile(Outfile.c_str(),"RECREATE");
0053   T = new TTree("T_inner","T_inner");
0054   T2 = new TTree("T_outer","T_outer");
0055   T -> Branch("hot_channels",&m_hot_channels_i);
0056   T2-> Branch("hot_channels",&m_hot_channels_o);
0057
0058   std::cout << "towerid::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0059   return Fun4AllReturnCodes::EVENT_OK;
0060
0061 }
0062 //_____________________________
0063 int hcal_towerid::InitRun(PHCompositeNode *topNode)
0064 {
0065   std::cout << "towerid::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0066   return Fun4AllReturnCodes::EVENT_OK;
0067 }
0068 //_____________________________
0069 int hcal_towerid::process_event(PHCompositeNode *topNode)
0070 {
0071
0072 //Get TowerInfoContainer
0073 TowerInfoContainer *hcalTowerContainer_i;
0074 TowerInfoContainer *hcalTowerContainer_o;
0075
0076 hcalTowerContainer_i =  findNode::getClass<TowerInfoContainer>(topNode,"TOWERS_HCALIN");
0077 hcalTowerContainer_o =  findNode::getClass<TowerInfoContainer>(topNode,"TOWERS_HCALOUT");
0078  if(!hcalTowerContainer_o)
0079     {
0080       std::cout << PHWHERE << "towerid::process_event Could not find node TOWERS_CEMC"  << std::endl;
0081       return Fun4AllReturnCodes::ABORTEVENT;
0082     }
0083 //iterate through all towers, incrementing their Frequency arrays if they record a hit
0084
0085  bool goodevent = false;
0086  int tower_range_i = hcalTowerContainer_i->size();
0087  int tower_range_o = hcalTowerContainer_o->size();
0088         for(int j = 0; j < tower_range_i; j++){
0089
0090         double energy = hcalTowerContainer_i -> get_tower_at_channel(j) -> get_energy();
0091         if(energy > adccut_i){
0092         itowerF[j]++;
0093         goodevent = true; //counter of events with nonzero EMCal energy
0094         }
0095
0096         //std::cout << energy << std::endl;
0097         itowerE[j]+=energy;
0098
0099     }
0100     for(int j = 0; j < tower_range_o; j++){
0101
0102             double energy = hcalTowerContainer_o -> get_tower_at_channel(j) -> get_energy();
0103                 if(energy > adccut_o){
0104                 otowerF[j]++;
0105                 goodevent = true; //counter of events with nonzero EMCal energy
0106                 }
0107
0108                 //std::cout << energy << std::endl;
0109                 otowerE[j]+=energy;
0110
0111         }
0112     if(goodevent == true){goodevents++;}
0113   return Fun4AllReturnCodes::EVENT_OK;
0114 }
0115 //__________________________
0116 int hcal_towerid::ResetEvent(PHCompositeNode *topNode)
0117 {
0118   return Fun4AllReturnCodes::EVENT_OK;
0119 }
0120
0121
0122 //__________________________
0123 int hcal_towerid::EndRun(const int runnumber)
0124 {
0125
0126  // initialize histograms: Initial frequency spectra, Energy spectrum, and Frequency spectrum with Fiducial cuts
0127
0128   Fspeci_i->SetBins(goodevents,0,goodevents);
0129   Fspeci_o->SetBins(goodevents,0,goodevents);
0130
0131   Fspec_i->SetBins(goodevents,0,inner_f*goodevents);
0132   Fspec_o->SetBins(goodevents,0,outer_f*goodevents);
0133
0134   Espec_i->SetBins(goodevents,0,goodevents);
0135   Espec_o->SetBins(goodevents,0,goodevents);
0136  //fill histograms
0137
0138  for(int i = 0; i < 1536; i++){
0139
0140     Fspec_i->Fill(itowerF[i]);
0141     Espec_i->Fill(itowerE[i]);
0142     Fspeci_i->Fill(itowerF[i]);
0143
0144     Fspec_o->Fill(otowerF[i]);
0145         Espec_o->Fill(otowerE[i]);
0146         Fspeci_o->Fill(otowerF[i]);
0147     }
0148
0149   //kill zeros:
0150
0151   Fspec_i->SetBinContent(1,0);
0152   Fspec_o->SetBinContent(1,0);
0153
0154     float cutoffFreq_i;
0155     float cutoffFreq_i_lo;
0156     float cutoffFreq_o;
0157     float cutoffFreq_o_lo;
0158
0159     cutoffFreq_i = Fspec_i->GetStdDev()*sigmas_hi + Fspec_i->GetXaxis()->GetBinCenter(Fspec_i->GetMaximumBin());
0160     cutoffFreq_i_lo = -1*Fspec_i->GetStdDev()*sigmas_lo + Fspec_i->GetXaxis()->GetBinCenter(Fspec_i->GetMaximumBin());
0161     cutoffFreq_o = Fspec_o->GetStdDev()*sigmas_hi + Fspec_o->GetXaxis()->GetBinCenter(Fspec_o->GetMaximumBin());
0162     cutoffFreq_o_lo = -1*Fspec_o->GetStdDev()*sigmas_lo + Fspec_o->GetXaxis()->GetBinCenter(Fspec_o->GetMaximumBin());
0163
0164     std::cout << "towerid::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0165     std::cout << "Inner hot tower cutoff: " << cutoffFreq_i << std::endl;
0166     std::cout << "Inner cold tower cutoff: " << cutoffFreq_i_lo << std::endl;
0167
0168     std::cout << "Outer hot tower cutoff: " << cutoffFreq_o << std::endl;
0169         std::cout << "OUter cold tower cutoff: " << cutoffFreq_o_lo << std::endl;
0170
0171   for(int i = 0; i < 1536; i++){
0172
0173     if(itowerF[i] > cutoffFreq_i){
0174         ihottowers[i]++;
0175     }
0176     if(otowerF[i]> cutoffFreq_o){
0177         ohottowers[i]++;
0178     }
0179     if(itowerF[i]==0){
0180         ideadtowers[i]++;
0181     }
0182     if(otowerF[i]==0){
0183         odeadtowers[i]++;
0184     }
0185     if(itowerF[i] < cutoffFreq_i_lo && itowerF[i]>0){
0186         icoldtowers[i]++;
0187     }
0188     if (otowerF[i] < cutoffFreq_o_lo && otowerF[i] > 0){
0189         ocoldtowers[i]++;
0190     }
0191     }
0192
0193   std::cout << "towerid::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0194   return Fun4AllReturnCodes::EVENT_OK;
0195 }
0196
0197 //____________________________________________________________________________..
0198 int hcal_towerid::End(PHCompositeNode *topNode)
0199 {
0200   std::cout << "towerid::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0201
0202    CDBTTree*icdbttree = new CDBTTree(cdbtreename_i);
0203    CDBTTree*ocdbttree = new CDBTTree(cdbtreename_o);
0204    std::string fieldname = "status";
0205    for(int i = 0; i<1536; i++){
0206
0207                 if(ihottowers[i] >= 0.5){
0208                         m_hot_channels_i = 2;
0209                         T->Fill();
0210             icdbttree->SetIntValue(i,fieldname,2);
0211                 }
0212            else if(ideadtowers[i] >= 0.5){
0213             m_hot_channels_i = 1;
0214                         T->Fill();
0215             icdbttree->SetIntValue(i,fieldname,1);
0216          }
0217         else if(icoldtowers[i] > 0.5){
0218             m_hot_channels_i = 3;
0219                         T->Fill();
0220             icdbttree->SetIntValue(i,fieldname,3);
0221         }
0222         else{
0223             m_hot_channels_i = 0;
0224                         T->Fill();
0225             icdbttree->SetIntValue(i,fieldname,0);
0226         }
0227         }
0228
0229   for(int j = 0; j<1536; j++){
0230
0231                 if(ohottowers[j] >= 0.5){
0232                         m_hot_channels_o = 2;
0233                         T2->Fill();
0234                         ocdbttree->SetIntValue(j,fieldname,2);
0235                 }
0236                else if(odeadtowers[j] >= 0.5){
0237                         m_hot_channels_o = 1;
0238                         T2->Fill();
0239                         ocdbttree->SetIntValue(j,fieldname,1);
0240                  }
0241                 else if(ocoldtowers[j] > 0.5){
0242                         m_hot_channels_i = 3;
0243                         T2->Fill();
0244                         ocdbttree->SetIntValue(j,fieldname,3);
0245                 }
0246                 else{
0247                         m_hot_channels_i = 0;
0248                         T2->Fill();
0249                         ocdbttree->SetIntValue(j,fieldname,0);
0250                 }
0251         }
0252
0253   out -> cd();
0254   Fspec_i->Write();
0255   Fspeci_i->Write();
0256   Espec_i->Write();
0257   Fspec_o->Write();
0258   Fspeci_o->Write();
0259   Espec_o->Write();
0260   T -> Write();
0261   T2->Write();
0262   out -> Close();
0263   delete out;
0264   icdbttree->Commit();
0265   //icdbttree->Print();
0266   icdbttree->WriteCDBTTree();
0267   delete icdbttree;
0268
0269   ocdbttree->Commit();
0270   //ocdbttree->Print();
0271   ocdbttree->WriteCDBTTree();
0272   delete ocdbttree;
0273
0274   return Fun4AllReturnCodes::EVENT_OK;
0275  }
0276
0277  //____________________________________________________________________________..
0278  int hcal_towerid::Reset(PHCompositeNode *topNode)
0279    {
0280  std::cout << "towerid::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0281  return Fun4AllReturnCodes::EVENT_OK;
0282     }
0283 //____________________________________________________
0284 void hcal_towerid::Print(const std::string &what) const
0285 {
0286   std::cout << "towerid::Print(const std::string &what) const Printing info for " << what << std::endl;
0287 }
0288
0289 //____________________________________________________
0290
0291
0292
0293
0294
0295
0296