HCal/TowerCalib/Proto4TowerCalib.C

0001 #include "Proto4TowerCalib.h"
0002
0003 #include <calobase/RawTowerContainer.h>
0004 #include <ffaobjects/EventHeader.h>
0005 #include <pdbcalbase/PdbParameterMap.h>
0006 #include <phparameter/PHParameters.h>
0007 #include <prototype4/RawTower_Prototype4.h>
0008 #include <prototype4/RawTower_Temperature.h>
0009
0010 #include <Event/EventTypes.h>
0011 #include <fun4all/Fun4AllHistoManager.h>
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <fun4all/Fun4AllServer.h>
0014 #include <fun4all/PHTFileServer.h>
0015 #include <fun4all/SubsysReco.h>
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/getClass.h>
0018
0019 #include <phool/PHCompositeNode.h>
0020
0021 #include <g4main/PHG4Particle.h>
0022 #include <g4main/PHG4TruthInfoContainer.h>
0023 #include <g4main/PHG4VtxPoint.h>
0024
0025 #include <TFile.h>
0026 #include <TString.h>
0027 #include <TH1F.h>
0028 #include <TH2F.h>
0029 #include <TLorentzVector.h>
0030 #include <TNtuple.h>
0031 #include <TVector3.h>
0032 #include <TChain.h>
0033
0034 #include <algorithm>
0035 #include <cassert>
0036 #include <cmath>
0037 #include <exception>
0038 #include <iostream>
0039 #include <set>
0040 #include <sstream>
0041 #include <stdexcept>
0042 #include <vector>
0043
0044 using namespace std;
0045
0046 ClassImp(Proto4TowerCalib::HCAL_Tower);
0047
0048 Proto4TowerCalib::Proto4TowerCalib(const std::string &filename)
0049   : SubsysReco("Proto4TowerCalib")
0050   , _is_sim(true)
0051   , _filename(filename)
0052   , _ievent(0)
0053 {
0054   Verbosity(1);
0055
0056   _tower.reset();
0057
0058   _mInPut_flag = 1;
0059   _mStartEvent = -1;
0060   _mStopEvent = -1;
0061   _mDet = "HCALIN";
0062   _mCol = 0;
0063   _mPedestal = 100.0;
0064   reset_pedestal();
0065 }
0066
0067 Proto4TowerCalib::~Proto4TowerCalib()
0068 {
0069 }
0070
0071 Fun4AllHistoManager *
0072 Proto4TowerCalib::get_HistoManager()
0073 {
0074   Fun4AllServer *se = Fun4AllServer::instance();
0075   Fun4AllHistoManager *hm = se->getHistoManager("Proto4TowerCalib_HISTOS");
0076
0077   if (not hm)
0078   {
0079     cout
0080         << "Proto4TowerCalib::get_HistoManager - Making Fun4AllHistoManager Proto4TowerCalib_HISTOS"
0081         << endl;
0082     hm = new Fun4AllHistoManager("Proto4TowerCalib_HISTOS");
0083     se->registerHistoManager(hm);
0084   }
0085
0086   assert(hm);
0087
0088   return hm;
0089 }
0090
0091 int Proto4TowerCalib::InitRun(PHCompositeNode *topNode)
0092 {
0093   if (Verbosity())
0094     cout << "Proto4TowerCalib::InitRun" << endl;
0095
0096   _ievent = 0;
0097
0098   PHNodeIterator iter(topNode);
0099   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst(
0100       "PHCompositeNode", "DST"));
0101   if (!dstNode)
0102   {
0103     std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0104     throw runtime_error(
0105         "Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
0106   }
0107
0108   return Fun4AllReturnCodes::EVENT_OK;
0109 }
0110
0111 int Proto4TowerCalib::End(PHCompositeNode *topNode)
0112 {
0113   cout << "Proto4TowerCalib::End - write to " << _filename << endl;
0114   PHTFileServer::get().cd(_filename);
0115
0116   Fun4AllHistoManager *hm = get_HistoManager();
0117   assert(hm);
0118   for (unsigned int i = 0; i < hm->nHistos(); i++)
0119     hm->getHisto(i)->Write();
0120
0121   //  if (_T_EMCalTrk)
0122   //    _T_EMCalTrk->Write();
0123
0124   return Fun4AllReturnCodes::EVENT_OK;
0125 }
0126
0127 int Proto4TowerCalib::Init(PHCompositeNode *topNode)
0128 {
0129   _ievent = 0;
0130
0131   cout << "Proto4TowerCalib::get_HistoManager - Making PHTFileServer "
0132        << _filename << endl;
0133   PHTFileServer::get().open(_filename, "RECREATE");
0134
0135   Fun4AllHistoManager *hm = get_HistoManager();
0136   assert(hm);
0137
0138   //! Envent nomalization
0139   TH1F *hNormalization = new TH1F("hNormalization", "hNormalization", 1, 0.5, 1.5);
0140   hNormalization->GetXaxis()->SetBinLabel(1, "ALL");
0141   hNormalization->GetXaxis()->SetTitle("Cuts");
0142   hNormalization->GetYaxis()->SetTitle("Event count");
0143   hm->registerHisto(hNormalization);
0144
0145   // SIM HCALIN/HCALOUT
0146   if(_is_sim)
0147   {
0148     TH1F *h_hcalin_tower_sim[16];
0149     TH1F *h_hcalout_tower_sim[16];
0150     for(int i_tower = 0; i_tower < 16; ++i_tower)
0151     {
0152       string HistName_sim;
0153       HistName_sim = Form("h_hcalin_tower_%d_sim",i_tower);
0154       h_hcalin_tower_sim[i_tower] = new TH1F(HistName_sim.c_str(),HistName_sim.c_str(),500,-0.5,99.5);
0155       hm->registerHisto(h_hcalin_tower_sim[i_tower]);
0156
0157       HistName_sim = Form("h_hcalout_tower_%d_sim",i_tower);
0158       h_hcalout_tower_sim[i_tower] = new TH1F(HistName_sim.c_str(),HistName_sim.c_str(),500,-0.5,99.5);
0159       hm->registerHisto(h_hcalout_tower_sim[i_tower]);
0160     }
0161   }
0162
0163   // HCALIN LG
0164   TH1F *h_hcalin_lg_tower_raw[16];
0165   TH1F *h_hcalin_lg_tower_calib[16];
0166   for(int i_tower = 0; i_tower < 16; ++i_tower)
0167   {
0168     string HistName_raw = Form("h_hcalin_lg_tower_%d_raw",i_tower);
0169     h_hcalin_lg_tower_raw[i_tower] = new TH1F(HistName_raw.c_str(),HistName_raw.c_str(),40,0.5,16000.5);
0170     hm->registerHisto(h_hcalin_lg_tower_raw[i_tower]);
0171
0172     string HistName_calib = Form("h_hcalin_lg_tower_%d_calib",i_tower);
0173     h_hcalin_lg_tower_calib[i_tower] = new TH1F(HistName_calib.c_str(),HistName_calib.c_str(),100,-0.5,19.5);
0174     hm->registerHisto(h_hcalin_lg_tower_calib[i_tower]);
0175   }
0176
0177   // HCALOUT LG
0178   TH1F *h_hcalout_lg_tower_raw[16];
0179   TH1F *h_hcalout_lg_tower_calib[16];
0180   for(int i_tower = 0; i_tower < 16; ++i_tower)
0181   {
0182     string HistName_raw = Form("h_hcalout_lg_tower_%d_raw",i_tower);
0183     h_hcalout_lg_tower_raw[i_tower] = new TH1F(HistName_raw.c_str(),HistName_raw.c_str(),40,0.5,16000.5);
0184     hm->registerHisto(h_hcalout_lg_tower_raw[i_tower]);
0185
0186     string HistName_calib = Form("h_hcalout_lg_tower_%d_calib",i_tower);
0187     h_hcalout_lg_tower_calib[i_tower] = new TH1F(HistName_calib.c_str(),HistName_calib.c_str(),100,-0.5,19.5);
0188     hm->registerHisto(h_hcalout_lg_tower_calib[i_tower]);
0189   }
0190
0191   // HCALOUT HG
0192   TH1F *h_hcalout_hg_tower_raw[16];
0193   TH1F *h_hcalout_hg_tower_calib[16];
0194   for(int i_tower = 0; i_tower < 16; ++i_tower)
0195   {
0196     string HistName_raw = Form("h_hcalout_hg_tower_%d_raw",i_tower);
0197     h_hcalout_hg_tower_raw[i_tower] = new TH1F(HistName_raw.c_str(),HistName_raw.c_str(),40,0.5,16000.5);
0198     hm->registerHisto(h_hcalout_hg_tower_raw[i_tower]);
0199
0200     string HistName_calib = Form("h_hcalout_hg_tower_%d_calib",i_tower);
0201     h_hcalout_hg_tower_calib[i_tower] = new TH1F(HistName_calib.c_str(),HistName_calib.c_str(),100,-0.5,19.5);
0202     hm->registerHisto(h_hcalout_hg_tower_calib[i_tower]);
0203   }
0204
0205   // help index files with TChain
0206   TTree *T = new TTree("HCAL_Info", "HCAL_Info");
0207   assert(T);
0208   hm->registerHisto(T);
0209
0210   T->Branch("TowerCalib", &_tower);
0211
0212   return Fun4AllReturnCodes::EVENT_OK;
0213 }
0214
0215 int Proto4TowerCalib::process_event(PHCompositeNode *topNode)
0216 {
0217   if (Verbosity() > 2)
0218     cout << "Proto4TowerCalib::process_event() entered" << endl;
0219
0220   // init eval objects
0221   _tower.reset();
0222
0223   Fun4AllHistoManager *hm = get_HistoManager();
0224   assert(hm);
0225
0226   // EventHeader *eventheader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0227
0228   /*
0229   if (eventheader->get_EvtType() != DATAEVENT)
0230   {
0231     cout << __PRETTY_FUNCTION__
0232          << " disgard this event: " << endl;
0233
0234     eventheader->identify();
0235
0236     return Fun4AllReturnCodes::DISCARDEVENT;
0237   }
0238   */
0239
0240   if (_is_sim)
0241   {
0242     PHG4TruthInfoContainer *truthInfoList = findNode::getClass<
0243         PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0244
0245     assert(truthInfoList);
0246
0247
0248     const PHG4Particle *p =
0249         truthInfoList->GetPrimaryParticleRange().first->second;
0250     assert(p);
0251
0252     const PHG4VtxPoint *v = truthInfoList->GetVtx(p->get_vtx_id());
0253     assert(v);
0254   }
0255
0256   // normalization
0257   TH1F *hNormalization = dynamic_cast<TH1F *>(hm->getHisto("hNormalization"));
0258   assert(hNormalization);
0259
0260   hNormalization->Fill("ALL", 1);
0261
0262   // get DST nodes
0263
0264   // HCALIN/HCALOUT SIM information
0265   RawTowerContainer *TOWER_SIM_HCALIN = findNode::getClass<
0266       RawTowerContainer>(topNode, "TOWER_SIM_HCALIN");
0267   if(_is_sim)
0268   {
0269     assert(TOWER_SIM_HCALIN);
0270   }
0271
0272   RawTowerContainer *TOWER_SIM_HCALOUT = findNode::getClass<
0273       RawTowerContainer>(topNode, "TOWER_SIM_HCALOUT");
0274   if(_is_sim)
0275   {
0276     assert(TOWER_SIM_HCALOUT);
0277   }
0278
0279   // HCALIN LG information
0280   RawTowerContainer *TOWER_RAW_LG_HCALIN = findNode::getClass<
0281       RawTowerContainer>(topNode, "TOWER_RAW_LG_HCALIN");
0282   assert(TOWER_RAW_LG_HCALIN);
0283
0284   RawTowerContainer *TOWER_CALIB_LG_HCALIN = findNode::getClass<
0285       RawTowerContainer>(topNode, "TOWER_CALIB_LG_HCALIN");
0286   assert(TOWER_CALIB_LG_HCALIN);
0287
0288   // HCALOUT LG information
0289   RawTowerContainer *TOWER_RAW_LG_HCALOUT = findNode::getClass<
0290       RawTowerContainer>(topNode, "TOWER_RAW_LG_HCALOUT");
0291   assert(TOWER_RAW_LG_HCALOUT);
0292
0293   RawTowerContainer *TOWER_CALIB_LG_HCALOUT = findNode::getClass<
0294       RawTowerContainer>(topNode, "TOWER_CALIB_LG_HCALOUT");
0295   assert(TOWER_CALIB_LG_HCALOUT);
0296
0297   // HCALOUT HG information
0298   RawTowerContainer *TOWER_RAW_HG_HCALOUT = findNode::getClass<
0299       RawTowerContainer>(topNode, "TOWER_RAW_HG_HCALOUT");
0300   assert(TOWER_RAW_HG_HCALOUT);
0301
0302   RawTowerContainer *TOWER_CALIB_HG_HCALOUT = findNode::getClass<
0303       RawTowerContainer>(topNode, "TOWER_CALIB_HG_HCALOUT");
0304   assert(TOWER_CALIB_HG_HCALOUT);
0305
0306   // simple clustering by matching to cluster of max energy
0307   // pair<int, int> max_tower = find_max(TOWER_CALIB_LG_HCALOUT, 1);
0308
0309   // process SIM HCALIN/HCALOUT
0310   if(_is_sim)
0311   {
0312     double hcalin_sum_e_sim = 0;
0313     double hcalout_sum_e_sim = 0;
0314     {
0315       auto range_hcalin_sim = TOWER_SIM_HCALIN->getTowers(); // sim
0316       for (auto it = range_hcalin_sim.first; it != range_hcalin_sim.second; ++it)
0317       {
0318     RawTower *tower = it->second;
0319     assert(tower);
0320
0321     const int col = tower->get_bineta();
0322     const int row = tower->get_binphi();
0323     const int chanNum = getChannelNumber(row,col);
0324
0325     // cout << "HCALIN SIM: ";
0326     // tower->identify();
0327     // cout << "getSimChannelNumber = " << chanNum << endl;
0328
0329     if (col < 0 or col >= 4)
0330       continue;
0331     if (row < 0 or row >= 4)
0332       continue;
0333
0334     const double energy_sim = tower->get_energy();
0335     hcalin_sum_e_sim += energy_sim;
0336     _tower.hcalin_twr_sim[chanNum] = energy_sim;
0337
0338     string HistName_sim = Form("h_hcalin_tower_%d_sim",chanNum);
0339     TH1F *h_hcalin_sim = dynamic_cast<TH1F *>(hm->getHisto(HistName_sim.c_str()));
0340     assert(h_hcalin_sim);
0341     h_hcalin_sim->Fill(_tower.hcalin_twr_sim[chanNum]*1000.0); // convert to MeV
0342       }
0343       _tower.hcalin_e_sim = hcalin_sum_e_sim;
0344
0345       auto range_hcalout_sim = TOWER_SIM_HCALOUT->getTowers(); // sim
0346       for (auto it = range_hcalout_sim.first; it != range_hcalout_sim.second; ++it)
0347       {
0348     RawTower *tower = it->second;
0349     assert(tower);
0350     // cout << "HCALOUT SIM: ";
0351     // tower->identify();
0352
0353     const int col = tower->get_bineta();
0354     const int row = tower->get_binphi();
0355     const int chanNum = getChannelNumber(row,col);
0356
0357     if (col < 0 or col >= 4)
0358       continue;
0359     if (row < 0 or row >= 4)
0360       continue;
0361
0362     const double energy_sim = tower->get_energy();
0363     hcalout_sum_e_sim += energy_sim;
0364     _tower.hcalout_twr_sim[chanNum] = energy_sim;
0365
0366     string HistName_sim = Form("h_hcalout_tower_%d_sim",chanNum);
0367     TH1F *h_hcalout_sim = dynamic_cast<TH1F *>(hm->getHisto(HistName_sim.c_str()));
0368     assert(h_hcalout_sim);
0369     h_hcalout_sim->Fill(_tower.hcalout_twr_sim[chanNum]*1000.0); // convert to MeV
0370       }
0371       _tower.hcalout_e_sim = hcalout_sum_e_sim;
0372     }
0373     _tower.hcal_total_sim = hcalin_sum_e_sim + hcalout_sum_e_sim;
0374     _tower.hcal_asym_sim = (hcalin_sum_e_sim - hcalout_sum_e_sim)/(hcalin_sum_e_sim + hcalout_sum_e_sim);
0375   }
0376
0377   // process HCALIN LG
0378   double hcalin_sum_lg_e_raw = 0;
0379   double hcalin_sum_lg_e_calib = 0;
0380   {
0381     auto range_hcalin_lg_raw = TOWER_RAW_LG_HCALIN->getTowers(); // raw
0382     for (auto it = range_hcalin_lg_raw.first; it != range_hcalin_lg_raw.second; ++it)
0383     {
0384       RawTower *tower = it->second;
0385       assert(tower);
0386       const int col = tower->get_bineta();
0387       const int row = tower->get_binphi();
0388       const int chanNum = getChannelNumber(row,col);
0389
0390       // cout << "HCALIN RAW: ";
0391       // tower->identify();
0392       // cout << "getDataChannelNumber = " << chanNum << endl;
0393
0394       if (col < 0 or col >= 4)
0395         continue;
0396       if (row < 0 or row >= 4)
0397         continue;
0398
0399       const double energy_raw = tower->get_energy();
0400       hcalin_sum_lg_e_raw += energy_raw;
0401       _tower.hcalin_lg_twr_raw[chanNum] = energy_raw;
0402
0403       string HistName_raw = Form("h_hcalin_lg_tower_%d_raw",chanNum);
0404       TH1F *h_hcalin_lg_raw = dynamic_cast<TH1F *>(hm->getHisto(HistName_raw.c_str()));
0405       assert(h_hcalin_lg_raw);
0406       h_hcalin_lg_raw->Fill(_tower.hcalin_lg_twr_raw[chanNum]);
0407     }
0408     _tower.hcalin_lg_e_raw = hcalin_sum_lg_e_raw;
0409
0410     auto range_hcalin_lg_calib = TOWER_CALIB_LG_HCALIN->getTowers(); // calib
0411     for (auto it = range_hcalin_lg_calib.first; it != range_hcalin_lg_calib.second; ++it)
0412     {
0413       RawTower *tower = it->second;
0414       assert(tower);
0415
0416       const int col = tower->get_bineta();
0417       const int row = tower->get_binphi();
0418       const int chanNum = getChannelNumber(row,col);
0419
0420       if (col < 0 or col >= 4)
0421         continue;
0422       if (row < 0 or row >= 4)
0423         continue;
0424
0425       const double energy_calib = tower->get_energy();
0426       hcalin_sum_lg_e_calib += energy_calib;
0427       _tower.hcalin_lg_twr_calib[chanNum] = energy_calib;
0428
0429       string HistName_calib = Form("h_hcalin_lg_tower_%d_calib",chanNum);
0430       TH1F *h_hcalin_lg_calib = dynamic_cast<TH1F *>(hm->getHisto(HistName_calib.c_str()));
0431       assert(h_hcalin_lg_calib);
0432       h_hcalin_lg_calib->Fill(_tower.hcalin_lg_twr_calib[chanNum]);
0433     }
0434     _tower.hcalin_lg_e_calib = hcalin_sum_lg_e_calib;
0435   }
0436
0437   // process HCALOUT LG
0438   double hcalout_sum_lg_e_raw = 0;
0439   double hcalout_sum_lg_e_calib = 0;
0440   {
0441     auto range_hcalout_lg_raw = TOWER_RAW_LG_HCALOUT->getTowers(); // raw
0442     for (auto it = range_hcalout_lg_raw.first; it != range_hcalout_lg_raw.second; ++it)
0443     {
0444       RawTower *tower = it->second;
0445       assert(tower);
0446
0447       const int col = tower->get_bineta();
0448       const int row = tower->get_binphi();
0449       const int chanNum = getChannelNumber(row,col);
0450
0451       // cout << "HCALOUT LG RAW: ";
0452       // tower->identify();
0453       // cout << "getDataChannelNumber = " << chanNum << endl;
0454
0455       if (col < 0 or col >= 4)
0456         continue;
0457       if (row < 0 or row >= 4)
0458         continue;
0459
0460       const double energy_raw = tower->get_energy();
0461       hcalout_sum_lg_e_raw += energy_raw;
0462       _tower.hcalout_lg_twr_raw[chanNum] = energy_raw;
0463
0464       string HistName_raw = Form("h_hcalout_lg_tower_%d_raw",chanNum);
0465       TH1F *h_hcalout_lg_raw = dynamic_cast<TH1F *>(hm->getHisto(HistName_raw.c_str()));
0466       assert(h_hcalout_lg_raw);
0467       h_hcalout_lg_raw->Fill(_tower.hcalout_lg_twr_raw[chanNum]);
0468     }
0469     _tower.hcalout_lg_e_raw = hcalout_sum_lg_e_raw;
0470
0471     auto range_hcalout_lg_calib = TOWER_CALIB_LG_HCALOUT->getTowers(); // calib
0472     for (auto it = range_hcalout_lg_calib.first; it != range_hcalout_lg_calib.second; ++it)
0473     {
0474       RawTower *tower = it->second;
0475       assert(tower);
0476
0477       const int col = tower->get_bineta();
0478       const int row = tower->get_binphi();
0479       const int chanNum = getChannelNumber(row,col);
0480
0481       if (col < 0 or col >= 4)
0482         continue;
0483       if (row < 0 or row >= 4)
0484         continue;
0485
0486       const double energy_calib = tower->get_energy();
0487       hcalout_sum_lg_e_calib += energy_calib;
0488       _tower.hcalout_lg_twr_calib[chanNum] = energy_calib;
0489
0490       string HistName_calib = Form("h_hcalout_lg_tower_%d_calib",chanNum);
0491       TH1F *h_hcalout_lg_calib = dynamic_cast<TH1F *>(hm->getHisto(HistName_calib.c_str()));
0492       assert(h_hcalout_lg_calib);
0493       h_hcalout_lg_calib->Fill(_tower.hcalout_lg_twr_calib[chanNum]);
0494     }
0495     _tower.hcalout_lg_e_calib = hcalout_sum_lg_e_calib;
0496   }
0497
0498   // process HCALOUT HG
0499   double hcalout_sum_hg_e_raw = 0;
0500   double hcalout_sum_hg_e_calib = 0;
0501   {
0502     auto range_hcalout_hg_raw = TOWER_RAW_HG_HCALOUT->getTowers(); // raw
0503     for (auto it = range_hcalout_hg_raw.first; it != range_hcalout_hg_raw.second; ++it)
0504     {
0505       RawTower *tower = it->second;
0506       assert(tower);
0507
0508       const int col = tower->get_bineta();
0509       const int row = tower->get_binphi();
0510       const int chanNum = getChannelNumber(row,col);
0511
0512       // cout << "HCALOUT HG RAW: ";
0513       // tower->identify();
0514       // cout << "getDataChannelNumber = " << chanNum << endl;
0515
0516       if (col < 0 or col >= 4)
0517         continue;
0518       if (row < 0 or row >= 4)
0519         continue;
0520
0521       const double energy_raw = tower->get_energy();
0522       hcalout_sum_hg_e_raw += energy_raw;
0523       _tower.hcalout_hg_twr_raw[chanNum] = energy_raw;
0524
0525       string HistName_raw = Form("h_hcalout_hg_tower_%d_raw",chanNum);
0526       TH1F *h_hcalout_hg_raw = dynamic_cast<TH1F *>(hm->getHisto(HistName_raw.c_str()));
0527       assert(h_hcalout_hg_raw);
0528       h_hcalout_hg_raw->Fill(_tower.hcalout_hg_twr_raw[chanNum]);
0529
0530     }
0531     _tower.hcalout_hg_e_raw = hcalout_sum_hg_e_raw;
0532
0533     auto range_hcalout_hg_calib = TOWER_CALIB_HG_HCALOUT->getTowers(); // calib
0534     for (auto it = range_hcalout_hg_calib.first; it != range_hcalout_hg_calib.second; ++it)
0535     {
0536       RawTower *tower = it->second;
0537       assert(tower);
0538
0539       const int col = tower->get_bineta();
0540       const int row = tower->get_binphi();
0541       const int chanNum = getChannelNumber(row,col);
0542
0543       if (col < 0 or col >= 4)
0544         continue;
0545       if (row < 0 or row >= 4)
0546         continue;
0547
0548       const double energy_calib = tower->get_energy();
0549       hcalout_sum_hg_e_calib += energy_calib;
0550       _tower.hcalout_hg_twr_calib[chanNum] = energy_calib;
0551
0552       string HistName_calib = Form("h_hcalout_hg_tower_%d_calib",chanNum);
0553       TH1F *h_hcalout_hg_calib = dynamic_cast<TH1F *>(hm->getHisto(HistName_calib.c_str()));
0554       assert(h_hcalout_hg_calib);
0555       h_hcalout_hg_calib->Fill(_tower.hcalout_hg_twr_calib[chanNum]);
0556     }
0557     _tower.hcalout_hg_e_calib = hcalout_sum_hg_e_calib;
0558   }
0559
0560   _tower.hcal_total_raw = hcalin_sum_lg_e_raw + hcalout_sum_lg_e_raw;
0561   _tower.hcal_total_calib = hcalin_sum_lg_e_calib + hcalout_sum_lg_e_calib;
0562
0563   _tower.hcal_asym_raw = (hcalin_sum_lg_e_raw - hcalout_sum_lg_e_raw)/(hcalin_sum_lg_e_raw + hcalout_sum_lg_e_raw);
0564   _tower.hcal_asym_calib = (hcalin_sum_lg_e_calib - hcalout_sum_lg_e_calib)/(hcalin_sum_lg_e_calib + hcalout_sum_lg_e_calib);
0565
0566   TTree *T = dynamic_cast<TTree *>(hm->getHisto("HCAL_Info"));
0567   assert(T);
0568   T->Fill();
0569
0570   return Fun4AllReturnCodes::EVENT_OK;
0571 }
0572
0573 pair<int, int>
0574 Proto4TowerCalib::find_max(RawTowerContainer *towers, int cluster_size)
0575 {
0576   const int clus_edge_size = (cluster_size - 1) / 2;
0577   assert(clus_edge_size >= 0);
0578
0579   pair<int, int> max(-1000, -1000);
0580   double max_e = 0;
0581
0582   for (int col = 0; col < n_size; ++col)
0583     for (int row = 0; row < n_size; ++row)
0584     {
0585       double energy = 0;
0586
0587       for (int dcol = col - clus_edge_size; dcol <= col + clus_edge_size;
0588            ++dcol)
0589         for (int drow = row - clus_edge_size; drow <= row + clus_edge_size;
0590              ++drow)
0591         {
0592           if (dcol < 0 or drow < 0)
0593             continue;
0594
0595           RawTower *t = towers->getTower(dcol, drow);
0596           if (t)
0597             energy += t->get_energy();
0598         }
0599
0600       if (energy > max_e)
0601       {
0602         max = make_pair(col, row);
0603         max_e = energy;
0604       }
0605     }
0606
0607   return max;
0608 }
0609
0610 int Proto4TowerCalib::getChannelNumber(int row, int column)
0611 {
0612   if(!_is_sim)
0613   {
0614     int hbdchanIHC[4][4] = {{4, 8, 12, 16},
0615                             {3, 7, 11, 15},
0616                             {2, 6, 10, 14},
0617                             {1, 5,  9, 13}};
0618
0619     return hbdchanIHC[row][column] - 1;
0620   }
0621
0622   if(_is_sim)
0623   {
0624     int hbdchanIHC[4][4] = {{13,  9, 5, 1},
0625                             {14, 10, 6, 2},
0626                             {15, 11, 7, 3},
0627                             {16, 12, 8, 4}};
0628
0629     return hbdchanIHC[row][column] - 1;
0630   }
0631
0632   return -1;
0633 }
0634
0635 int Proto4TowerCalib::InitAna()
0636 {
0637   if(_is_sim) _mList = "SIM";
0638   if(!_is_sim) _mList = "RAW";
0639   string inputdir = "/sphenix/user/xusun/TestBeam/TowerCalib/";
0640   string InPutList = Form("/direct/phenix+u/xusun/WorkSpace/sPHENIX/analysis/Prototype4/HCal/macros/list/Proto4TowerInfo%s_%s_%d.list",_mList.c_str(),_mDet.c_str(),_mCol);
0641
0642   mChainInPut = new TChain("HCAL_Info");
0643
0644   if (!InPutList.empty())   // if input file is ok
0645   {
0646     cout << "Open input database file list: " << InPutList.c_str() << endl;
0647     ifstream in(InPutList.c_str());  // input stream
0648     if(in)
0649     {
0650       cout << "input database file list is ok" << endl;
0651       char str[255];       // char array for each file name
0652       Long64_t entries_save = 0;
0653       while(in)
0654       {
0655     in.getline(str,255);  // take the lines of the file list
0656     if(str[0] != 0)
0657     {
0658       string addfile;
0659       addfile = str;
0660       addfile = inputdir+addfile;
0661       mChainInPut->AddFile(addfile.c_str(),-1,"HCAL_Info");
0662       long file_entries = mChainInPut->GetEntries();
0663       cout << "File added to data chain: " << addfile.c_str() << " with " << (file_entries-entries_save) << " entries" << endl;
0664       entries_save = file_entries;
0665     }
0666       }
0667     }
0668     else
0669     {
0670       cout << "WARNING: input database file input is problemtic" << endl;
0671       _mInPut_flag = 0;
0672     }
0673   }
0674
0675   // Set the input tree
0676   if (_mInPut_flag == 1 && !mChainInPut->GetBranch( "TowerCalib" ))
0677   {
0678     cerr << "ERROR: Could not find branch 'HCAL_Tower' in tree!" << endl;
0679   }
0680
0681   _mTower = new Proto4TowerCalib::HCAL_Tower();
0682   if(_mInPut_flag == 1)
0683   {
0684     mChainInPut->SetBranchAddress("TowerCalib", &_mTower);
0685
0686     long NumOfEvents = (long)mChainInPut->GetEntries();
0687     cout << "total number of events: " << NumOfEvents << endl;
0688     _mStartEvent = 0;
0689     _mStopEvent = NumOfEvents;
0690   }
0691
0692   // initialize TH1Fs for energy/adc spectra
0693   for(int i_tower = 0; i_tower < 16; ++i_tower)
0694   {
0695     std::string HistName = Form("h_m%s_%s_twr_%d",_mDet.c_str(),_mList.c_str(),i_tower);
0696     if(_is_sim) h_mHCAL[i_tower] = new TH1F(HistName.c_str(),HistName.c_str(),500,0.5,100.5);
0697     if(!_is_sim) h_mHCAL[i_tower] = new TH1F(HistName.c_str(),HistName.c_str(),80,0.5,16000.5);
0698     // if(!_is_sim) h_mHCAL[i_tower] = new TH1F(HistName.c_str(),HistName.c_str(),500,0.5,100.5);
0699   }
0700
0701   return 0;
0702 }
0703
0704 int Proto4TowerCalib::MakeAna()
0705 {
0706   cout << "Make()" << endl;
0707   unsigned long start_event_use = _mStartEvent;
0708   unsigned long stop_event_use = _mStopEvent;
0709
0710   mChainInPut->SetBranchAddress("TowerCalib", &_mTower);
0711   mChainInPut->GetEntry(0);
0712
0713   for(unsigned long i_event = start_event_use; i_event < stop_event_use; ++i_event)
0714   // for(unsigned long i_event = 20; i_event < 40; ++i_event)
0715   {
0716     if (!mChainInPut->GetEntry( i_event )) // take the event -> information is stored in event
0717       break;  // end of data chunk
0718
0719     if(_is_sim) // cosmic simulation
0720     {
0721       for(int i_tower = 0; i_tower < 16; ++i_tower)
0722       {
0723     if(_mDet == "HCALIN")  // HCALIN
0724     {
0725       hcal_twr[i_tower] = _mTower->hcalin_twr_sim[i_tower];
0726       h_mHCAL[i_tower]->Fill(hcal_twr[i_tower]*1000.0); // convert to MeV
0727     }
0728     if(_mDet == "HCALOUT")  // HCALOUT
0729     {
0730       hcal_twr[i_tower] = _mTower->hcalout_twr_sim[i_tower];
0731       h_mHCAL[i_tower]->Fill(hcal_twr[i_tower]*1000.0); // convert to MeV
0732     }
0733       }
0734     }
0735     if(!_is_sim) // cosmic data
0736     {
0737       for(int i_tower = 0; i_tower < 16; ++i_tower) // 1st loop to decide pedestal
0738       {
0739     if(_mDet == "HCALIN") // HCALIN
0740     {
0741       hcal_twr[i_tower] = _mTower->hcalin_lg_twr_raw[i_tower];
0742       if( hcal_twr[i_tower] > _mPedestal ) _is_sig[i_tower] = true;
0743     }
0744     if(_mDet == "HCALOUT") // HCALOUT
0745     {
0746       hcal_twr[i_tower] = _mTower->hcalout_hg_twr_raw[i_tower];
0747       if(hcal_twr[i_tower] > _mPedestal) _is_sig[i_tower] = true;
0748     }
0749       }
0750       if( is_sig(_mCol) )
0751       {
0752     // cout << "i_event = " << i_event << ", is_sig " << is_sig(_mCol) << endl;
0753     fill_sig(_mCol);
0754       }
0755       reset_pedestal();
0756     }
0757   }
0758
0759   return 1;
0760 }
0761
0762 int Proto4TowerCalib::FinishAna()
0763 {
0764   cout << "Finish()" << endl;
0765   mFile_OutPut = new TFile(_filename.c_str(),"RECREATE");
0766   mFile_OutPut->cd();
0767   for(int i_tower = 0; i_tower < 16; ++i_tower)
0768   {
0769     h_mHCAL[i_tower]->Write();
0770   }
0771   mFile_OutPut->Close();
0772
0773   return 1;
0774 }
0775
0776 bool Proto4TowerCalib::is_sig(int colID)
0777 {
0778   int twr_0 = getChannelNumber(0,colID);
0779   int twr_1 = getChannelNumber(1,colID);
0780   int twr_2 = getChannelNumber(2,colID);
0781   int twr_3 = getChannelNumber(3,colID);
0782
0783   if(_is_sig[twr_0] && _is_sig[twr_1] && _is_sig[twr_2]) return true;
0784   if(_is_sig[twr_0] && _is_sig[twr_1] && _is_sig[twr_3]) return true;
0785   if(_is_sig[twr_0] && _is_sig[twr_2] && _is_sig[twr_3]) return true;
0786   if(_is_sig[twr_1] && _is_sig[twr_2] && _is_sig[twr_3]) return true;
0787
0788   return false;
0789 }
0790
0791 int Proto4TowerCalib::fill_sig(int colID)
0792 {
0793   for(int i_row = 0; i_row < 4; ++i_row)
0794   {
0795     // fill adc spectra for cosmic
0796     int i_tower = getChannelNumber(i_row,colID);
0797     h_mHCAL[i_tower]->Fill(hcal_twr[i_tower]);
0798   }
0799
0800   return 1;
0801 }
0802
0803 int Proto4TowerCalib::reset_pedestal()
0804 {
0805   for(int i_tower = 0; i_tower < 16; ++i_tower)
0806   {
0807     hcal_twr[i_tower] = -1.0; // set tower energy
0808     _is_sig[i_tower] = false; // set pedestal cut
0809   }
0810
0811   return 1;
0812 }