Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:14:49 

0001 #include "Proto4ShowerCalib.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 <TH1F.h>
0027 #include <TH2F.h>
0028 #include <TLorentzVector.h>
0029 #include <TNtuple.h>
0030 #include <TVector3.h>
0031 #include <TChain.h>
0032 
0033 #include <algorithm>
0034 #include <cassert>
0035 #include <cmath>
0036 #include <exception>
0037 #include <iostream>
0038 #include <set>
0039 #include <sstream>
0040 #include <stdexcept>
0041 #include <vector>
0042 
0043 using namespace std;
0044 
0045 ClassImp(Proto4ShowerCalib::HCAL_Tower);
0046 ClassImp(Proto4ShowerCalib::Eval_Run);
0047 
0048 Proto4ShowerCalib::Proto4ShowerCalib(const std::string &filename)
0049   : SubsysReco("Proto4ShowerCalib")
0050   , _is_sim(false)
0051   , _filename(filename)
0052   , _ievent(0)
0053 {
0054   Verbosity(1);
0055 
0056   _eval_run.reset();
0057   _tower.reset();
0058 
0059   _mInPut_flag = 1;
0060   _mStartEvent = -1;
0061   _mStopEvent = -1;
0062 
0063   setTowerCalibParas();
0064   _mRunID = "0422";
0065 }
0066 
0067 Proto4ShowerCalib::~Proto4ShowerCalib()
0068 {
0069 }
0070 
0071 Fun4AllHistoManager *
0072 Proto4ShowerCalib::get_HistoManager()
0073 {
0074   Fun4AllServer *se = Fun4AllServer::instance();
0075   Fun4AllHistoManager *hm = se->getHistoManager("Proto4ShowerCalib_HISTOS");
0076 
0077   if (not hm)
0078   {
0079     cout
0080         << "Proto4ShowerCalib::get_HistoManager - Making Fun4AllHistoManager Proto4ShowerCalib_HISTOS"
0081         << endl;
0082     hm = new Fun4AllHistoManager("Proto4ShowerCalib_HISTOS");
0083     se->registerHistoManager(hm);
0084   }
0085 
0086   assert(hm);
0087 
0088   return hm;
0089 }
0090 
0091 int Proto4ShowerCalib::InitRun(PHCompositeNode *topNode)
0092 {
0093   if (Verbosity())
0094     cout << "Proto4ShowerCalib::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 Proto4ShowerCalib::End(PHCompositeNode *topNode)
0112 {
0113   cout << "Proto4ShowerCalib::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 Proto4ShowerCalib::Init(PHCompositeNode *topNode)
0128 {
0129   _ievent = 0;
0130 
0131   cout << "Proto4ShowerCalib::get_HistoManager - Making PHTFileServer "
0132        << _filename << endl;
0133   PHTFileServer::get().open(_filename, "RECREATE");
0134 
0135   Fun4AllHistoManager *hm = get_HistoManager();
0136   assert(hm);
0137 
0138   //! Histogram of Cherenkov counters
0139   TH2F *hCheck_Cherenkov = new TH2F("hCheck_Cherenkov", "hCheck_Cherenkov",
0140                                     110, -200, 20000, 5, .5, 5.5);
0141   hCheck_Cherenkov->GetYaxis()->SetBinLabel(1, "C1");
0142   hCheck_Cherenkov->GetYaxis()->SetBinLabel(2, "C2 in");
0143   hCheck_Cherenkov->GetYaxis()->SetBinLabel(3, "C2 out");
0144   hCheck_Cherenkov->GetYaxis()->SetBinLabel(4, "C2 sum");
0145   hCheck_Cherenkov->GetXaxis()->SetTitle("Amplitude");
0146   hCheck_Cherenkov->GetYaxis()->SetTitle("Cherenkov Channel");
0147   hm->registerHisto(hCheck_Cherenkov);
0148 
0149   //! Envent nomalization
0150   TH1F *hNormalization = new TH1F("hNormalization", "hNormalization", 10, .5,
0151                                   10.5);
0152   hNormalization->GetXaxis()->SetBinLabel(1, "ALL");
0153   hNormalization->GetXaxis()->SetBinLabel(2, "C2-e");
0154   hNormalization->GetXaxis()->SetBinLabel(3, "trigger_veto_pass");
0155   hNormalization->GetXaxis()->SetBinLabel(4, "valid_hodo_h");
0156   hNormalization->GetXaxis()->SetBinLabel(5, "valid_hodo_v");
0157   hNormalization->GetXaxis()->SetBinLabel(6, "good_e");
0158   hNormalization->GetXaxis()->SetBinLabel(7, "good_anti_e");
0159   hNormalization->GetXaxis()->SetTitle("Cuts");
0160   hNormalization->GetYaxis()->SetTitle("Event count");
0161   hm->registerHisto(hNormalization);
0162 
0163   hm->registerHisto(new TH1F("hCheck_Veto", "hCheck_Veto", 1000, -500, 1500));
0164   hm->registerHisto(
0165       new TH1F("hCheck_Hodo_H", "hCheck_Hodo_H", 1000, -500, 1500));
0166   hm->registerHisto(
0167       new TH1F("hCheck_Hodo_V", "hCheck_Hodo_V", 1000, -500, 1500));
0168 
0169   hm->registerHisto(new TH1F("hBeam_Mom", "hBeam_Mom", 1200, -120, 120));
0170 
0171   float histo_range = find_range();
0172   // cout << "histo_range = " << histo_range << endl;
0173   // float histo_range = 49.5;
0174 
0175   // SIM HCALIN/HCALOUT
0176   if(_is_sim)
0177   {
0178     TH1F *h_hcalin_tower_sim[16];
0179     TH1F *h_hcalout_tower_sim[16];
0180     for(int i_tower = 0; i_tower < 16; ++i_tower)
0181     {
0182       string HistName_sim;
0183       HistName_sim = Form("h_hcalin_tower_%d_sim",i_tower);
0184       // h_hcalin_tower_sim[i_tower] = new TH1F(HistName_sim.c_str(),HistName_sim.c_str(),100,-0.5,9.5);
0185       h_hcalin_tower_sim[i_tower] = new TH1F(HistName_sim.c_str(),HistName_sim.c_str(),100,-0.5,histo_range);
0186       hm->registerHisto(h_hcalin_tower_sim[i_tower]);
0187 
0188       HistName_sim = Form("h_hcalout_tower_%d_sim",i_tower);
0189       // h_hcalout_tower_sim[i_tower] = new TH1F(HistName_sim.c_str(),HistName_sim.c_str(),100,-0.5,9.5);
0190       h_hcalout_tower_sim[i_tower] = new TH1F(HistName_sim.c_str(),HistName_sim.c_str(),100,-0.5,histo_range);
0191       hm->registerHisto(h_hcalout_tower_sim[i_tower]);
0192     }
0193 
0194     TH2F *h_tower_energy_sim = new TH2F("h_tower_energy_sim","h_tower_energy_sim",100,-1.0,1.0,100,-0.5,histo_range);
0195     hm->registerHisto(h_tower_energy_sim);
0196   }
0197 
0198   // HCALIN LG
0199   TH1F *h_hcalin_lg_tower_raw[16];
0200   TH1F *h_hcalin_lg_tower_calib[16];
0201   for(int i_tower = 0; i_tower < 16; ++i_tower)
0202   {
0203     string HistName_raw = Form("h_hcalin_lg_tower_%d_raw",i_tower);
0204     h_hcalin_lg_tower_raw[i_tower] = new TH1F(HistName_raw.c_str(),HistName_raw.c_str(),100,-0.5,histo_range*20.0);
0205     hm->registerHisto(h_hcalin_lg_tower_raw[i_tower]);
0206 
0207     string HistName_calib = Form("h_hcalin_lg_tower_%d_calib",i_tower);
0208     h_hcalin_lg_tower_calib[i_tower] = new TH1F(HistName_calib.c_str(),HistName_calib.c_str(),100,-0.5,histo_range);
0209     hm->registerHisto(h_hcalin_lg_tower_calib[i_tower]);
0210   }
0211 
0212   // HCALOUT LG
0213   TH1F *h_hcalout_lg_tower_raw[16];
0214   TH1F *h_hcalout_lg_tower_calib[16];
0215   for(int i_tower = 0; i_tower < 16; ++i_tower)
0216   {
0217     string HistName_raw = Form("h_hcalout_lg_tower_%d_raw",i_tower);
0218     h_hcalout_lg_tower_raw[i_tower] = new TH1F(HistName_raw.c_str(),HistName_raw.c_str(),100,-0.5,histo_range*20.0);
0219     hm->registerHisto(h_hcalout_lg_tower_raw[i_tower]);
0220 
0221     string HistName_calib = Form("h_hcalout_lg_tower_%d_calib",i_tower);
0222     h_hcalout_lg_tower_calib[i_tower] = new TH1F(HistName_calib.c_str(),HistName_calib.c_str(),100,-0.5,histo_range);
0223     hm->registerHisto(h_hcalout_lg_tower_calib[i_tower]);
0224   }
0225 
0226   TH2F *h_tower_energy_raw = new TH2F("h_tower_energy_raw","h_tower_energy_raw",110,-1.1,1.1,100,-0.5,histo_range*640.0);
0227   hm->registerHisto(h_tower_energy_raw);
0228 
0229   TH2F *h_tower_energy_calib = new TH2F("h_tower_energy_calib","h_tower_energy_calib",110,-1.1,1.1,100,-0.5,histo_range);
0230   hm->registerHisto(h_tower_energy_calib);
0231 
0232   // HCALOUT HG
0233   TH1F *h_hcalout_hg_tower_raw[16];
0234   TH1F *h_hcalout_hg_tower_calib[16];
0235   for(int i_tower = 0; i_tower < 16; ++i_tower)
0236   {
0237     string HistName_raw = Form("h_hcalout_hg_tower_%d_raw",i_tower);
0238     h_hcalout_hg_tower_raw[i_tower] = new TH1F(HistName_raw.c_str(),HistName_raw.c_str(),100,-0.5,histo_range*20.0);
0239     hm->registerHisto(h_hcalout_hg_tower_raw[i_tower]);
0240 
0241     string HistName_calib = Form("h_hcalout_hg_tower_%d_calib",i_tower);
0242     h_hcalout_hg_tower_calib[i_tower] = new TH1F(HistName_calib.c_str(),HistName_calib.c_str(),100,-0.5,histo_range);
0243     hm->registerHisto(h_hcalout_hg_tower_calib[i_tower]);
0244   }
0245 
0246 
0247   // help index files with TChain
0248   TTree *T = new TTree("HCAL_Info", "HCAL_Info");
0249   assert(T);
0250   hm->registerHisto(T);
0251 
0252   T->Branch("info", &_eval_run);
0253   T->Branch("tower", &_tower);
0254 
0255   return Fun4AllReturnCodes::EVENT_OK;
0256 }
0257 
0258 int Proto4ShowerCalib::process_event(PHCompositeNode *topNode)
0259 {
0260   if (Verbosity() > 2)
0261     cout << "Proto4ShowerCalib::process_event() entered" << endl;
0262 
0263   // init eval objects
0264   _eval_run.reset();
0265   _tower.reset();
0266 
0267   Fun4AllHistoManager *hm = get_HistoManager();
0268   assert(hm);
0269 
0270   if (not _is_sim)
0271   {
0272     PdbParameterMap *info = findNode::getClass<PdbParameterMap>(topNode,
0273                                                                 "RUN_INFO");
0274 
0275     assert(info);
0276 
0277     PHParameters run_info_copy("RunInfo");
0278     run_info_copy.FillFrom(info);
0279 
0280     _eval_run.beam_mom = run_info_copy.get_double_param("beam_MTNRG_GeV");
0281 
0282     TH1F *hBeam_Mom = dynamic_cast<TH1F *>(hm->getHisto("hBeam_Mom"));
0283     assert(hBeam_Mom);
0284 
0285     hBeam_Mom->Fill(_eval_run.beam_mom);
0286 
0287     _eval_run.beam_2CH_mm = run_info_copy.get_double_param("beam_2CH_mm");
0288     _eval_run.beam_2CV_mm = run_info_copy.get_double_param("beam_2CV_mm");
0289   }
0290 
0291   if (not _is_sim)
0292   {
0293     EventHeader *eventheader = findNode::getClass<EventHeader>(topNode,
0294     "EventHeader");
0295 
0296     if (eventheader->get_EvtType() != DATAEVENT)
0297     {
0298       cout << __PRETTY_FUNCTION__
0299     << " disgard this event: " << endl;
0300 
0301       eventheader->identify();
0302 
0303       return Fun4AllReturnCodes::DISCARDEVENT;
0304     }
0305 
0306     assert(eventheader);
0307 
0308     _eval_run.run = eventheader->get_RunNumber();
0309     if (Verbosity() > 4)
0310       cout << __PRETTY_FUNCTION__ << _eval_run.run << endl;
0311 
0312     _eval_run.event = eventheader->get_EvtSequence();
0313   }
0314 
0315   if (_is_sim)
0316   {
0317     PHG4TruthInfoContainer *truthInfoList = findNode::getClass<
0318         PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0319 
0320     assert(truthInfoList);
0321 
0322     _eval_run.run = -1;
0323 
0324     const PHG4Particle *p =
0325         truthInfoList->GetPrimaryParticleRange().first->second;
0326     assert(p);
0327 
0328     const PHG4VtxPoint *v = truthInfoList->GetVtx(p->get_vtx_id());
0329     assert(v);
0330 
0331     _eval_run.beam_mom = sqrt(
0332         p->get_px() * p->get_px() + p->get_py() * p->get_py() + p->get_pz() * p->get_pz());
0333     _eval_run.truth_y = v->get_y();
0334     _eval_run.truth_z = v->get_z();
0335   }
0336 
0337   // normalization
0338   TH1F *hNormalization = dynamic_cast<TH1F *>(hm->getHisto("hNormalization"));
0339   assert(hNormalization);
0340 
0341   hNormalization->Fill("ALL", 1);
0342 
0343   // get DST nodes
0344 
0345   // HCALIN/HCALOUT SIM information
0346   RawTowerContainer *TOWER_SIM_HCALIN = findNode::getClass<
0347       RawTowerContainer>(topNode, "TOWER_SIM_HCALIN");
0348   if(_is_sim)
0349   {
0350     assert(TOWER_SIM_HCALIN);
0351   }
0352 
0353   RawTowerContainer *TOWER_SIM_HCALOUT = findNode::getClass<
0354       RawTowerContainer>(topNode, "TOWER_SIM_HCALOUT");
0355   if(_is_sim)
0356   {
0357     assert(TOWER_SIM_HCALOUT);
0358   }
0359 
0360   // HCAL information
0361   RawTowerContainer *TOWER_RAW_LG_HCALIN = findNode::getClass<
0362       RawTowerContainer>(topNode, "TOWER_RAW_LG_HCALIN");
0363   assert(TOWER_RAW_LG_HCALIN);
0364 
0365   RawTowerContainer *TOWER_RAW_LG_HCALOUT = findNode::getClass<
0366       RawTowerContainer>(topNode, "TOWER_RAW_LG_HCALOUT");
0367   assert(TOWER_RAW_LG_HCALOUT);
0368 
0369   RawTowerContainer *TOWER_RAW_HG_HCALOUT = findNode::getClass<
0370       RawTowerContainer>(topNode, "TOWER_RAW_HG_HCALOUT");
0371   assert(TOWER_RAW_HG_HCALOUT);
0372 
0373   RawTowerContainer *TOWER_CALIB_LG_HCALIN = findNode::getClass<
0374       RawTowerContainer>(topNode, "TOWER_CALIB_LG_HCALIN");
0375   assert(TOWER_CALIB_LG_HCALIN);
0376 
0377   RawTowerContainer *TOWER_CALIB_LG_HCALOUT = findNode::getClass<
0378       RawTowerContainer>(topNode, "TOWER_CALIB_LG_HCALOUT");
0379   assert(TOWER_CALIB_LG_HCALOUT);
0380 
0381   RawTowerContainer *TOWER_CALIB_HG_HCALOUT = findNode::getClass<
0382       RawTowerContainer>(topNode, "TOWER_CALIB_HG_HCALOUT");
0383   assert(TOWER_CALIB_HG_HCALOUT);
0384 
0385   // other nodes
0386   RawTowerContainer *TOWER_CALIB_TRIGGER_VETO = findNode::getClass<
0387       RawTowerContainer>(topNode, "TOWER_CALIB_TRIGGER_VETO");
0388   if (not _is_sim)
0389   {
0390     assert(TOWER_CALIB_TRIGGER_VETO);
0391   }
0392 
0393   RawTowerContainer *TOWER_CALIB_HODO_HORIZONTAL = findNode::getClass<
0394       RawTowerContainer>(topNode, "TOWER_CALIB_HODO_HORIZONTAL");
0395   if (not _is_sim)
0396   {
0397     assert(TOWER_CALIB_HODO_HORIZONTAL);
0398   }
0399   RawTowerContainer *TOWER_CALIB_HODO_VERTICAL = findNode::getClass<
0400       RawTowerContainer>(topNode, "TOWER_CALIB_HODO_VERTICAL");
0401   if (not _is_sim)
0402   {
0403     assert(TOWER_CALIB_HODO_VERTICAL);
0404   }
0405 
0406   RawTowerContainer *TOWER_CALIB_C1 = findNode::getClass<RawTowerContainer>(
0407       topNode, "TOWER_CALIB_C1");
0408   if (not _is_sim)
0409   {
0410     assert(TOWER_CALIB_C1);
0411   }
0412   RawTowerContainer *TOWER_CALIB_C2 = findNode::getClass<RawTowerContainer>(
0413       topNode, "TOWER_CALIB_C2");
0414   if (not _is_sim)
0415   {
0416     assert(TOWER_CALIB_C2);
0417   }
0418 
0419   // Cherenkov
0420   bool cherekov_e = false;
0421   bool cherekov_anti_e = false;
0422   if (not _is_sim)
0423   {
0424     RawTower *t_c2_in = NULL;
0425     RawTower *t_c2_out = NULL;
0426 
0427     t_c2_in = TOWER_CALIB_C2->getTower(0);
0428     t_c2_out = TOWER_CALIB_C2->getTower(1);
0429 
0430     assert(t_c2_in);
0431     assert(t_c2_out);
0432 
0433     const double c2_in = t_c2_in->get_energy();
0434     const double c2_out = t_c2_out->get_energy();
0435     const double c1 = TOWER_CALIB_C1->getTower(0)->get_energy();
0436 
0437     _eval_run.C2_sum = c2_in + c2_out;
0438     _eval_run.C1 = c1;
0439 
0440     cherekov_e = (_eval_run.C2_sum) > (abs(_eval_run.beam_mom) >= 10 ? 1000 : 1500);
0441     cherekov_anti_e = (_eval_run.C2_sum) < 100;
0442     hNormalization->Fill("C2-e", cherekov_e);
0443 
0444     TH2F *hCheck_Cherenkov = dynamic_cast<TH2F *>(hm->getHisto(
0445         "hCheck_Cherenkov"));
0446     assert(hCheck_Cherenkov);
0447     hCheck_Cherenkov->Fill(c1, "C1", 1);
0448     hCheck_Cherenkov->Fill(c2_in, "C2 in", 1);
0449     hCheck_Cherenkov->Fill(c2_out, "C2 out", 1);
0450     hCheck_Cherenkov->Fill(c2_in + c2_out, "C2 sum", 1);
0451   }
0452 
0453   // veto
0454   TH1F *hCheck_Veto = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Veto"));
0455   assert(hCheck_Veto);
0456   bool trigger_veto_pass = true;
0457   if (not _is_sim)
0458   {
0459     auto range = TOWER_CALIB_TRIGGER_VETO->getTowers();
0460     for (auto it = range.first; it != range.second; ++it)
0461     {
0462       RawTower *tower = it->second;
0463       assert(tower);
0464 
0465       hCheck_Veto->Fill(tower->get_energy());
0466 
0467       if (abs(tower->get_energy()) > 0.2)
0468         trigger_veto_pass = false;
0469     }
0470   }
0471   hNormalization->Fill("trigger_veto_pass", trigger_veto_pass);
0472   _eval_run.trigger_veto_pass = trigger_veto_pass;
0473 
0474   // hodoscope
0475   TH1F *hCheck_Hodo_H = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_H"));
0476   assert(hCheck_Hodo_H);
0477   int hodo_h_count = 0;
0478   if (not _is_sim)
0479   {
0480     auto range = TOWER_CALIB_HODO_HORIZONTAL->getTowers();
0481     for (auto it = range.first; it != range.second; ++it)
0482     {
0483       RawTower *tower = it->second;
0484       assert(tower);
0485 
0486       hCheck_Hodo_H->Fill(tower->get_energy());
0487 
0488       if (abs(tower->get_energy()) >  0.5)
0489       {
0490         hodo_h_count++;
0491         _eval_run.hodo_h = tower->get_id();
0492       }
0493     }
0494   }
0495   const bool valid_hodo_h = hodo_h_count == 1;
0496   hNormalization->Fill("valid_hodo_h", valid_hodo_h);
0497   _eval_run.valid_hodo_h = valid_hodo_h;
0498 
0499   TH1F *hCheck_Hodo_V = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_V"));
0500   assert(hCheck_Hodo_V);
0501   int hodo_v_count = 0;
0502   if (not _is_sim)
0503   {
0504     auto range = TOWER_CALIB_HODO_VERTICAL->getTowers();
0505     for (auto it = range.first; it != range.second; ++it)
0506     {
0507       RawTower *tower = it->second;
0508       assert(tower);
0509 
0510       hCheck_Hodo_V->Fill(tower->get_energy());
0511 
0512       if (abs(tower->get_energy()) >  0.5)
0513       {
0514         hodo_v_count++;
0515         _eval_run.hodo_v = tower->get_id();
0516       }
0517     }
0518   }
0519   const bool valid_hodo_v = hodo_v_count == 1;
0520   _eval_run.valid_hodo_v = valid_hodo_v;
0521   hNormalization->Fill("valid_hodo_v", valid_hodo_v);
0522 
0523   const bool good_e = (valid_hodo_v and valid_hodo_h and cherekov_e and trigger_veto_pass) and (not _is_sim);
0524   const bool good_anti_e = (valid_hodo_v and valid_hodo_h and cherekov_anti_e and trigger_veto_pass) and (not _is_sim);
0525   hNormalization->Fill("good_e", good_e);
0526   hNormalization->Fill("good_anti_e", good_anti_e);
0527 
0528   // simple clustering by matching to cluster of max energy
0529   // pair<int, int> max_tower = find_max(TOWER_CALIB_LG_HCALOUT, 1);
0530 
0531   // process SIM HCALIN/HCALOUT
0532   if(_is_sim)
0533   {
0534     double hcalin_sum_e_sim = 0;
0535     double hcalout_sum_e_sim = 0;
0536     {
0537       auto range_hcalin_sim = TOWER_SIM_HCALIN->getTowers(); // sim 
0538       for (auto it = range_hcalin_sim.first; it != range_hcalin_sim.second; ++it)
0539       {
0540     RawTower *tower = it->second;
0541     assert(tower);
0542     // cout << "HCALIN SIM: ";
0543     // tower->identify();
0544 
0545     const int col = tower->get_bineta();
0546     const int row = tower->get_binphi();
0547     const int chanNum = getChannelNumber(row,col);
0548 
0549     if (col < 0 or col >= 4)
0550       continue;
0551     if (row < 0 or row >= 4)
0552       continue;
0553 
0554     const double energy_sim = tower->get_energy()/samplefrac_in; // apply samplefrac_in
0555     hcalin_sum_e_sim += energy_sim;
0556     _tower.hcalin_twr_sim[chanNum] = energy_sim;
0557 
0558     string HistName_sim = Form("h_hcalin_tower_%d_sim",chanNum);
0559     TH1F *h_hcalin_sim = dynamic_cast<TH1F *>(hm->getHisto(HistName_sim.c_str()));
0560     assert(h_hcalin_sim);
0561     h_hcalin_sim->Fill(_tower.hcalin_twr_sim[chanNum]);
0562       }
0563       _tower.hcalin_e_sim = hcalin_sum_e_sim;
0564 
0565       auto range_hcalout_sim = TOWER_SIM_HCALOUT->getTowers(); // sim 
0566       for (auto it = range_hcalout_sim.first; it != range_hcalout_sim.second; ++it)
0567       {
0568     RawTower *tower = it->second;
0569     assert(tower);
0570     // cout << "HCALOUT SIM: ";
0571     // tower->identify();
0572 
0573     const int col = tower->get_bineta();
0574     const int row = tower->get_binphi();
0575     const int chanNum = getChannelNumber(row,col);
0576 
0577     if (col < 0 or col >= 4)
0578       continue;
0579     if (row < 0 or row >= 4)
0580       continue;
0581 
0582     const double energy_sim = tower->get_energy()/samplefrac_out; // apply samplefrac_out
0583     hcalout_sum_e_sim += energy_sim;
0584     _tower.hcalout_twr_sim[chanNum] = energy_sim;
0585 
0586     string HistName_sim = Form("h_hcalout_tower_%d_sim",chanNum);
0587     TH1F *h_hcalout_sim = dynamic_cast<TH1F *>(hm->getHisto(HistName_sim.c_str()));
0588     assert(h_hcalout_sim);
0589     h_hcalout_sim->Fill(_tower.hcalout_twr_sim[chanNum]);
0590       }
0591       _tower.hcalout_e_sim = hcalout_sum_e_sim;
0592     }
0593     _tower.hcal_total_sim = hcalin_sum_e_sim + hcalout_sum_e_sim;
0594     _tower.hcal_asym_sim = (hcalin_sum_e_sim - hcalout_sum_e_sim)/(hcalin_sum_e_sim + hcalout_sum_e_sim);
0595 
0596     TH2F *h_tower_energy_sim = dynamic_cast<TH2F *>(hm->getHisto("h_tower_energy_sim"));
0597     assert(h_tower_energy_sim);
0598     h_tower_energy_sim->Fill(_tower.hcal_asym_sim,_tower.hcal_total_sim);
0599   }
0600 
0601   // process HCALIN LG
0602   double hcalin_sum_lg_e_raw = 0;
0603   double hcalin_sum_lg_e_calib = 0;
0604   {
0605     auto range_hcalin_lg_raw = TOWER_RAW_LG_HCALIN->getTowers(); // raw
0606     for (auto it = range_hcalin_lg_raw.first; it != range_hcalin_lg_raw.second; ++it)
0607     {
0608       RawTower *tower = it->second;
0609       assert(tower);
0610       // cout << "HCALIN RAW: ";
0611       // tower->identify();
0612 
0613       const int col = tower->get_bineta();
0614       const int row = tower->get_binphi();
0615       const int chanNum = getChannelNumber(row,col);
0616 
0617       if (col < 0 or col >= 4)
0618         continue;
0619       if (row < 0 or row >= 4)
0620         continue;
0621 
0622       // const double energy_raw = tower->get_energy()*towercalib_lg_in[chanNum]; // manual calibration
0623       const double energy_raw = tower->get_energy(); // raw
0624       hcalin_sum_lg_e_raw += energy_raw;
0625       _tower.hcalin_lg_twr_raw[chanNum] = energy_raw;
0626 
0627       string HistName_raw = Form("h_hcalin_lg_tower_%d_raw",chanNum);
0628       TH1F *h_hcalin_lg_raw = dynamic_cast<TH1F *>(hm->getHisto(HistName_raw.c_str()));
0629       assert(h_hcalin_lg_raw);
0630       h_hcalin_lg_raw->Fill(_tower.hcalin_lg_twr_raw[chanNum]);
0631     }
0632     _tower.hcalin_lg_e_raw = hcalin_sum_lg_e_raw;
0633 
0634     auto range_hcalin_lg_calib = TOWER_CALIB_LG_HCALIN->getTowers(); // calib
0635     for (auto it = range_hcalin_lg_calib.first; it != range_hcalin_lg_calib.second; ++it)
0636     {
0637       RawTower *tower = it->second;
0638       assert(tower);
0639 
0640       const int col = tower->get_bineta();
0641       const int row = tower->get_binphi();
0642       const int chanNum = getChannelNumber(row,col);
0643 
0644       if (col < 0 or col >= 4)
0645         continue;
0646       if (row < 0 or row >= 4)
0647         continue;
0648 
0649       const double energy_calib = tower->get_energy();
0650       hcalin_sum_lg_e_calib += energy_calib;
0651       _tower.hcalin_lg_twr_calib[chanNum] = energy_calib;
0652 
0653       string HistName_calib = Form("h_hcalin_lg_tower_%d_calib",chanNum);
0654       TH1F *h_hcalin_lg_calib = dynamic_cast<TH1F *>(hm->getHisto(HistName_calib.c_str()));
0655       assert(h_hcalin_lg_calib);
0656       h_hcalin_lg_calib->Fill(_tower.hcalin_lg_twr_calib[chanNum]);
0657     }
0658     _tower.hcalin_lg_e_calib = hcalin_sum_lg_e_calib;
0659   }
0660 
0661   // process HCALOUT LG
0662   double hcalout_sum_lg_e_raw = 0;
0663   double hcalout_sum_lg_e_calib = 0;
0664   {
0665     auto range_hcalout_lg_raw = TOWER_RAW_LG_HCALOUT->getTowers(); // raw
0666     for (auto it = range_hcalout_lg_raw.first; it != range_hcalout_lg_raw.second; ++it)
0667     {
0668       RawTower *tower = it->second;
0669       assert(tower);
0670       // cout << "HCALOUT RAM: ";
0671       // tower->identify();
0672 
0673       const int col = tower->get_bineta();
0674       const int row = tower->get_binphi();
0675       const int chanNum = getChannelNumber(row,col);
0676 
0677       if (col < 0 or col >= 4)
0678         continue;
0679       if (row < 0 or row >= 4)
0680         continue;
0681 
0682       // const double energy_raw = tower->get_energy()*towercalib_lg_out[chanNum]; // manual calibration
0683       const double energy_raw = tower->get_energy(); // raw
0684       hcalout_sum_lg_e_raw += energy_raw;
0685       _tower.hcalout_lg_twr_raw[chanNum] = energy_raw;
0686 
0687       string HistName_raw = Form("h_hcalout_lg_tower_%d_raw",chanNum);
0688       TH1F *h_hcalout_lg_raw = dynamic_cast<TH1F *>(hm->getHisto(HistName_raw.c_str()));
0689       assert(h_hcalout_lg_raw);
0690       h_hcalout_lg_raw->Fill(_tower.hcalout_lg_twr_raw[chanNum]);
0691     }
0692     _tower.hcalout_lg_e_raw = hcalout_sum_lg_e_raw;
0693 
0694     auto range_hcalout_lg_calib = TOWER_CALIB_LG_HCALOUT->getTowers(); // calib
0695     for (auto it = range_hcalout_lg_calib.first; it != range_hcalout_lg_calib.second; ++it)
0696     {
0697       RawTower *tower = it->second;
0698       assert(tower);
0699 
0700       const int col = tower->get_bineta();
0701       const int row = tower->get_binphi();
0702       const int chanNum = getChannelNumber(row,col);
0703 
0704       if (col < 0 or col >= 4)
0705         continue;
0706       if (row < 0 or row >= 4)
0707         continue;
0708 
0709       const double energy_calib = tower->get_energy();
0710       hcalout_sum_lg_e_calib += energy_calib;
0711       _tower.hcalout_lg_twr_calib[chanNum] = energy_calib;
0712 
0713       string HistName_calib = Form("h_hcalout_lg_tower_%d_calib",chanNum);
0714       TH1F *h_hcalout_lg_calib = dynamic_cast<TH1F *>(hm->getHisto(HistName_calib.c_str()));
0715       assert(h_hcalout_lg_calib);
0716       h_hcalout_lg_calib->Fill(_tower.hcalout_lg_twr_calib[chanNum]);
0717     }
0718     _tower.hcalout_lg_e_calib = hcalout_sum_lg_e_calib;
0719   }
0720 
0721   _tower.hcal_total_raw = hcalin_sum_lg_e_raw + hcalout_sum_lg_e_raw;
0722   _tower.hcal_asym_raw = (hcalin_sum_lg_e_raw - hcalout_sum_lg_e_raw)/(hcalin_sum_lg_e_raw + hcalout_sum_lg_e_raw);
0723   TH2F *h_tower_energy_raw = dynamic_cast<TH2F *>(hm->getHisto("h_tower_energy_raw"));
0724   assert(h_tower_energy_raw);
0725   h_tower_energy_raw->Fill(_tower.hcal_asym_raw,_tower.hcal_total_raw); 
0726 
0727   _tower.hcal_total_calib = hcalin_sum_lg_e_calib + hcalout_sum_lg_e_calib;
0728   _tower.hcal_asym_calib = (hcalin_sum_lg_e_calib - hcalout_sum_lg_e_calib)/(hcalin_sum_lg_e_calib + hcalout_sum_lg_e_calib);
0729   TH2F *h_tower_energy_calib = dynamic_cast<TH2F *>(hm->getHisto("h_tower_energy_calib"));
0730   assert(h_tower_energy_calib);
0731   h_tower_energy_calib->Fill(_tower.hcal_asym_calib,_tower.hcal_total_calib); 
0732 
0733   // process HCALOUT HG
0734   {
0735     double hcalout_sum_hg_e_raw = 0;
0736     double hcalout_sum_hg_e_calib = 0;
0737     auto range_hcalout_hg_raw = TOWER_RAW_HG_HCALOUT->getTowers(); // raw
0738     for (auto it = range_hcalout_hg_raw.first; it != range_hcalout_hg_raw.second; ++it)
0739     {
0740       RawTower *tower = it->second;
0741       assert(tower);
0742 
0743       const int col = tower->get_bineta();
0744       const int row = tower->get_binphi();
0745       const int chanNum = getChannelNumber(row,col);
0746 
0747       if (col < 0 or col >= 4)
0748         continue;
0749       if (row < 0 or row >= 4)
0750         continue;
0751 
0752       // const double energy_raw = tower->get_energy()*towercalib_hg_out[chanNum]; // manual calibration
0753       const double energy_raw = tower->get_energy(); // raw
0754       hcalout_sum_hg_e_raw += energy_raw;
0755       _tower.hcalout_hg_twr_raw[chanNum] = energy_raw;
0756 
0757       string HistName_raw = Form("h_hcalout_hg_tower_%d_raw",chanNum);
0758       TH1F *h_hcalout_hg_raw = dynamic_cast<TH1F *>(hm->getHisto(HistName_raw.c_str()));
0759       assert(h_hcalout_hg_raw);
0760       h_hcalout_hg_raw->Fill(_tower.hcalout_hg_twr_raw[chanNum]);
0761     }
0762     _tower.hcalout_hg_e_raw = hcalout_sum_hg_e_raw;
0763 
0764     auto range_hcalout_hg_calib = TOWER_CALIB_HG_HCALOUT->getTowers(); // calib
0765     for (auto it = range_hcalout_hg_calib.first; it != range_hcalout_hg_calib.second; ++it)
0766     {
0767       RawTower *tower = it->second;
0768       assert(tower);
0769 
0770       const int col = tower->get_bineta();
0771       const int row = tower->get_binphi();
0772       const int chanNum = getChannelNumber(row,col);
0773 
0774       if (col < 0 or col >= 4)
0775         continue;
0776       if (row < 0 or row >= 4)
0777         continue;
0778 
0779       const double energy_calib = tower->get_energy();
0780       hcalout_sum_hg_e_calib += energy_calib;
0781       _tower.hcalout_hg_twr_calib[chanNum] = energy_calib;
0782 
0783       string HistName_calib = Form("h_hcalout_hg_tower_%d_calib",chanNum);
0784       TH1F *h_hcalout_hg_calib = dynamic_cast<TH1F *>(hm->getHisto(HistName_calib.c_str()));
0785       assert(h_hcalout_hg_calib);
0786       h_hcalout_hg_calib->Fill(_tower.hcalout_hg_twr_calib[chanNum]);
0787     }
0788     _tower.hcalout_hg_e_calib = hcalout_sum_hg_e_calib;
0789   }
0790 
0791   _eval_run.sum_E_HCAL_IN = hcalin_sum_lg_e_calib;
0792   _eval_run.sum_E_HCAL_OUT = hcalout_sum_lg_e_calib;
0793 
0794   _eval_run.good_e = good_e;
0795   _eval_run.good_anti_e = good_anti_e;
0796 
0797   TTree *T = dynamic_cast<TTree *>(hm->getHisto("HCAL_Info"));
0798   assert(T);
0799   T->Fill();
0800 
0801   return Fun4AllReturnCodes::EVENT_OK;
0802 }
0803 
0804 pair<int, int>
0805 Proto4ShowerCalib::find_max(RawTowerContainer *towers, int cluster_size)
0806 {
0807   const int clus_edge_size = (cluster_size - 1) / 2;
0808   assert(clus_edge_size >= 0);
0809 
0810   pair<int, int> max(-1000, -1000);
0811   double max_e = 0;
0812 
0813   for (int col = 0; col < n_size; ++col)
0814     for (int row = 0; row < n_size; ++row)
0815     {
0816       double energy = 0;
0817 
0818       for (int dcol = col - clus_edge_size; dcol <= col + clus_edge_size;
0819            ++dcol)
0820         for (int drow = row - clus_edge_size; drow <= row + clus_edge_size;
0821              ++drow)
0822         {
0823           if (dcol < 0 or drow < 0)
0824             continue;
0825 
0826           RawTower *t = towers->getTower(dcol, drow);
0827           if (t)
0828             energy += t->get_energy();
0829         }
0830 
0831       if (energy > max_e)
0832       {
0833         max = make_pair(col, row);
0834         max_e = energy;
0835       }
0836     }
0837 
0838   return max;
0839 }
0840 
0841 int Proto4ShowerCalib::InitAna()
0842 {
0843   if(_is_sim) _mList = "SIM";
0844   if(!_is_sim) _mList = "RAW";
0845   string inputdir = "/sphenix/user/xusun/TestBeam/ShowerCalib/";
0846   string InPutList = Form("/direct/phenix+u/xusun/WorkSpace/sPHENIX/analysis/Prototype4/HCal/macros/list/ShowerCalib/Proto4TowerInfo%s_%s.list",_mList.c_str(),_mRunID.c_str());
0847 
0848   mChainInPut = new TChain("HCAL_Info");
0849 
0850   if (!InPutList.empty())   // if input file is ok
0851   { 
0852     cout << "Open input database file list: " << InPutList.c_str() << endl;
0853     ifstream in(InPutList.c_str());  // input stream
0854     if(in)
0855     { 
0856       cout << "input database file list is ok" << endl;
0857       char str[255];       // char array for each file name
0858       Long64_t entries_save = 0;
0859       while(in)
0860       { 
0861     in.getline(str,255);  // take the lines of the file list
0862     if(str[0] != 0)
0863     { 
0864       string addfile;
0865       addfile = str;
0866       addfile = inputdir+addfile;
0867       mChainInPut->AddFile(addfile.c_str(),-1,"HCAL_Info");
0868       long file_entries = mChainInPut->GetEntries();
0869       cout << "File added to data chain: " << addfile.c_str() << " with " << (file_entries-entries_save) << " entries" << endl;
0870       entries_save = file_entries;
0871     }
0872       }
0873     }
0874     else
0875     { 
0876       cout << "WARNING: input database file input is problemtic" << endl;
0877       _mInPut_flag = 0;
0878     }
0879   }
0880 
0881   // Set the input tree
0882   if (_mInPut_flag == 1 && !mChainInPut->GetBranch( "info" ))
0883   {
0884     cerr << "ERROR: Could not find branch 'info' in tree!" << endl;
0885   }
0886 
0887   _mInfo = new Proto4ShowerCalib::Eval_Run();
0888   _mTower = new Proto4ShowerCalib::HCAL_Tower();
0889   if(_mInPut_flag == 1)
0890   {
0891     mChainInPut->SetBranchAddress("info", &_mInfo);
0892     mChainInPut->SetBranchAddress("tower", &_mTower);
0893 
0894     long NumOfEvents = (long)mChainInPut->GetEntries();
0895     cout << "total number of events: " << NumOfEvents << endl;
0896     _mStartEvent = 0;
0897     _mStopEvent = NumOfEvents;
0898 
0899     cout << "New nStartEvent = " << _mStartEvent << ", new nStopEvent = " << _mStopEvent << endl;
0900   }
0901 
0902   float histo_range = find_range();
0903 
0904   if(_is_sim)
0905   {
0906     // h_mAsymmEnergy_mixed_sim_wo_cut = new TH2F("h_mAsymmEnergy_mixed_sim_wo_cut","h_mAsymmEnergy_mixed_sim_wo_cut",100,-1.0,1.0,100,-0.5,histo_range);
0907     h_mAsymmEnergy_pion_sim_wo_cut = new TH2F("h_mAsymmEnergy_pion_sim_wo_cut","h_mAsymmEnergy_pion_sim_wo_cut",100,-1.0,1.0,100,-0.5,histo_range);
0908 
0909     // h_mAsymmEnergy_mixed_sim = new TH2F("h_mAsymmEnergy_mixed_sim","h_mAsymmEnergy_mixed_sim",100,-1.0,1.0,100,-0.5,histo_range);
0910     h_mAsymmEnergy_pion_sim = new TH2F("h_mAsymmEnergy_pion_sim","h_mAsymmEnergy_pion_sim",100,-1.0,1.0,100,-0.5,histo_range);
0911   }
0912   if(!_is_sim)
0913   {
0914     h_mAsymmAdc_mixed = new TH2F("h_mAsymmAdc_mixed","h_mAsymmAdc_mixed",110,-1.1,1.1,100,-0.5,histo_range*640.0);
0915     h_mAsymmAdc_electron = new TH2F("h_mAsymmAdc_electron","h_mAsymmAdc_electron",110,-1.1,1.1,100,-0.5,histo_range*640.0);
0916     h_mAsymmAdc_pion = new TH2F("h_mAsymmAdc_pion","h_mAsymmAdc_pion",110,-1.1,1.1,100,-0.5,histo_range*640.0);
0917 
0918     h_mAsymmEnergy_mixed_wo_cut = new TH2F("h_mAsymmEnergy_mixed_wo_cut","h_mAsymmEnergy_mixed_wo_cut",110,-1.1,1.1,100,-0.5,0.5*histo_range);
0919     h_mAsymmEnergy_electron_wo_cut = new TH2F("h_mAsymmEnergy_electron_wo_cut","h_mAsymmEnergy_electron_wo_cut",110,-1.1,1.1,100,-0.5,0.5*histo_range);
0920     h_mAsymmEnergy_pion_wo_cut = new TH2F("h_mAsymmEnergy_pion_wo_cut","h_mAsymmEnergy_pion_wo_cut",110,-1.1,1.1,100,-0.5,0.5*histo_range);
0921 
0922     h_mAsymmEnergy_mixed = new TH2F("h_mAsymmEnergy_mixed","h_mAsymmEnergy_mixed",110,-1.1,1.1,100,-0.5,0.5*histo_range);
0923     h_mAsymmEnergy_electron = new TH2F("h_mAsymmEnergy_electron","h_mAsymmEnergy_electron",110,-1.1,1.1,100,-0.5,0.5*histo_range);
0924     h_mAsymmEnergy_pion = new TH2F("h_mAsymmEnergy_pion","h_mAsymmEnergy_pion",110,-1.1,1.1,100,-0.5,0.5*histo_range);
0925 
0926     h_mAsymmEnergy_mixed_balancing = new TH2F("h_mAsymmEnergy_mixed_balancing","h_mAsymmEnergy_mixed_balancing",110,-1.1,1.1,100,-0.5,0.5*histo_range);
0927     h_mAsymmEnergy_electron_balancing = new TH2F("h_mAsymmEnergy_electron_balancing","h_mAsymmEnergy_electron_balancing",110,-1.1,1.1,100,-0.5,0.5*histo_range);
0928     h_mAsymmEnergy_pion_balancing = new TH2F("h_mAsymmEnergy_pion_balancing","h_mAsymmEnergy_pion_balancing",110,-1.1,1.1,100,-0.5,0.5*histo_range);
0929 
0930     h_mAsymmEnergy_mixed_leveling = new TH2F("h_mAsymmEnergy_mixed_leveling","h_mAsymmEnergy_mixed_leveling",110,-1.1,1.1,100,-0.5,0.5*histo_range);
0931     h_mAsymmEnergy_electron_leveling = new TH2F("h_mAsymmEnergy_electron_leveling","h_mAsymmEnergy_electron_leveling",110,-1.1,1.1,100,-0.5,0.5*histo_range);
0932     h_mAsymmEnergy_pion_leveling = new TH2F("h_mAsymmEnergy_pion_leveling","h_mAsymmEnergy_pion_leveling",110,-1.1,1.1,100,-0.5,0.5*histo_range);
0933 
0934     h_mAsymmEnergy_mixed_ShowerCalib = new TH2F("h_mAsymmEnergy_mixed_ShowerCalib","h_mAsymmEnergy_mixed_ShowerCalib",110,-1.1,1.1,100,-0.5,histo_range+15.0);
0935     h_mAsymmEnergy_electron_ShowerCalib = new TH2F("h_mAsymmEnergy_electron_ShowerCalib","h_mAsymmEnergy_electron_ShowerCalib",110,-1.1,1.1,100,-0.5,histo_range+15.0);
0936     h_mAsymmEnergy_pion_ShowerCalib = new TH2F("h_mAsymmEnergy_pion_ShowerCalib","h_mAsymmEnergy_pion_ShowerCalib",110,-1.1,1.1,100,-0.5,histo_range+15.0);
0937 
0938     // Outer HCal only study
0939     h_mAsymmEnergy_mixed_MIP = new TH2F("h_mAsymmEnergy_mixed_MIP","h_mAsymmEnergy_mixed_MIP",110,-1.1,1.1,100,-0.5,0.5*histo_range);
0940     h_mEnergyOut_electron = new TH1F("h_mEnergyOut_electron","h_mEnergyOut_electron",100,-0.5,0.5*histo_range);
0941     h_mEnergyOut_pion = new TH1F("h_mEnergyOut_pion","h_mEnergyOut_pion",100,-0.5,0.5*histo_range);
0942 
0943     h_mEnergyOut_electron_ShowerCalib = new TH1F("h_mEnergyOut_electron_ShowerCalib","h_mEnergyOut_electron_ShowerCalib",100,-0.5,histo_range);
0944     h_mEnergyOut_pion_ShowerCalib = new TH1F("h_mEnergyOut_pion_ShowerCalib","h_mEnergyOut_pion_ShowerCalib",100,-0.5,histo_range);
0945   }
0946 
0947   return 0;
0948 }
0949 
0950 int Proto4ShowerCalib::MakeAna()
0951 {
0952   cout << "Make()" << endl;
0953 
0954   // const int mBeamEnergy[6] = {4, 8, 12, 16, 24, 30};
0955   const float c_in[5]  = {0.905856, 0.943406, 0.964591, 1.0, 0.962649};
0956   const float c_out[5] = {1.115980, 1.063820, 1.038110, 1.0, 1.040370};
0957   // const float c_in[5]  = {0.942372, 0.932490, 0.955448, 0.962649, 0.962649};
0958   // const float c_out[5] = {1.065140, 1.078050, 1.048910, 1.040370, 1.040370};
0959   const int mEnergyBin = find_energy();
0960   const float c_in_leveling = c_in[mEnergyBin];
0961   const float c_out_leveling = c_out[mEnergyBin];
0962   // cout << "mEnergyBin = " << mEnergyBin << ", c_in = " << c_in[mEnergyBin] << ", c_out = " << c_out[mEnergyBin] << endl;
0963 
0964   unsigned long start_event_use = _mStartEvent;
0965   unsigned long stop_event_use = _mStopEvent;
0966 
0967   mChainInPut->SetBranchAddress("info", &_mInfo);
0968   mChainInPut->SetBranchAddress("tower", &_mTower);
0969   mChainInPut->GetEntry(0);
0970 
0971   for(unsigned long i_event = start_event_use; i_event < stop_event_use; ++i_event)
0972   // for(unsigned long i_event = 20; i_event < 40; ++i_event)
0973   {
0974     if (!mChainInPut->GetEntry( i_event )) // take the event -> information is stored in event
0975       break;  // end of data chunk
0976 
0977     if (i_event != 0  &&  i_event % 1000 == 0)
0978       cout << "." << flush;
0979     if (i_event != 0  &&  i_event % 10000 == 0)
0980     {
0981       if((stop_event_use-start_event_use) > 0)
0982       {
0983     double event_percent = 100.0*((double)(i_event-start_event_use))/((double)(stop_event_use-start_event_use));
0984     cout << " " << i_event-start_event_use << " (" << event_percent << "%) " << "\n" << "==> Processing data (ShowerCalib) " << flush;
0985       }
0986     }
0987 
0988     if(_is_sim) // beam test simulation
0989     {
0990       float energy_sim = _mTower->hcal_total_sim;
0991       float asymm_sim = _mTower->hcal_asym_sim;
0992       const float energy_hcalin_sim = _mTower->hcalin_e_sim; // sim
0993       const float energy_hcalout_sim = _mTower->hcalout_e_sim;
0994 
0995       // h_mAsymmEnergy_mixed_sim_wo_cut->Fill(asymm_sim,energy_sim);
0996       h_mAsymmEnergy_pion_sim_wo_cut->Fill(asymm_sim,energy_sim);
0997       if(energy_hcalin_sim > 0.2 && energy_hcalout_sim > 0/2)
0998       {
0999     // h_mAsymmEnergy_mixed_sim->Fill(asymm_sim,energy_sim);
1000     h_mAsymmEnergy_pion_sim->Fill(asymm_sim,energy_sim);
1001       }
1002     }
1003     if(!_is_sim) // beam test data
1004     {
1005       const bool good_electron = _mInfo->good_e;
1006       const bool good_pion = _mInfo->good_anti_e;
1007 
1008       float energy_raw = _mTower->hcal_total_raw; // raw
1009       float asymm_raw = _mTower->hcal_asym_raw;
1010 
1011       h_mAsymmAdc_mixed->Fill(asymm_raw,energy_raw);
1012       if(good_electron) h_mAsymmAdc_electron->Fill(asymm_raw,energy_raw);
1013       if(good_pion) h_mAsymmAdc_pion->Fill(asymm_raw,energy_raw);
1014 
1015       float energy_calib = _mTower->hcal_total_calib; // calib
1016       float asymm_calib = _mTower->hcal_asym_calib;
1017       const float energy_hcalin_calib = _mTower->hcalin_lg_e_calib;
1018       const float energy_hcalout_calib = _mTower->hcalout_lg_e_calib;
1019 
1020       h_mAsymmEnergy_mixed_wo_cut->Fill(asymm_calib,energy_calib);
1021       if(good_electron) h_mAsymmEnergy_electron_wo_cut->Fill(asymm_calib,energy_calib);
1022       if(good_pion) h_mAsymmEnergy_pion_wo_cut->Fill(asymm_calib,energy_calib);
1023 
1024       if(energy_hcalin_calib > 0.001 && energy_hcalout_calib > 0.001) // remove ped
1025       { // extract MIP
1026     h_mAsymmEnergy_mixed->Fill(asymm_calib,energy_calib);
1027     if(good_electron) h_mAsymmEnergy_electron->Fill(asymm_calib,energy_calib);
1028     if(good_pion) h_mAsymmEnergy_pion->Fill(asymm_calib,energy_calib);
1029       }
1030 
1031       const double MIP_cut = MIP_mean+MIP_width; // 1 sigma MIP energy cut
1032       if(energy_hcalin_calib <= MIP_cut && energy_hcalout_calib > 0.001 && energy_calib > MIP_cut)
1033       { // OHCal with MIP event in IHCal
1034     h_mAsymmEnergy_mixed_MIP->Fill(asymm_calib,energy_calib);
1035     if(good_electron) h_mEnergyOut_electron->Fill(energy_hcalout_calib);
1036     if(good_pion) h_mEnergyOut_pion->Fill(energy_hcalout_calib);
1037 
1038     // showercalib
1039     if(good_electron) h_mEnergyOut_electron_ShowerCalib->Fill(energy_hcalout_calib*showercalib_ohcal);
1040     if(good_pion) h_mEnergyOut_pion_ShowerCalib->Fill(energy_hcalout_calib*showercalib_ohcal);
1041       }
1042       const double MIP_energy_cut = MIP_mean+3.0*MIP_width; // 3 sigma MIP energy cut
1043       if(energy_hcalin_calib > 0.001 && energy_hcalout_calib > 0.001 && energy_calib > MIP_energy_cut)
1044       { // balancing without muon
1045     h_mAsymmEnergy_mixed_balancing->Fill(asymm_calib,energy_calib);
1046     if(good_electron) h_mAsymmEnergy_electron_balancing->Fill(asymm_calib,energy_calib);
1047     if(good_pion) h_mAsymmEnergy_pion_balancing->Fill(asymm_calib,energy_calib);
1048 
1049     // apply leveling
1050     const float energy_leveling = c_in_leveling*energy_hcalin_calib + c_out_leveling*energy_hcalout_calib;
1051     const float asymm_leveling = (c_in_leveling*energy_hcalin_calib - c_out_leveling*energy_hcalout_calib)/energy_leveling;
1052     h_mAsymmEnergy_mixed_leveling->Fill(asymm_leveling,energy_leveling);
1053     if(good_electron) h_mAsymmEnergy_electron_leveling->Fill(asymm_leveling,energy_leveling);
1054     if(good_pion) h_mAsymmEnergy_pion_leveling->Fill(asymm_leveling,energy_leveling);
1055 
1056     // apply shower calibration
1057     const float energy_showercalib = showercalib*c_in_leveling*energy_hcalin_calib + showercalib*c_out_leveling*energy_hcalout_calib;
1058     const float asymm_showercalib = (showercalib*c_in_leveling*energy_hcalin_calib - showercalib*c_out_leveling*energy_hcalout_calib)/energy_showercalib;
1059     // cout << "energy_leveling = " << energy_leveling << ", energy_showercalib = " << energy_showercalib << endl;
1060     // cout << "asymm_leveling = " << asymm_leveling << ", asymm_showercalib = " << asymm_showercalib << endl;
1061     h_mAsymmEnergy_mixed_ShowerCalib->Fill(asymm_showercalib,energy_showercalib);
1062     if(good_electron) h_mAsymmEnergy_electron_ShowerCalib->Fill(asymm_showercalib,energy_showercalib);
1063     if(good_pion) h_mAsymmEnergy_pion_ShowerCalib->Fill(asymm_showercalib,energy_showercalib);
1064       }
1065     }
1066   }
1067 
1068   cout << "." << flush;
1069   cout << " " << stop_event_use-start_event_use << "(" << 100 << "%)";
1070   cout << endl;
1071 
1072   return 1;
1073 }
1074 
1075 int Proto4ShowerCalib::FinishAna()
1076 {
1077   cout << "Finish()" << endl;
1078   mFile_OutPut = new TFile(_filename.c_str(),"RECREATE");
1079   mFile_OutPut->cd();
1080   if(_is_sim) // beam test simulation
1081   {
1082     // h_mAsymmEnergy_mixed_sim_wo_cut->Write();
1083     h_mAsymmEnergy_pion_sim_wo_cut->Write();
1084     // h_mAsymmEnergy_mixed_sim->Write();
1085     h_mAsymmEnergy_pion_sim->Write();
1086   }
1087   if(!_is_sim) // beam test data
1088   {
1089     h_mAsymmAdc_mixed->Write();
1090     h_mAsymmAdc_electron->Write();
1091     h_mAsymmAdc_pion->Write();
1092 
1093     h_mAsymmEnergy_mixed_wo_cut->Write();
1094     h_mAsymmEnergy_electron_wo_cut->Write();
1095     h_mAsymmEnergy_pion_wo_cut->Write();
1096 
1097     h_mAsymmEnergy_mixed->Write();
1098     h_mAsymmEnergy_electron->Write();
1099     h_mAsymmEnergy_pion->Write();
1100 
1101     h_mAsymmEnergy_mixed_balancing->Write();
1102     h_mAsymmEnergy_electron_balancing->Write();
1103     h_mAsymmEnergy_pion_balancing->Write();
1104 
1105     h_mAsymmEnergy_mixed_leveling->Write();
1106     h_mAsymmEnergy_electron_leveling->Write();
1107     h_mAsymmEnergy_pion_leveling->Write();
1108 
1109     h_mAsymmEnergy_mixed_ShowerCalib->Write();
1110     h_mAsymmEnergy_electron_ShowerCalib->Write();
1111     h_mAsymmEnergy_pion_ShowerCalib->Write();
1112 
1113     h_mAsymmEnergy_mixed_MIP->Write();
1114     h_mEnergyOut_electron->Write();
1115     h_mEnergyOut_pion->Write();
1116 
1117     h_mEnergyOut_electron_ShowerCalib->Write();
1118     h_mEnergyOut_pion_ShowerCalib->Write();
1119   }
1120   mFile_OutPut->Close();
1121 
1122   return 1;
1123 }
1124 
1125 int Proto4ShowerCalib::getChannelNumber(int row, int column)
1126 {
1127   if(!_is_sim)
1128   {
1129     int hbdchanIHC[4][4] = {{4, 8, 12, 16},
1130                             {3, 7, 11, 15},
1131                             {2, 6, 10, 14},
1132                             {1, 5,  9, 13}};
1133 
1134     return hbdchanIHC[row][column] - 1;
1135   }
1136 
1137   if(_is_sim)
1138   {
1139     int hbdchanIHC[4][4] = {{13,  9, 5, 1},
1140                             {14, 10, 6, 2},
1141                             {15, 11, 7, 3},
1142                             {16, 12, 8, 4}};
1143 
1144     return hbdchanIHC[row][column] - 1;
1145   }
1146 
1147   return -1;
1148 }
1149 
1150 int Proto4ShowerCalib::setTowerCalibParas()
1151 {
1152   const double energy_in[16] = {0.00763174, 0.00692298, 0.00637355, 0.0059323, 0.00762296, 0.00691832, 0.00636611, 0.0059203, 0.00762873, 0.00693594, 0.00637791, 0.00592433, 0.00762898, 0.00691679, 0.00636373, 0.00592433};
1153   const double adc_in[16] = {2972.61, 2856.43, 2658.19, 2376.10, 3283.39, 2632.81, 2775.77, 2491.68, 2994.11, 3385.70, 3258.01, 2638.31, 3479.97, 3081.41, 2768.36, 2626.77};
1154   const double gain_factor_in = 32.0;
1155 
1156   const double energy_out[16] = {0.00668176, 0.00678014, 0.00687082, 0.00706854, 0.00668973, 0.00678279, 0.00684794, 0.00705448, 0.00668976, 0.0068013, 0.00685931, 0.00704985, 0.0066926, 0.00678282, 0.00684403, 0.00704143};
1157   // const double adc_out[16] = {666.378, 1488.15, 1493.36, 1816.82, 666.378, 1488.15, 1493.36, 1816.82, 666.378, 1488.15, 1493.36, 1816.82, 666.378, 1488.15, 1493.36, 1816.82}; // use 1st column for whole HCALOUT
1158   const double adc_out[16] = {276.9, 290.0, 280.7, 272.1, 309.5, 304.8, 318.5, 289.6, 289.9, 324.2, 297.9, 294.6, 292.7, 310.5, 302.3, 298.5}; // Songkyo's number
1159   const double adc_amp[16] = {2.505, 5.330, 5.330, 6.965, 2.505, 5.330, 5.330, 6.965, 2.505, 5.330, 5.330, 6.965, 2.505, 5.330, 5.330, 6.965}; // amplify from 2017 to 2018
1160   const double gain_factor_out = 16.0;
1161 
1162   // const double tower_in[16] = {3.98902E-05,3.76295E-05,3.75111E-05,3.8131E-05,3.84078E-05,3.66107E-05,4.20291E-05,3.61503E-05,3.67743E-05,3.22359E-05,3.27909E-05,3.89949E-05,3.55663E-05,3.62068E-05,3.55327E-05,3.34781E-05}; // from Songkyo
1163   // const double tower_out[16] ={0.000144539,0.000136337,0.000129978,0.00013414,0.000135812,0.000120164,0.000123676,0.000119989,0.000133497,0.000128479,0.000116275,0.000128108,0.000131765,0.000126611,0.000121401,0.000126879};
1164 
1165   for(int i_tower = 0; i_tower < 16; ++i_tower)
1166   {
1167     towercalib_lg_in[i_tower] = gain_factor_in*energy_in[i_tower]/(adc_in[i_tower]*samplefrac_in);
1168     towercalib_lg_out[i_tower] = gain_factor_out*energy_out[i_tower]/(adc_amp[i_tower]*adc_out[i_tower]*samplefrac_out);
1169     towercalib_hg_out[i_tower] = energy_out[i_tower]/(adc_amp[i_tower]*adc_out[i_tower]*samplefrac_out);
1170     // towercalib_lg_out[i_tower] = gain_factor_out*energy_out[i_tower]/(adc_out[i_tower]*samplefrac_out);
1171     // towercalib_hg_out[i_tower] = energy_out[i_tower]/(adc_out[i_tower]*samplefrac_out);
1172 
1173     // towercalib_lg_in[i_tower] = tower_in[i_tower]*gain_factor_in;
1174     // towercalib_lg_out[i_tower] = tower_out[i_tower]*gain_factor_out;
1175     // towercalib_hg_out[i_tower] = tower_out[i_tower];
1176 
1177     // cout << "i_tower = " << i_tower << ", towercalib_lg_in = " << towercalib_lg_in[i_tower] << ", towercalib_lg_out = " << towercalib_lg_out[i_tower] << endl;
1178   }
1179 
1180   return 1;
1181 }
1182 
1183 float Proto4ShowerCalib::find_range()
1184 {
1185   float range = 99.5;
1186 
1187   std::string hcal_120GeV[2] = {"0498","0762"};
1188   std::string hcal_60GeV[2] = {"0820","0821"};
1189   std::string hcal_12GeV[2] = {"0570","0571"};
1190   std::string hcal_8GeV[3] = {"0421","0422","1245"};
1191   std::string hcal_5GeV[3] = {"1087","1089","1091"};
1192 
1193   // cout << "_mRunID = " << _mRunID.c_str() << ", compare = " << _mRunID.compare(hcal_4GeV[0]) << endl;
1194 
1195   for(int i_run = 0; i_run < 2; ++i_run) // 120 GeV
1196   {
1197     if(_mRunID.compare(hcal_120GeV[i_run]) == 0) 
1198     {
1199       range = 99.5;
1200       return range;
1201     }
1202   }
1203   for(int i_run = 0; i_run < 2; ++i_run) // 60 GeV
1204   {
1205     if(_mRunID.compare(hcal_60GeV[i_run]) == 0) 
1206     {
1207       range = 79.5;
1208       return range;
1209     }
1210   }
1211   for(int i_run = 0; i_run < 2; ++i_run) // 12 GeV
1212   {
1213     if(_mRunID.compare(hcal_12GeV[i_run]) == 0)
1214     {
1215       range = 24.5;
1216       return range;
1217     }
1218   }
1219   for(int i_run = 0; i_run < 3; ++i_run) // 8 GeV
1220   {
1221     if(_mRunID.compare(hcal_8GeV[i_run]) == 0)
1222     {
1223       range = 16.5;
1224       return range;
1225     }
1226   }
1227   for(int i_run = 0; i_run < 3; ++i_run) // 5 GeV
1228   {
1229     if(_mRunID.compare(hcal_5GeV[i_run]) == 0) 
1230     {
1231       range = 9.5;
1232       return range;
1233     }
1234   }
1235 
1236   return range;
1237 }
1238 
1239 int Proto4ShowerCalib::find_energy()
1240 {
1241   int energy = -1;
1242 
1243   std::string hcal_120GeV[2] = {"0498","0762"};
1244   std::string hcal_60GeV[2] = {"0820","0821"};
1245   std::string hcal_12GeV[2] = {"0570","0571"};
1246   std::string hcal_8GeV[3] = {"0421","0422","1245"};
1247   std::string hcal_5GeV[3] = {"1087","1089","1091"};
1248 
1249   // cout << "_mRunID = " << _mRunID.c_str() << ", compare = " << _mRunID.compare(hcal_4GeV[0]) << endl;
1250 
1251   for(int i_run = 0; i_run < 2; ++i_run) // 120 GeV
1252   {
1253     if(_mRunID.compare(hcal_120GeV[i_run]) == 0) 
1254     {
1255       energy = 4;
1256       return energy;
1257     }
1258   }
1259   for(int i_run = 0; i_run < 2; ++i_run) // 60 GeV
1260   {
1261     if(_mRunID.compare(hcal_60GeV[i_run]) == 0) 
1262     {
1263       energy = 3;
1264       return energy;
1265     }
1266   }
1267   for(int i_run = 0; i_run < 2; ++i_run) // 12 GeV
1268   {
1269     if(_mRunID.compare(hcal_12GeV[i_run]) == 0)
1270     {
1271       energy = 2;
1272       return energy;
1273     }
1274   }
1275   for(int i_run = 0; i_run < 3; ++i_run) // 8 GeV
1276   {
1277     if(_mRunID.compare(hcal_8GeV[i_run]) == 0)
1278     {
1279       energy = 1;
1280       return energy;
1281     }
1282   }
1283   for(int i_run = 0; i_run < 3; ++i_run) // 5 GeV
1284   {
1285     if(_mRunID.compare(hcal_5GeV[i_run]) == 0) 
1286     {
1287       energy = 0;
1288       return energy;
1289     }
1290   }
1291 
1292   return energy;
1293 }
[ 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!