Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:14:14 

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 #include <TStyle.h>
0020 #include <TSystem.h>
0021 #include <TPad.h>
0022 #include <TLine.h>
0023 
0024 //Tower stuff
0025 #include <calobase/TowerInfoContainer.h>
0026 #include <calobase/TowerInfo.h>
0027 #include <calobase/TowerInfoDefs.h>
0028 
0029 #include <cdbobjects/CDBTTree.h>
0030 
0031 //________________________________
0032 hcal_towerid::hcal_towerid(const std::string &name, const std::string &outfilename, const std::string &icdbtreename, const std::string &ocdbtreename, float adccut_i, float adccut_o, float sigmas_lo, float sigmas_hi, float inner_f,float outer_f): 
0033 SubsysReco(name)
0034 ,Tinner(nullptr)
0035 ,Touter(nullptr)
0036 ,Outfile(outfilename)   
0037 ,icdbtreename(icdbtreename)
0038   ,ocdbtreename(ocdbtreename)
0039   ,adccut_i(adccut_i)
0040 ,adccut_o(adccut_o)     // Minimum ADC for Kurary fibers to register a hit
0041   ,sigmas_lo(sigmas_lo)     // # of standard deviations from the mode for a cold tower
0042   ,sigmas_hi(sigmas_hi)     // # of standard deviations from the mode for a hot tower
0043   ,inner_f(inner_f)
0044   ,outer_f(outer_f)     // Fiducial cut for Sectors and IBs 
0045 {
0046   std::cout << "towerid::towerid(const std::string &name) Calling ctor" << std::endl;
0047 }
0048 //__________________________________
0049 hcal_towerid::~hcal_towerid()
0050 {
0051   std::cout << "towerid::~towerid() Calling dtor" << std::endl;
0052 }
0053 //_____________________________
0054 int hcal_towerid::Init(PHCompositeNode *topNode)
0055 {
0056     //initialize output file with hot channels
0057   std::cout << "towerid::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0058   return Fun4AllReturnCodes::EVENT_OK;
0059 
0060 }
0061 //_____________________________
0062 int hcal_towerid::InitRun(PHCompositeNode *topNode)
0063 {
0064   std::cout << "towerid::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0065   return Fun4AllReturnCodes::EVENT_OK;
0066 }
0067 //_____________________________
0068 int hcal_towerid::process_event(PHCompositeNode *topNode)
0069 {
0070 
0071 //Get TowerInfoContainer
0072   TowerInfoContainer *hcalTowerContainer_i;
0073 //TowerInfoContainer *hcalTowerContainer_o;
0074 
0075   hcalTowerContainer_i =  findNode::getClass<TowerInfoContainer>(topNode,"TOWERS_HCALIN");
0076 // adding HCALIN
0077   if(!hcalTowerContainer_i)
0078     {
0079       std::cout << PHWHERE << "towerid::process_event Could not find node TOWERS_HCALIN"  << std::endl;
0080       return Fun4AllReturnCodes::ABORTEVENT;
0081     }
0082   
0083   TowerInfoContainer *hcalTowerContainer_o;
0084   hcalTowerContainer_o =  findNode::getClass<TowerInfoContainer>(topNode,"TOWERS_HCALOUT");
0085   if(!hcalTowerContainer_o)
0086     {
0087       std::cout << PHWHERE << "towerid::process_event Could not find node TOWERS_HCALOUT"  << std::endl;
0088       return Fun4AllReturnCodes::ABORTEVENT;
0089     }
0090   //iterate through all towers, incrementing their Frequency arrays if they record a hit
0091 
0092   bool goodevent = false;
0093   int tower_range_i = hcalTowerContainer_i->size();
0094   int tower_range_o = hcalTowerContainer_o->size();
0095   iHCal_hist = new TH1F("iHCal_energy", "Inner HCal Energy Distribution", 100, 0, 100);
0096   oHCal_hist = new TH1F("oHCal_energy", "Outer HCal Energy Distribution", 100, 0, 100);
0097   // Return EVENT_OK if at least one tower had non-zero energy
0098    for(int j = 0; j < tower_range_i; j++){
0099 
0100     double energy = hcalTowerContainer_i -> get_tower_at_channel(j) -> get_energy();
0101     iHCal_hist->Fill(energy);
0102     if(energy > adccut_i){
0103       itowerF[j]++;
0104       goodevent = true; //counter of events with nonzero IHCal energy
0105     }               
0106     std::cout << energy << std::endl;
0107     itowerE[j]+=energy;
0108    }
0109   
0110   //    int tower_range_o = hcalTowerContainer_o->size();
0111   for(int j = 0; j < tower_range_o; j++){
0112     
0113     double energy = hcalTowerContainer_o -> get_tower_at_channel(j) -> get_energy();
0114     oHCal_hist->Fill(energy);
0115     if(energy > adccut_o){
0116       otowerF[j]++;
0117       goodevent = true; //counter of events with nonzero OHCal energy
0118     }
0119     std::cout << energy << std::endl;
0120     otowerE[j]+=energy;
0121   }
0122   if(goodevent == true){goodevents++;}
0123   return Fun4AllReturnCodes::EVENT_OK;
0124 }
0125 //__________________________
0126 int hcal_towerid::ResetEvent(PHCompositeNode *topNode)
0127 {
0128   return Fun4AllReturnCodes::EVENT_OK;
0129 }
0130 
0131 //__________________________ 
0132 void hcal_towerid::FillHistograms(const int runnumber, const int segment)
0133 {
0134   const int nBins = 101;
0135   
0136   out = new TFile(Outfile.c_str(),"RECREATE");
0137  
0138   Fspec_i = new TH2F("Fspec_i","Fspec_i",nTowers+1,-0.5,nTowers+0.5,nBins,0.-inner_f/(2*nBins),inner_f+inner_f/(2*nBins));
0139   Fspec_o = new TH2F("Fspec_o","Fspec_o",nTowers+1,-0.5,nTowers+0.5,nBins,0.-outer_f/(2*nBins),outer_f+outer_f/(2*nBins));
0140   
0141   Fspeci_i = new TH2F("Fspeci_i","Fspeci_i",nTowers+1,-0.5,nTowers+0.5,nBins,0.-1/(2*nBins),1);
0142   Fspeci_o = new TH2F("Fspeci_o","Fspeci_o",nTowers+1,-0.5,nTowers+0.5,nBins,0.-1/(2*nBins),1);
0143   
0144   Espec_i = new TH2F("Espec_i","Espec_i",nTowers+1,-0.5,nTowers-0.5,nBins,0.-1/(2*nBins),1);
0145   Espec_o = new TH2F("Espec_o","Espec_o",nTowers+1,-0.5,nTowers-0.5,nBins,0.-1/(2*nBins),1);
0146   for(int i = 0; i < nTowers; i++){
0147     
0148     Fspec_i->Fill(i,itowerF[i]/goodevents);
0149     Espec_i->Fill(i,itowerE[i]/goodevents);
0150     Fspeci_i->Fill(i,itowerF[i]/goodevents);
0151     
0152     Fspec_o->Fill(i,otowerF[i]/goodevents);
0153     Espec_o->Fill(i,otowerE[i]/goodevents);
0154     Fspeci_o->Fill(i,otowerF[i]/goodevents);
0155   }
0156 
0157   out ->cd();
0158   
0159   Fspec_i->Write();
0160   delete Fspec_i;
0161   
0162   Fspeci_i->Write();
0163   delete Fspeci_i;
0164   
0165   Espec_i->Write();
0166   delete Espec_i;
0167 
0168   Fspec_o->Write();
0169   delete Fspec_o;
0170 
0171   Fspeci_o->Write();
0172   delete Fspeci_o;
0173   
0174   Espec_o->Write();
0175   delete Espec_o;
0176 
0177   iHCal_hist->Write();
0178   delete iHCal_hist;
0179 
0180   oHCal_hist->Write();
0181   delete oHCal_hist;
0182    
0183 }
0184 //__________________________
0185 
0186 void hcal_towerid::CalculateCutOffs(const int runnumber)
0187 {
0188   TFile *histsIn = new TFile(Form("output/DST_CALOR-%08d_badTowerMapTree.root",runnumber));
0189   //  TFile *histsIn = new TFile(Form("DST_CALOR_badTowerMapTree.root"));
0190   
0191   float cutoffFreq_i;
0192   float cutoffFreq_i_lo;
0193  
0194   float cutoffFreq_o;
0195   float cutoffFreq_o_lo;
0196   
0197   TH1F *dummyProj = NULL;
0198   gStyle -> SetOptStat(0);
0199   Fspec_i = (TH2F*)histsIn -> Get("Fspec_i");
0200   dummyProj = (TH1F*)Fspec_i -> ProjectionY("dummy");
0201   dummyProj -> SetBinContent(1,0);
0202   cutoffFreq_i = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0203   cutoffFreq_i_lo = -1*dummyProj->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0204   
0205   Fspec_o = (TH2F*)histsIn -> Get("Fspec_o");
0206   dummyProj = (TH1F*)Fspec_o -> ProjectionY("dummy");
0207   dummyProj -> SetBinContent(1,0);
0208   cutoffFreq_o = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0209   cutoffFreq_o_lo = -1*Fspec_o->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0210   
0211   Fspeci_i = (TH2F*)histsIn -> Get("Fspeci_i");
0212   Fspeci_i -> SetTitle("Tower ID; Hit Frequency");
0213   
0214   Fspeci_o = (TH2F*)histsIn -> Get("Fspeci_o");
0215   Fspeci_o -> SetTitle("Tower ID; Hit Frequency");
0216   
0217   for(int i = 0; i < nTowers; i++){ 
0218 
0219     TH1F *fSpecProj_i = (TH1F*)Fspeci_i->ProjectionY("dummyi",i+1,i+1);
0220   //  TH1F *fSpecProj_o = (TH1F*)Fspeci_o->ProjectionY("dummyo",i+1,i+1);
0221     
0222     if(fSpecProj_i->GetBinCenter(fSpecProj_i->GetMaximumBin()) > cutoffFreq_i){
0223       ihottowers[i]++;
0224     }
0225     else if(fSpecProj_i -> GetBinCenter(fSpecProj_i->GetMaximumBin()) == 0){
0226       ideadtowers[i]++;
0227   }
0228     else if (fSpecProj_i->GetBinCenter(fSpecProj_i->GetMaximumBin()) < cutoffFreq_i_lo &&  fSpecProj_i->GetBinCenter(fSpecProj_i->GetMaximumBin()) >0 ){
0229       icoldtowers[i]++;
0230     }
0231   }
0232 
0233   for(int i = 0; i < nTowers; i++){
0234     TH1F *fSpecProj_o = (TH1F*)Fspeci_o->ProjectionY("dummyo",i+1,i+1); 
0235     if(fSpecProj_o->GetBinCenter(fSpecProj_o->GetMaximumBin()) > cutoffFreq_o){
0236       ohottowers[i]++;
0237   }
0238     else if(fSpecProj_o -> GetBinCenter(fSpecProj_o->GetMaximumBin()) == 0){
0239       odeadtowers[i]++;
0240     }
0241     else if (fSpecProj_o->GetBinCenter(fSpecProj_o->GetMaximumBin()) < cutoffFreq_o_lo &&  fSpecProj_o->GetBinCenter(fSpecProj_o->GetMaximumBin()) >0 ){
0242       ocoldtowers[i]++;
0243     }
0244   }
0245   //adding
0246   hot_regions = 0;
0247   cold_regions = 0;
0248   std::cout << "removing hot/cold regions" << std::endl;
0249   while(hot_regions == 1 || cold_regions ==1){
0250     std::cout << "hot/cold tower detected. Running another pass for hot towers" << std::endl;
0251     dummyProj = (TH1F*)Fspec_i -> ProjectionY("dummy");
0252     dummyProj -> SetBinContent(1,0);
0253     cutoffFreq_i = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0254     cutoffFreq_i_lo = -1*dummyProj->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0255     
0256     dummyProj = (TH1F*)Fspec_o -> ProjectionY("dummy");
0257     dummyProj -> SetBinContent(1,0);
0258     cutoffFreq_o = dummyProj->GetStdDev()*sigmas_hi + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0259     cutoffFreq_o_lo = -1*Fspec_o->GetStdDev()*sigmas_lo + dummyProj->GetBinCenter(dummyProj->GetMaximumBin());
0260     for(int i = 0; i < nTowers; i++){
0261       
0262       TH1F *fSpecProj_i = (TH1F*)Fspeci_i->ProjectionY("dummyi",i+1,i+1);
0263     //    TH1F *fSpecProj_o = (TH1F*)Fspeci_o->ProjectionY("dummyo",i+1,i+1);
0264 
0265       if(fSpecProj_i->GetBinCenter(fSpecProj_i->GetMaximumBin()) > cutoffFreq_i){
0266     ihottowers[i]++;
0267       }
0268       else if(fSpecProj_i -> GetBinCenter(fSpecProj_i->GetMaximumBin()) == 0){
0269     ideadtowers[i]++;
0270       }
0271       else if (fSpecProj_i->GetBinCenter(fSpecProj_i->GetMaximumBin()) < cutoffFreq_i_lo &&  fSpecProj_i->GetBinCenter(fSpecProj_i->GetMaximumBin()) >0 ){
0272     icoldtowers[i]++;
0273       }
0274     }
0275     for(int i = 0; i < nTowers; i++){
0276       TH1F *fSpecProj_o = (TH1F*)Fspeci_o->ProjectionY("dummyo",i+1,i+1);
0277       if(fSpecProj_o->GetBinCenter(fSpecProj_o->GetMaximumBin()) > cutoffFreq_o){
0278     ohottowers[i]++;
0279       }
0280       else if(fSpecProj_o -> GetBinCenter(fSpecProj_o->GetMaximumBin()) == 0){
0281     odeadtowers[i]++;
0282       }
0283       else if (fSpecProj_o->GetBinCenter(fSpecProj_o->GetMaximumBin()) < cutoffFreq_o_lo &&  fSpecProj_o->GetBinCenter(fSpecProj_o->GetMaximumBin()) >0 ){
0284     ocoldtowers[i]++;
0285       }
0286     }
0287     
0288     hot_regions = 0;
0289     cold_regions = 0;
0290   } 
0291     std::cout << "towerid::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0292     std::cout << "Inner Cutoff Frequency: " << cutoffFreq_i << std::endl;
0293     std::cout << "Outer Cutoff Frequency: " << cutoffFreq_o << std::endl;
0294 }
0295 //__________________________
0296   void hcal_towerid::WriteCDBTree(const int runnumber)
0297   {
0298     //    TFile *histsIn = new TFile(Form("output/%d/DST_CALOR-%08d_badTowerMapTree.root",runnumber,runnumber));
0299     TFile *cdbOut = new TFile(Form("cdbMaps/DST_CALOR-%08d_badTowerMapCDBTTree.root",runnumber),"RECREATE");
0300     //  TFile *cdbOut = new TFile(Form("badTowerMapCDBTTree.root"),"RECREATE");
0301     //TFile *cdbOut = new TFile(Form("/sphenix/u/shumail0001/analysis/JS-Jet/HCalHotTowerFinder/macro/HCAL-%08d_badTowerMapCDBTTree.root",runnumber),"RECREATE");
0302     
0303     Tinner = new TTree("T_inner","T_inner");
0304     Tinner -> Branch("hot_channels_i",&hot_channels_i);
0305     Touter = new TTree("T_outer","T_outer");
0306     Touter -> Branch("hot_channels_o",&hot_channels_o);
0307 
0308     CDBTTree*icdbttree = new CDBTTree(icdbtreename);
0309     CDBTTree*ocdbttree = new CDBTTree(ocdbtreename);
0310     std::string fieldname = "status";
0311     for(int i = 0; i<nTowers; i++){
0312 
0313     unsigned int key = TowerInfoDefs::encode_hcal(i);
0314     
0315       if(ihottowers[i] >= 0.5){
0316     hot_channels_i = 2;
0317     Tinner->Fill();
0318         icdbttree->SetIntValue(key,fieldname,2);
0319       }
0320       else if(ideadtowers[i] >= 0.5){
0321     hot_channels_i = 1;
0322     Tinner->Fill();
0323     icdbttree->SetIntValue(key,fieldname,1);
0324       }
0325       else if(icoldtowers[i] > 0.5){
0326     hot_channels_i = 3;
0327     Tinner->Fill();
0328         icdbttree->SetIntValue(key,fieldname,3);
0329       }
0330       else{
0331     hot_channels_i = 0;
0332     Tinner->Fill();
0333         icdbttree->SetIntValue(key,fieldname,0);
0334       }
0335     }
0336       for(int i = 0; i<nTowers; i++){
0337 
0338     unsigned int key1 = TowerInfoDefs::encode_hcal(i);
0339 
0340     if(ohottowers[i] >= 0.5){
0341       hot_channels_o = 2;
0342       Touter->Fill();
0343       ocdbttree->SetIntValue(key1,fieldname,2);
0344     }
0345     else if(odeadtowers[i] >= 0.5){
0346       hot_channels_o = 1;
0347       Touter->Fill();
0348       ocdbttree->SetIntValue(key1,fieldname,1);
0349     }
0350     else if(ocoldtowers[i] > 0.5){
0351       hot_channels_o = 3;
0352       Touter->Fill();
0353           ocdbttree->SetIntValue(key1,fieldname,3);
0354     }
0355     else{
0356       hot_channels_o = 0;
0357       Touter->Fill();
0358       ocdbttree->SetIntValue(key1,fieldname,0);
0359     }
0360       }
0361 
0362 
0363   
0364     cdbOut->cd();
0365     Tinner -> Write();
0366     icdbttree->Commit();
0367     //icdbttree->Print();
0368     icdbttree->WriteCDBTTree();
0369     delete icdbttree;
0370     //cdbOut->Close();
0371 
0372     Touter -> Write();
0373     ocdbttree->Commit();
0374     //ocdbttree->Print()
0375     ocdbttree->WriteCDBTTree();
0376     delete ocdbttree;
0377     cdbOut->Close();
0378 
0379   }
0380 
0381  //____________________________________________________________________________..
0382  int hcal_towerid::Reset(PHCompositeNode *topNode)
0383    {
0384  std::cout << "towerid::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0385  return Fun4AllReturnCodes::EVENT_OK;
0386     }
0387 //____________________________________________________
0388 void hcal_towerid::Print(const std::string &what) const
0389 {
0390   std::cout << "towerid::Print(const std::string &what) const Printing info for " << what << std::endl;
0391 }
0392 
0393 //____________________________________________________
0394  int hcal_towerid::End(PHCompositeNode *topNode)
0395  {
0396    std::cout << "towerid::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0397    return Fun4AllReturnCodes::EVENT_OK;
0398  }
0399  //____________________________________________________
0400  int hcal_towerid::EndRun(int runnumber)
0401  {
0402    std::cout << "towerid::EndRun: this is the end of run: "<< runnumber << std::endl;
0403    return Fun4AllReturnCodes::EVENT_OK;
0404  }
0405  //____________________________________________________
0406
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!