analysis/HCalib/HCalib.C

0001 //HCAL Calibration
0002 //Abhisek Sen
0003 //
0004 #include "HCalib.h"
0005 #include <fun4all/SubsysReco.h>
0006 #include <fun4all/Fun4AllServer.h>
0007 #include <fun4all/PHTFileServer.h>
0008 #include <phool/PHCompositeNode.h>
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010 #include <phool/getClass.h>
0011 #include <fun4all/Fun4AllHistoManager.h>
0012 #include <phool/getClass.h>
0013
0014 #include <g4detectors/PHG4HcalDefs.h>
0015 #include <phparameter/PHParameters.h>
0016 #include <g4detectors/PHG4Cell.h>
0017 #include <g4detectors/PHG4CellContainer.h>
0018 #include <g4detectors/PHG4CellDefs.h>
0019
0020 #include <g4detectors/PHG4ScintillatorSlat.h>
0021 #include <g4detectors/PHG4ScintillatorSlatContainer.h>
0022 #include <g4detectors/PHG4ScintillatorSlatDefs.h>
0023
0024 #include <calobase/RawTowerContainer.h>
0025 #include <calobase/RawTowerGeomContainer.h>
0026 #include <calobase/RawTower.h>
0027 #include <TH1F.h>
0028
0029 using namespace std;
0030
0031 HCalib::HCalib():
0032 SubsysReco("HCal_Calibration")
0033 {
0034  towers=0;
0035  slats = 0;
0036  threshold=0.000; //GeV
0037  outfile=0;
0038  fill_towers = 1;
0039  fill_slats = 1;
0040  is_proto = true;
0041  return;
0042 }
0043
0044 int HCalib::Init(PHCompositeNode *topNode)
0045 {
0046  outfile = new TFile("out_calib.root","RECREATE");
0047  return Fun4AllReturnCodes::EVENT_OK;
0048 }
0049
0050 int HCalib::process_event(PHCompositeNode *topNode)
0051 {
0052  //cout << "Processing event:: HCalib " << endl;
0053  string detnames[] = {"HCALIN","HCALOUT"};
0054  for(int seg=0; seg<2; seg++)
0055  {
0056   GetNodes(topNode,detnames[seg]);
0057
0058   if(fill_towers and towers)
0059   {
0060    auto twr_range = towers->getTowers();
0061    for (auto twr_iter = twr_range.first; twr_iter != twr_range.second; ++twr_iter)
0062     {
0063       RawTower* tower = twr_iter->second;
0064       int etabin = tower->get_bineta();
0065       int phibin = tower->get_binphi();
0066       int key = genkey(seg, etabin, phibin);
0067       auto it = tower_hists.find(key);
0068       if(it==tower_hists.end())
0069       {
0070        string hname=detnames[seg]+"_Tower_Energy_ieta_"+to_string(etabin)+"_iphi_"+to_string(phibin);
0071        TH1F *h = new TH1F(hname.c_str(),hname.c_str(),500,0,50);
0072        tower_hists[key] = h;
0073        h->Fill( 1e3*tower->get_energy());
0074       }
0075       else {
0076        tower_hists[key]->Fill( 1e3*tower->get_energy());
0077      }
0078     }
0079
0080    }
0081
0082
0083   /*
0084   for(unsigned int ieta=1; ieta<23; ieta++)
0085   {
0086    for(unsigned int iphi=1; iphi<63; iphi++)
0087    {
0088    RawTower *twr = towers->getTower(ieta,iphi);
0089    if(!twr) continue;
0090
0091
0092    float twr_e = twr->get_energy();
0093
0094    //cout << "Energy " << twr->get_energy() << endl;
0095    if(twr_e<threshold) continue;
0096
0097    bool good_tower = true; //false;
0098    if(!good_tower)
0099    {
0100    // check for ieta-1, ieta+1
0101    RawTower *twr_p = towers->getTower(ieta-1,iphi);
0102    RawTower *twr_n = towers->getTower(ieta+1,iphi);
0103    if(twr_p && twr_n)
0104    {
0105     if(twr_p->get_energy()>threshold &&
0106       twr_n->get_energy()>threshold ) good_tower=true;
0107    }
0108
0109    //check for iphi-1, iphi+1
0110    twr_p = towers->getTower(ieta,iphi-1);
0111    twr_n = towers->getTower(ieta,iphi-1);
0112    if(twr_p && twr_n)
0113    {
0114     if(twr_p->get_energy()>threshold &&
0115       twr_n->get_energy()>threshold ) good_tower=true;
0116    }
0117
0118    //Diagonal
0119    twr_p = towers->getTower(ieta-1,iphi-1);
0120    twr_n = towers->getTower(ieta+1,iphi+1);
0121    if(twr_p && twr_n)
0122    {
0123     if(twr_p->get_energy()>threshold &&
0124        twr_n->get_energy()>threshold ) good_tower=true;
0125    }
0126
0127    //Diagonal
0128    twr_p = towers->getTower(ieta-1,iphi+1);
0129    twr_n = towers->getTower(ieta+1,iphi-1);
0130    if(twr_p && twr_n)
0131    {
0132     if(twr_p->get_energy()>threshold &&
0133        twr_n->get_energy()>threshold ) good_tower=true;
0134    }
0135    }
0136
0137    if(good_tower) Fill(twr,seg);
0138
0139    }
0140    }
0141   }
0142    */
0143  if(fill_slats)
0144  {
0145    if(is_proto && proto_slats)
0146    {
0147     auto cell_range = proto_slats->getScintillatorSlats();
0148     for (auto cell_iter = cell_range.first; cell_iter != cell_range.second; ++cell_iter)
0149     {
0150       PHG4ScintillatorSlat *cell = cell_iter->second;
0151       int etabin = cell->get_column();
0152       int phibin = cell->get_row();
0153       int key = genkey(seg, etabin, phibin);
0154       auto it = slat_hists.find(key);
0155       if(it==slat_hists.end())
0156       {
0157        string hname=detnames[seg]+"_Slat_Energy_ieta_"+to_string(etabin)+"_iphi_"+to_string(phibin);
0158        TH1F *h = new TH1F(hname.c_str(),hname.c_str(),500,0,50);
0159        slat_hists[key] = h;
0160        h->Fill( 1e3*cell->get_light_yield() );
0161       }
0162       else {
0163        slat_hists[key]->Fill( 1e3*cell->get_light_yield() );
0164      }
0165     }
0166    }
0167
0168    if(slats)
0169    {
0170     auto cell_range = slats->getCells();
0171     for (auto cell_iter = cell_range.first; cell_iter != cell_range.second; ++cell_iter)
0172     {
0173       PHG4Cell *cell = cell_iter->second;
0174       int etabin = PHG4CellDefs::ScintillatorSlatBinning::get_column(cell->get_cellid());
0175       int phibin = PHG4CellDefs::ScintillatorSlatBinning::get_row(cell->get_cellid());
0176       int key = genkey(seg, etabin, phibin);
0177       auto it = slat_hists.find(key);
0178       if(it==slat_hists.end())
0179       {
0180        string hname=detnames[seg]+"_Slat_Energy_ieta_"+to_string(etabin)+"_iphi_"+to_string(phibin);
0181        TH1F *h = new TH1F(hname.c_str(),hname.c_str(),500,0,50);
0182        slat_hists[key] = h;
0183        h->Fill( 1e3*cell->get_light_yield() );
0184       }
0185       else {
0186        slat_hists[key]->Fill( 1e3*cell->get_light_yield() );
0187      }
0188     }
0189    }
0190  }
0191
0192  }
0193
0194  return Fun4AllReturnCodes::EVENT_OK;
0195 }
0196
0197 int HCalib::End(PHCompositeNode *topNode)
0198 {
0199  outfile->cd();
0200  for( auto it: tower_hists)
0201  {
0202   if(it.second->GetEntries()>0)
0203    it.second->Write();
0204  }
0205
0206  for(auto it : slat_hists)
0207  {
0208   if(it.second->GetEntries()>0)
0209    it.second->Write();
0210  }
0211
0212  outfile->Close();
0213  return Fun4AllReturnCodes::EVENT_OK;
0214 }
0215
0216 void HCalib::GetNodes(PHCompositeNode *topNode, const string &det)
0217 {
0218     static int ncalls = 0;
0219     //string towernode = "TOWER_CALIB_" + det;
0220     string towernode = "TOWER_SIM_" + det;
0221     towers = findNode::getClass<RawTowerContainer>(topNode,towernode.c_str());
0222   if (!towers) {
0223       if(ncalls<5) cerr << PHWHERE << " ERROR: Can't find " << towernode << endl;
0224       ncalls++;
0225       //exit(-1);
0226     }
0227
0228   //Slats
0229   std::string cellnodename = "G4CELL_" + det;
0230   slats = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
0231   if(!slats)
0232   {
0233     if(ncalls<5) cerr << PHWHERE << " Node missing " << cellnodename << endl;
0234   }
0235
0236   //Proto slats
0237   if( is_proto )
0238   {
0239    proto_slats = findNode::getClass<PHG4ScintillatorSlatContainer>(topNode, cellnodename);
0240    if(!proto_slats)
0241   {
0242     if(ncalls<5) cerr << PHWHERE << " Prototype Node missing " << cellnodename << endl;
0243   }
0244   }
0245
0246 }
0247
0248 int HCalib::genkey(const int detid, const int etabin, const int phibin)
0249 {
0250  //Assuming
0251  //At most detid : 2, etabin: 24, phibin:256
0252  return 30000*detid + 500*etabin + phibin;
0253 }