Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:12:16 

0001 #include "EMCalCalib.h"
0002 #include "PhotonPair.h"
0003 
0004 #include <fun4all/SubsysReco.h>
0005 #include <fun4all/Fun4AllServer.h>
0006 #include <fun4all/PHTFileServer.h>
0007 #include <phool/PHCompositeNode.h>
0008 #include <fun4all/Fun4AllReturnCodes.h>
0009 #include <phool/getClass.h>
0010 #include <fun4all/Fun4AllHistoManager.h>
0011 
0012 #include <phool/PHCompositeNode.h>
0013 
0014 #include <g4main/PHG4TruthInfoContainer.h>
0015 #include <g4main/PHG4VtxPoint.h>
0016 #include <g4main/PHG4Particle.h>
0017 
0018 #include <g4hough/SvtxVertexMap.h>
0019 
0020 #include <calobase/RawTowerContainer.h>
0021 #include <calobase/RawTowerGeomContainer.h>
0022 #include <calobase/RawTower.h>
0023 #include <calobase/RawClusterContainer.h>
0024 #include <calobase/RawCluster.h>
0025 
0026 #include <g4eval/CaloEvalStack.h>
0027 #include <g4eval/CaloRawClusterEval.h>
0028 #include <g4eval/CaloRawTowerEval.h>
0029 #include <g4eval/CaloTruthEval.h>
0030 #include <g4eval/SvtxEvalStack.h>
0031 
0032 #include <TNtuple.h>
0033 #include <TFile.h>
0034 #include <TH1F.h>
0035 #include <TH2F.h>
0036 #include <TVector3.h>
0037 #include <TLorentzVector.h>
0038 
0039 #include <exception>
0040 #include <stdexcept>
0041 #include <iostream>
0042 #include <vector>
0043 #include <set>
0044 #include <algorithm>
0045 #include <cassert>
0046 #include <cmath>
0047 
0048 using namespace std;
0049 
0050 EMCalCalib::EMCalCalib(const std::string &filename, EMCalCalib::enu_flags flags) :
0051     SubsysReco("EMCalCalib"), _eval_stack(NULL), _T_EMCalTrk(NULL), _mc_photon(
0052         NULL), //
0053     _magfield(1.5), //
0054     _filename(filename), _flags(flags), _ievent(0)
0055 {
0056 
0057   verbosity = 1;
0058 
0059   _hcalout_hit_container = NULL;
0060   _hcalin_hit_container = NULL;
0061   _cemc_hit_container = NULL;
0062   _hcalout_abs_hit_container = NULL;
0063   _hcalin_abs_hit_container = NULL;
0064   _cemc_abs_hit_container = NULL;
0065   _magnet_hit_container = NULL;
0066   _bh_hit_container = NULL;
0067   _bh_plus_hit_container = NULL;
0068   _bh_minus_hit_container = NULL;
0069   _cemc_electronics_hit_container = NULL;
0070   _hcalin_spt_hit_container = NULL;
0071   _svtx_hit_container = NULL;
0072   _svtx_support_hit_container = NULL;
0073 
0074 }
0075 
0076 EMCalCalib::~EMCalCalib()
0077 {
0078   if (_eval_stack)
0079     {
0080       delete _eval_stack;
0081     }
0082 }
0083 
0084 int
0085 EMCalCalib::InitRun(PHCompositeNode *topNode)
0086 {
0087   _ievent = 0;
0088 
0089   PHNodeIterator iter(topNode);
0090   PHCompositeNode *dstNode = static_cast<PHCompositeNode*>(iter.findFirst(
0091       "PHCompositeNode", "DST"));
0092   if (!dstNode)
0093     {
0094       std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0095       throw runtime_error(
0096           "Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
0097     }
0098 
0099   // get DST objects
0100   _hcalout_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
0101       "G4HIT_HCALOUT");
0102   _hcalin_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
0103       "G4HIT_HCALIN");
0104 
0105   _cemc_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
0106       "G4HIT_CEMC");
0107 
0108   _hcalout_abs_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
0109       "G4HIT_ABSORBER_HCALOUT");
0110 
0111   _hcalin_abs_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
0112       "G4HIT_ABSORBER_HCALIN");
0113 
0114   _cemc_abs_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
0115       "G4HIT_ABSORBER_CEMC");
0116 
0117   if (flag(kProcessUpslisonTrig))
0118     {
0119 //      PhotonPair * photon_pair = findNode::getClass<PhotonPair>(dstNode,
0120 //          "PhotonPair");
0121 //
0122 //      if (not photon_pair)
0123 //        {
0124 //          photon_pair = new PhotonPair();
0125 //
0126 //          PHIODataNode<PHObject> *node = new PHIODataNode<PHObject>(
0127 //              photon_pair, "PhotonPair", "PHObject");
0128 //          dstNode->addNode(node);
0129 //        }
0130 //
0131 //      assert(photon_pair);
0132     }
0133   if (flag(kProcessMCPhoton))
0134     {
0135       _mc_photon = findNode::getClass<MCPhoton>(dstNode, "MCPhoton");
0136 
0137       if (not _mc_photon)
0138         {
0139           _mc_photon = new MCPhoton();
0140 
0141           PHIODataNode<PHObject> *node = new PHIODataNode<PHObject>(_mc_photon,
0142               "MCPhoton", "PHObject");
0143           dstNode->addNode(node);
0144         }
0145 
0146       assert(_mc_photon);
0147     }
0148 
0149   return Fun4AllReturnCodes::EVENT_OK;
0150 }
0151 
0152 int
0153 EMCalCalib::End(PHCompositeNode *topNode)
0154 {
0155   cout << "EMCalCalib::End - write to " << _filename << endl;
0156   PHTFileServer::get().cd(_filename);
0157 
0158   Fun4AllHistoManager *hm = get_HistoManager();
0159   assert(hm);
0160   for (unsigned int i = 0; i < hm->nHistos(); i++)
0161     hm->getHisto(i)->Write();
0162 
0163   if (_T_EMCalTrk)
0164     _T_EMCalTrk->Write();
0165 
0166   // help index files with TChain
0167   TTree * T_Index = new TTree("T_Index", "T_Index");
0168   assert(T_Index);
0169   T_Index->Write();
0170 
0171   return Fun4AllReturnCodes::EVENT_OK;
0172 }
0173 
0174 int
0175 EMCalCalib::Init(PHCompositeNode *topNode)
0176 {
0177 
0178   _ievent = 0;
0179 
0180   cout << "EMCalCalib::get_HistoManager - Making PHTFileServer " << _filename
0181       << endl;
0182   PHTFileServer::get().open(_filename, "RECREATE");
0183 
0184   if (flag(kProcessSF))
0185     {
0186       cout << "EMCalCalib::Init - Process sampling fraction" << endl;
0187       Init_SF(topNode);
0188     }
0189   if (flag(kProcessTower))
0190     {
0191       cout << "EMCalCalib::Init - Process tower occupancies" << endl;
0192       Init_Tower(topNode);
0193     }
0194   if (flag(kProcessMCPhoton))
0195     {
0196       cout << "EMCalCalib::Init - Process trakcs" << endl;
0197       Init_MCPhoton(topNode);
0198     }
0199   if (flag(kProcessUpslisonTrig))
0200     {
0201       cout << "EMCalCalib::Init - Process kProcessUpslisonTrig" << endl;
0202       Init_UpslisonTrig(topNode);
0203     }
0204 
0205   return Fun4AllReturnCodes::EVENT_OK;
0206 }
0207 
0208 int
0209 EMCalCalib::process_event(PHCompositeNode *topNode)
0210 {
0211 
0212   if (verbosity > 2)
0213     cout << "EMCalCalib::process_event() entered" << endl;
0214 
0215   if (_eval_stack)
0216     _eval_stack->next_event(topNode);
0217 
0218   if (flag(kProcessUpslisonTrig))
0219     {
0220       int ret = process_event_UpslisonTrig(topNode);
0221       if (ret != Fun4AllReturnCodes::EVENT_OK)
0222         return ret;
0223     }
0224   if (flag(kProcessSF))
0225     process_event_SF(topNode);
0226   if (flag(kProcessTower))
0227     process_event_Tower(topNode);
0228   if (flag(kProcessMCPhoton))
0229     process_event_MCPhoton(topNode);
0230 
0231   return Fun4AllReturnCodes::EVENT_OK;
0232 }
0233 
0234 int
0235 EMCalCalib::Init_UpslisonTrig(PHCompositeNode *topNode)
0236 {
0237 
0238   return Fun4AllReturnCodes::EVENT_OK;
0239 }
0240 
0241 int
0242 EMCalCalib::process_event_UpslisonTrig(PHCompositeNode *topNode)
0243 {
0244 
0245   if (verbosity > 2)
0246     cout << "EMCalCalib::process_event_UpslisonTrig() entered" << endl;
0247 
0248 //  PhotonPair * photon_pair = findNode::getClass<PhotonPair>(topNode,
0249 //      "PhotonPair");
0250 //  assert(photon_pair);
0251 //  static const double upsilon_mass = 9460.30e-3;
0252 
0253 //  if (!_eval_stack)
0254 //    {
0255 //      _eval_stack = new SvtxEvalStack(topNode);
0256 //      _eval_stack->set_strict(false);
0257 //    }
0258 //
0259 //  SvtxVertexEval* vertexeval = svtxevalstack.get_vertex_eval();
0260 //  SvtxTrackEval* trackeval = svtxevalstack.get_track_eval();
0261 //  SvtxClusterEval* clustereval = svtxevalstack.get_cluster_eval();
0262 //  SvtxHitEval* hiteval = svtxevalstack.get_hit_eval();
0263 //  SvtxTruthEval* trutheval = svtxevalstack.get_truth_eval();
0264 //
0265 //  SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,
0266 //      "SvtxTrackMap");
0267 //  assert(trackmap);
0268 //  PHG4TruthInfoContainer* truthinfo =
0269 //      findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0270 //  assert(truthinfo);
0271 //  PHG4TruthInfoContainer::Map particle_map = truthinfo->GetMap();
0272 //  PHG4TruthInfoContainer::ConstRange primary_range =
0273 //      truthinfo->GetPrimaryParticleRange();
0274 //
0275 //  set<int> e_pos_candidate;
0276 //  set<int> e_neg_candidate;
0277 //
0278 //  for (PHG4TruthInfoContainer::ConstIterator iter = primary_range.first;
0279 //      iter != primary_range.second; ++iter)
0280 //    {
0281 //
0282 //      PHG4Particle * particle = iter->second;
0283 //      assert(particle);
0284 //
0285 //      if (particle->get_pid() == 11)
0286 //        e_neg_candidate.insert(particle->get_track_id());
0287 //      if (particle->get_pid() == -11)
0288 //        e_pos_candidate.insert(particle->get_track_id());
0289 //    }
0290 //
0291 //  map<double, pair<int, int> > mass_diff_id_map;
0292 //  for (set<int>::const_iterator i_e_pos = e_pos_candidate.begin();
0293 //      i_e_pos != e_pos_candidate.end(); ++i_e_pos)
0294 //    for (set<int>::const_iterator i_e_neg = e_neg_candidate.begin();
0295 //        i_e_neg != e_neg_candidate.end(); ++i_e_neg)
0296 //      {
0297 //        TLorentzVector vpos(particle_map[*i_e_pos]->get_px(),
0298 //            particle_map[*i_e_pos]->get_py(), particle_map[*i_e_pos]->get_pz(),
0299 //            particle_map[*i_e_pos]->get_e());
0300 //        TLorentzVector vneg(particle_map[*i_e_neg]->get_px(),
0301 //            particle_map[*i_e_neg]->get_py(), particle_map[*i_e_neg]->get_pz(),
0302 //            particle_map[*i_e_neg]->get_e());
0303 //
0304 //        TLorentzVector invar = vpos + vneg;
0305 //
0306 //        const double mass = invar.M();
0307 //
0308 //        const double mass_diff = fabs(mass - upsilon_mass);
0309 //
0310 //        mass_diff_id_map[mass_diff] = make_pair<int, int>(*i_e_pos, *i_e_neg);
0311 //      }
0312 //
0313 //  if (mass_diff_id_map.size() <= 0)
0314 //    return Fun4AllReturnCodes::DISCARDEVENT;
0315 //  else
0316 //    {
0317 //      const pair<int, int> best_upsilon_pair = mass_diff_id_map.begin()->second;
0318 //
0319 //      //truth mass
0320 //      TLorentzVector gvpos(particle_map[best_upsilon_pair.first]->get_px(),
0321 //          particle_map[best_upsilon_pair.first]->get_py(),
0322 //          particle_map[best_upsilon_pair.first]->get_pz(),
0323 //          particle_map[best_upsilon_pair.first]->get_e());
0324 //      TLorentzVector gvneg(particle_map[best_upsilon_pair.second]->get_px(),
0325 //          particle_map[best_upsilon_pair.second]->get_py(),
0326 //          particle_map[best_upsilon_pair.second]->get_pz(),
0327 //          particle_map[best_upsilon_pair.second]->get_e());
0328 //      TLorentzVector ginvar = gvpos + gvneg;
0329 //      photon_pair->gmass = ginvar.M();
0330 //
0331 //      eval_photon(particle_map[best_upsilon_pair.first], *photon_pair->get_trk(0),
0332 //          topNode);
0333 //      eval_photon(particle_map[best_upsilon_pair.second],
0334 //          *photon_pair->get_trk(1), topNode);
0335 //
0336 //      //calculated mass
0337 //      TLorentzVector vpos(photon_pair->get_trk(0)->px,
0338 //          photon_pair->get_trk(0)->py, photon_pair->get_trk(0)->pz, 0);
0339 //      TLorentzVector vneg(photon_pair->get_trk(1)->px,
0340 //          photon_pair->get_trk(1)->py, photon_pair->get_trk(1)->pz, 0);
0341 //      vpos.SetE(vpos.P());
0342 //      vneg.SetE(vneg.P());
0343 //      TLorentzVector invar = vpos + vneg;
0344 //      photon_pair->mass = invar.M();
0345 //
0346 //      // cuts
0347 //      const bool eta_pos_cut = fabs(gvpos.Eta()) < 1.0;
0348 //      const bool eta_neg_cut = fabs(gvneg.Eta()) < 1.0;
0349 //      const bool reco_upsilon_mass = fabs(photon_pair->mass - upsilon_mass)
0350 //          < 0.2; // two sigma cuts
0351 //      photon_pair->good_upsilon = eta_pos_cut and eta_neg_cut
0352 //          and reco_upsilon_mass;
0353 //
0354 //      if (not photon_pair->good_upsilon)
0355 //        {
0356 //          return Fun4AllReturnCodes::DISCARDEVENT;
0357 //        }
0358 //
0359 //    }
0360 
0361   return Fun4AllReturnCodes::EVENT_OK;
0362 }
0363 
0364 //! project a photon to calorimeters and collect towers around it
0365 void
0366 EMCalCalib::eval_photon(PHG4Particle * g4particle, MCPhoton & mc_photon,
0367     PHCompositeNode *topNode)
0368 {
0369   assert(g4particle);
0370 
0371   if (!_eval_stack)
0372     {
0373       _eval_stack = new SvtxEvalStack(topNode);
0374       _eval_stack->set_strict(false);
0375     }
0376 //  SvtxVertexEval* vertexeval = _eval_stack->get_vertex_eval();
0377   SvtxTrackEval* trackeval = _eval_stack->get_track_eval();
0378 //  SvtxClusterEval* clustereval = _eval_stack->get_cluster_eval();
0379 //  SvtxHitEval* hiteval = _eval_stack->get_hit_eval();
0380   SvtxTruthEval* trutheval = _eval_stack->get_truth_eval();
0381 
0382   PHG4TruthInfoContainer* truthinfo =
0383       findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0384   assert(truthinfo);
0385 
0386   SvtxClusterMap* clustermap = findNode::getClass<SvtxClusterMap>(topNode,
0387       "SvtxClusterMap");
0388   assert(clustermap);
0389 
0390   int gtrackID = g4particle->get_track_id();
0391   int gflavor = g4particle->get_pid();
0392 
0393   std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits(g4particle);
0394   int ng4hits = g4hits.size();
0395   float gpx = g4particle->get_px();
0396   float gpy = g4particle->get_py();
0397   float gpz = g4particle->get_pz();
0398 
0399   PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
0400   float gvx = vtx->get_x();
0401   float gvy = vtx->get_y();
0402   float gvz = vtx->get_z();
0403 
0404   float gfpx = NULL;
0405   float gfpy = NULL;
0406   float gfpz = NULL;
0407   float gfx = NULL;
0408   float gfy = NULL;
0409   float gfz = NULL;
0410 
0411   PHG4Hit* outerhit = trutheval->get_outermost_truth_hit(g4particle);
0412   if (outerhit)
0413     {
0414       gfpx = outerhit->get_px(1);
0415       gfpy = outerhit->get_py(1);
0416       gfpz = outerhit->get_pz(1);
0417       gfx = outerhit->get_x(1);
0418       gfy = outerhit->get_y(1);
0419       gfz = outerhit->get_z(1);
0420     }
0421 
0422   int gembed = trutheval->get_embed(g4particle);
0423 
0424   SvtxTrack* track = trackeval->best_track_from(g4particle);
0425 
0426   float trackID = NAN;
0427   float charge = NAN;
0428   float quality = NAN;
0429   float chisq = NAN;
0430   float ndf = NAN;
0431   float nhits = NAN;
0432   unsigned int layers = 0x0;
0433   float dca = NAN;
0434   float dca2d = NAN;
0435   float dca2dsigma = NAN;
0436   float px = NAN;
0437   float py = NAN;
0438   float pz = NAN;
0439   float pcax = NAN;
0440   float pcay = NAN;
0441   float pcaz = NAN;
0442 
0443   int nfromtruth = -1;
0444 
0445   if (track)
0446     {
0447       trackID = track->get_id();
0448       charge = track->get_charge();
0449       quality = track->get_quality();
0450       chisq = track->get_chisq();
0451       ndf = track->get_ndf();
0452       nhits = track->size_clusters();
0453 
0454       for (SvtxTrack::ConstClusterIter iter = track->begin_clusters();
0455           iter != track->end_clusters(); ++iter)
0456         {
0457           unsigned int cluster_id = *iter;
0458           SvtxCluster* cluster = clustermap->get(cluster_id);
0459           unsigned int layer = cluster->get_layer();
0460           if (layer < 32)
0461             layers |= (0x1 << layer);
0462         }
0463 
0464       dca = track->get_dca();
0465       dca2d = track->get_dca2d();
0466       dca2dsigma = track->get_dca2d_error();
0467       px = track->get_px();
0468       py = track->get_py();
0469       pz = track->get_pz();
0470       pcax = track->get_x();
0471       pcay = track->get_y();
0472       pcaz = track->get_z();
0473 
0474       nfromtruth = trackeval->get_nclusters_contribution(track, g4particle);
0475 
0476       mc_photon.presdphi = track->get_cal_dphi(SvtxTrack::PRES);
0477       mc_photon.presdeta = track->get_cal_deta(SvtxTrack::PRES);
0478       mc_photon.prese3x3 = track->get_cal_energy_3x3(SvtxTrack::PRES);
0479       mc_photon.prese = track->get_cal_cluster_e(SvtxTrack::PRES);
0480 
0481       mc_photon.cemcdphi = track->get_cal_dphi(SvtxTrack::CEMC);
0482       mc_photon.cemcdeta = track->get_cal_deta(SvtxTrack::CEMC);
0483       mc_photon.cemce3x3 = track->get_cal_energy_3x3(SvtxTrack::CEMC);
0484       mc_photon.cemce = track->get_cal_cluster_e(SvtxTrack::CEMC);
0485 
0486       mc_photon.hcalindphi = track->get_cal_dphi(SvtxTrack::HCALIN);
0487       mc_photon.hcalindeta = track->get_cal_deta(SvtxTrack::HCALIN);
0488       mc_photon.hcaline3x3 = track->get_cal_energy_3x3(SvtxTrack::HCALIN);
0489       mc_photon.hcaline = track->get_cal_cluster_e(SvtxTrack::HCALIN);
0490 
0491       mc_photon.hcaloutdphi = track->get_cal_dphi(SvtxTrack::HCALOUT);
0492       mc_photon.hcaloutdeta = track->get_cal_deta(SvtxTrack::HCALOUT);
0493       mc_photon.hcaloute3x3 = track->get_cal_energy_3x3(SvtxTrack::HCALOUT);
0494       mc_photon.hcaloute = track->get_cal_cluster_e(SvtxTrack::HCALOUT);
0495 
0496     }
0497 
0498   eval_photon_proj( //
0499       /*PHG4Particle **/g4particle, //
0500       /*MCPhoton &*/mc_photon,
0501       /*enu_calo*/kCEMC, //
0502       /*const double */gvz,
0503       /*PHCompositeNode **/topNode //
0504       );
0505   eval_photon_proj( //
0506       /*PHG4Particle **/g4particle, //
0507       /*MCPhoton &*/mc_photon,
0508       /*enu_calo*/kHCALIN, //
0509       /*const double */gvz,
0510       /*PHCompositeNode **/topNode //
0511       );
0512 
0513   mc_photon.gtrackID = gtrackID;
0514   mc_photon.gflavor = gflavor;
0515   mc_photon.ng4hits = ng4hits;
0516   mc_photon.gpx = gpx;
0517   mc_photon.gpy = gpy;
0518   mc_photon.gpz = gpz;
0519   mc_photon.gvx = gvx;
0520   mc_photon.gvy = gvy;
0521   mc_photon.gvz = gvz;
0522   mc_photon.gfpx = gfpx;
0523   mc_photon.gfpy = gfpy;
0524   mc_photon.gfpz = gfpz;
0525   mc_photon.gfx = gfx;
0526   mc_photon.gfy = gfy;
0527   mc_photon.gfz = gfz;
0528   mc_photon.gembed = gembed;
0529 //  trk.gprimary = gprimary;
0530   mc_photon.trackID = trackID;
0531   mc_photon.px = px;
0532   mc_photon.py = py;
0533   mc_photon.pz = pz;
0534   mc_photon.charge = charge;
0535   mc_photon.quality = quality;
0536   mc_photon.chisq = chisq;
0537   mc_photon.ndf = ndf;
0538   mc_photon.nhits = nhits;
0539   mc_photon.layers = layers;
0540   mc_photon.dca2d = dca2d;
0541   mc_photon.dca2dsigma = dca2dsigma;
0542   mc_photon.pcax = pcax;
0543   mc_photon.pcay = pcay;
0544   mc_photon.pcaz = pcaz;
0545   mc_photon.nfromtruth = nfromtruth;
0546 
0547 }
0548 
0549 void
0550 EMCalCalib::eval_photon_proj( //
0551     PHG4Particle * g4particle, //
0552     MCPhoton & trk, enu_calo calo_id, //
0553     const double gvz, PHCompositeNode *topNode //
0554     )
0555 // Track projections
0556 {
0557   // hard coded radius
0558   string detector = "InvalidDetector";
0559   double radius = 110;
0560   if (calo_id == kCEMC)
0561     {
0562       detector = "CEMC";
0563       radius = 105;
0564     }
0565   if (calo_id == kHCALIN)
0566     {
0567       detector = "HCALIN";
0568       radius = 125;
0569     }
0570 
0571   if (!_eval_stack)
0572     {
0573       _eval_stack = new SvtxEvalStack(topNode);
0574       _eval_stack->set_strict(false);
0575     }
0576 
0577   // pull the tower geometry
0578   string towergeonodename = "TOWERGEOM_" + detector;
0579   RawTowerGeomContainer *cemc_towergeo = findNode::getClass<
0580       RawTowerGeomContainer>(topNode, towergeonodename.c_str());
0581   assert(cemc_towergeo);
0582 
0583   if (verbosity > 1)
0584     {
0585       cemc_towergeo->identify();
0586     }
0587   // pull the towers
0588   string towernodename = "TOWER_CALIB_" + detector;
0589   RawTowerContainer *cemc_towerList = findNode::getClass<RawTowerContainer>(
0590       topNode, towernodename.c_str());
0591   assert(cemc_towerList);
0592 
0593   // curved tracks inside mag field
0594   // straight projections thereafter
0595   std::vector<double> point;
0596 //  point.assign(3, -9999.);
0597 //  _hough.projectToRadius(track, _magfield, radius, point);
0598 
0599   const double pt = sqrt(
0600       g4particle->get_px() * g4particle->get_px()
0601           + g4particle->get_py() * g4particle->get_py());
0602 
0603   double x = g4particle->get_px() / pt * radius;
0604   double y = g4particle->get_py() / pt * radius;
0605   double z = g4particle->get_pz() / pt * radius + gvz;
0606 
0607   double phi = atan2(y, x);
0608   double eta = asinh(z / sqrt(x * x + y * y));
0609 
0610   // calculate 3x3 tower energy
0611   int binphi = cemc_towergeo->get_phibin(phi);
0612   int bineta = cemc_towergeo->get_etabin(eta);
0613 
0614   double etabin_width = cemc_towergeo->get_etabounds(bineta).second
0615       - cemc_towergeo->get_etabounds(bineta).first;
0616   if (bineta > 1 and bineta < cemc_towergeo->get_etabins() - 1)
0617     etabin_width = (cemc_towergeo->get_etacenter(bineta + 1)
0618         - cemc_towergeo->get_etacenter(bineta - 1)) / 2.;
0619 
0620   double phibin_width = cemc_towergeo->get_phibounds(binphi).second
0621       - cemc_towergeo->get_phibounds(binphi).first;
0622 
0623   assert(etabin_width>0);
0624   assert(phibin_width>0);
0625 
0626   const double etabin_shift = (cemc_towergeo->get_etacenter(bineta) - eta)
0627       / etabin_width;
0628   const double phibin_shift = (fmod(
0629       cemc_towergeo->get_phicenter(binphi) - phi + 5 * M_PI, 2 * M_PI) - M_PI)
0630       / phibin_width;
0631 
0632   if (verbosity > 1)
0633     cout << __PRETTY_FUNCTION__ << " - info - handling track proj. (" << x
0634         << ", " << y << ", " << z << ")" << ", eta " << eta << ", phi" << phi
0635         << ", bineta " << bineta << " - ["
0636         << cemc_towergeo->get_etabounds(bineta).first << ", "
0637         << cemc_towergeo->get_etabounds(bineta).second << "], etabin_shift = "
0638         << etabin_shift << ", binphi " << binphi << " - ["
0639         << cemc_towergeo->get_phibounds(binphi).first << ", "
0640         << cemc_towergeo->get_phibounds(binphi).second << "], phibin_shift = "
0641         << phibin_shift << endl;
0642 
0643   const int bin_search_range = (trk.Max_N_Tower - 1) / 2;
0644   for (int iphi = binphi - bin_search_range; iphi <= binphi + bin_search_range;
0645       ++iphi)
0646 
0647     for (int ieta = bineta - bin_search_range;
0648         ieta <= bineta + bin_search_range; ++ieta)
0649       {
0650 
0651         // wrap around
0652         int wrapphi = iphi;
0653         if (wrapphi < 0)
0654           {
0655             wrapphi = cemc_towergeo->get_phibins() + wrapphi;
0656           }
0657         if (wrapphi >= cemc_towergeo->get_phibins())
0658           {
0659             wrapphi = wrapphi - cemc_towergeo->get_phibins();
0660           }
0661 
0662         // edges
0663         if (ieta < 0)
0664           continue;
0665         if (ieta >= cemc_towergeo->get_etabins())
0666           continue;
0667 
0668         if (verbosity > 1)
0669           cout << __PRETTY_FUNCTION__ << " - info - processing tower"
0670               << ", bineta " << ieta << " - ["
0671               << cemc_towergeo->get_etabounds(ieta).first << ", "
0672               << cemc_towergeo->get_etabounds(ieta).second << "]" << ", binphi "
0673               << wrapphi << " - ["
0674               << cemc_towergeo->get_phibounds(wrapphi).first << ", "
0675               << cemc_towergeo->get_phibounds(wrapphi).second << "]" << endl;
0676 
0677         RawTower* tower = cemc_towerList->getTower(ieta, wrapphi);
0678         double energy = 0;
0679         if (tower)
0680           {
0681             if (verbosity > 1)
0682               cout << __PRETTY_FUNCTION__ << " - info - tower " << ieta << " "
0683                   << wrapphi << " energy = " << tower->get_energy() << endl;
0684 
0685             energy = tower->get_energy();
0686           }
0687 
0688         if (calo_id == kCEMC)
0689           {
0690             trk.cemc_ieta //
0691             [ieta - (bineta - bin_search_range)] //
0692             [iphi - (binphi - bin_search_range)] //
0693             = ieta - bineta + etabin_shift;
0694             trk.cemc_iphi //
0695             [ieta - (bineta - bin_search_range)] //
0696             [iphi - (binphi - bin_search_range)] //
0697             = iphi - binphi + phibin_shift;
0698             trk.cemc_energy //
0699             [ieta - (bineta - bin_search_range)] //
0700             [iphi - (binphi - bin_search_range)] //
0701             = energy;
0702           }
0703         if (calo_id == kHCALIN)
0704           {
0705             trk.hcalin_ieta //
0706             [ieta - (bineta - bin_search_range)] //
0707             [iphi - (binphi - bin_search_range)] //
0708             = ieta - bineta + etabin_shift;
0709             trk.hcalin_iphi //
0710             [ieta - (bineta - bin_search_range)] //
0711             [iphi - (binphi - bin_search_range)] //
0712             = iphi - binphi + phibin_shift;
0713             trk.hcalin_energy //
0714             [ieta - (bineta - bin_search_range)] //
0715             [iphi - (binphi - bin_search_range)] //
0716             = energy;
0717           }
0718 
0719       } //            for (int ieta = bineta-1; ieta < bineta+2; ++ieta) {
0720 
0721 } //       // Track projections
0722 
0723 int
0724 EMCalCalib::Init_MCPhoton(PHCompositeNode *topNode)
0725 {
0726 //  static const int BUFFER_SIZE = 32000;
0727 //  _T_EMCalTrk = new TTree("T_EMCalTrk", "T_EMCalTrk");
0728 //
0729 //  _T_EMCalTrk->Branch("trk.", _trk.ClassName(), &_trk, BUFFER_SIZE);
0730 //  _T_EMCalTrk->Branch("event", &_ievent, "event/I", BUFFER_SIZE);
0731 
0732   return Fun4AllReturnCodes::EVENT_OK;
0733 }
0734 
0735 int
0736 EMCalCalib::process_event_MCPhoton(PHCompositeNode *topNode)
0737 {
0738 
0739   if (verbosity > 2)
0740     cout << "EMCalCalib::process_event_MCPhoton() entered" << endl;
0741 
0742   if (!_eval_stack)
0743     {
0744       _eval_stack = new SvtxEvalStack(topNode);
0745       _eval_stack->set_strict(false);
0746     }
0747 
0748   _mc_photon = findNode::getClass<MCPhoton>(topNode, "MCPhoton");
0749   assert(_mc_photon);
0750 
0751   PHG4TruthInfoContainer* truthinfo =
0752       findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0753   assert(truthinfo);
0754   PHG4TruthInfoContainer::ConstRange range =
0755       truthinfo->GetPrimaryParticleRange();
0756 
0757   assert(range.first != range.second);
0758   // make sure there is at least one primary partcle
0759 
0760   eval_photon(range.first->second, *_mc_photon, topNode);
0761 
0762   return Fun4AllReturnCodes::EVENT_OK;
0763 }
0764 
0765 Fun4AllHistoManager *
0766 EMCalCalib::get_HistoManager()
0767 {
0768 
0769   Fun4AllServer *se = Fun4AllServer::instance();
0770   Fun4AllHistoManager *hm = se->getHistoManager("EMCalAna_HISTOS");
0771 
0772   if (not hm)
0773     {
0774       cout
0775           << "EMCalCalib::get_HistoManager - Making Fun4AllHistoManager EMCalAna_HISTOS"
0776           << endl;
0777       hm = new Fun4AllHistoManager("EMCalAna_HISTOS");
0778       se->registerHistoManager(hm);
0779     }
0780 
0781   assert(hm);
0782 
0783   return hm;
0784 }
0785 
0786 int
0787 EMCalCalib::Init_SF(PHCompositeNode *topNode)
0788 {
0789 
0790   Fun4AllHistoManager *hm = get_HistoManager();
0791   assert(hm);
0792 
0793   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_SF", //
0794       "_h_CEMC_SF", 1000, 0, .2));
0795   hm->registerHisto(new TH1F("EMCalAna_h_HCALIN_SF", //
0796       "_h_HCALIN_SF", 1000, 0, .2));
0797   hm->registerHisto(new TH1F("EMCalAna_h_HCALOUT_SF", //
0798       "_h_HCALOUT_SF", 1000, 0, .2));
0799   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_VSF", //
0800       "_h_CEMC_VSF", 1000, 0, .2));
0801   hm->registerHisto(new TH1F("EMCalAna_h_HCALIN_VSF", //
0802       "_h_HCALIN_VSF", 1000, 0, .2));
0803   hm->registerHisto(new TH1F("EMCalAna_h_HCALOUT_VSF", //
0804       "_h_HCALOUT_VSF", 1000, 0, .2));
0805 
0806   hm->registerHisto(new TH2F("EMCalAna_h_CEMC_RZ", //
0807       "EMCalAna_h_CEMC_RZ;Z (cm);R (cm)", 1200, -300, 300, 600, -000, 300));
0808   hm->registerHisto(new TH2F("EMCalAna_h_HCALIN_RZ", //
0809       "EMCalAna_h_HCALIN_RZ;Z (cm);R (cm)", 1200, -300, 300, 600, -000, 300));
0810   hm->registerHisto(new TH2F("EMCalAna_h_HCALOUT_RZ", //
0811       "EMCalAna_h_HCALOUT_RZ;Z (cm);R (cm)", 1200, -300, 300, 600, -000, 300));
0812 
0813   hm->registerHisto(new TH2F("EMCalAna_h_CEMC_XY", //
0814       "EMCalAna_h_CEMC_XY;X (cm);Y (cm)", 1200, -300, 300, 1200, -300, 300));
0815   hm->registerHisto(new TH2F("EMCalAna_h_HCALIN_XY", //
0816       "EMCalAna_h_HCALIN_XY;X (cm);Y (cm)", 1200, -300, 300, 1200, -300, 300));
0817   hm->registerHisto(new TH2F("EMCalAna_h_HCALOUT_XY", //
0818       "EMCalAna_h_HCALOUT_XY;X (cm);Y (cm)", 1200, -300, 300, 1200, -300, 300));
0819 
0820   return Fun4AllReturnCodes::EVENT_OK;
0821 }
0822 
0823 int
0824 EMCalCalib::process_event_SF(PHCompositeNode *topNode)
0825 {
0826 
0827   if (verbosity > 2)
0828     cout << "EMCalCalib::process_event_SF() entered" << endl;
0829 
0830   TH1F* h = NULL;
0831 
0832   Fun4AllHistoManager *hm = get_HistoManager();
0833   assert(hm);
0834 
0835   double e_hcin = 0.0, e_hcout = 0.0, e_cemc = 0.0;
0836   double ev_hcin = 0.0, ev_hcout = 0.0, ev_cemc = 0.0;
0837   double ea_hcin = 0.0, ea_hcout = 0.0, ea_cemc = 0.0;
0838 
0839   if (_hcalout_hit_container)
0840     {
0841       TH2F * hrz = (TH2F*) hm->getHisto("EMCalAna_h_HCALOUT_RZ");
0842       assert(hrz);
0843       TH2F * hxy = (TH2F*) hm->getHisto("EMCalAna_h_HCALOUT_XY");
0844       assert(hxy);
0845 
0846       PHG4HitContainer::ConstRange hcalout_hit_range =
0847           _hcalout_hit_container->getHits();
0848       for (PHG4HitContainer::ConstIterator hit_iter = hcalout_hit_range.first;
0849           hit_iter != hcalout_hit_range.second; hit_iter++)
0850         {
0851 
0852           PHG4Hit *this_hit = hit_iter->second;
0853           assert(this_hit);
0854 
0855           e_hcout += this_hit->get_edep();
0856           ev_hcout += this_hit->get_light_yield();
0857 
0858           const TVector3 hit(this_hit->get_avg_x(), this_hit->get_avg_y(),
0859               this_hit->get_avg_z());
0860           hrz->Fill(hit.Z(), hit.Perp(), this_hit->get_light_yield());
0861           hxy->Fill(hit.X(), hit.Y(), this_hit->get_light_yield());
0862         }
0863     }
0864 
0865   if (_hcalin_hit_container)
0866     {
0867       TH2F * hrz = (TH2F*) hm->getHisto("EMCalAna_h_HCALIN_RZ");
0868       assert(hrz);
0869       TH2F * hxy = (TH2F*) hm->getHisto("EMCalAna_h_HCALIN_XY");
0870       assert(hxy);
0871 
0872       PHG4HitContainer::ConstRange hcalin_hit_range =
0873           _hcalin_hit_container->getHits();
0874       for (PHG4HitContainer::ConstIterator hit_iter = hcalin_hit_range.first;
0875           hit_iter != hcalin_hit_range.second; hit_iter++)
0876         {
0877 
0878           PHG4Hit *this_hit = hit_iter->second;
0879           assert(this_hit);
0880 
0881           e_hcin += this_hit->get_edep();
0882           ev_hcin += this_hit->get_light_yield();
0883 
0884           const TVector3 hit(this_hit->get_avg_x(), this_hit->get_avg_y(),
0885               this_hit->get_avg_z());
0886           hrz->Fill(hit.Z(), hit.Perp(), this_hit->get_light_yield());
0887           hxy->Fill(hit.X(), hit.Y(), this_hit->get_light_yield());
0888 
0889         }
0890     }
0891 
0892   if (_cemc_hit_container)
0893     {
0894       TH2F * hrz = (TH2F*) hm->getHisto("EMCalAna_h_CEMC_RZ");
0895       assert(hrz);
0896       TH2F * hxy = (TH2F*) hm->getHisto("EMCalAna_h_CEMC_XY");
0897       assert(hxy);
0898 
0899       PHG4HitContainer::ConstRange cemc_hit_range =
0900           _cemc_hit_container->getHits();
0901       for (PHG4HitContainer::ConstIterator hit_iter = cemc_hit_range.first;
0902           hit_iter != cemc_hit_range.second; hit_iter++)
0903         {
0904 
0905           PHG4Hit *this_hit = hit_iter->second;
0906           assert(this_hit);
0907 
0908           e_cemc += this_hit->get_edep();
0909           ev_cemc += this_hit->get_light_yield();
0910 
0911           const TVector3 hit(this_hit->get_avg_x(), this_hit->get_avg_y(),
0912               this_hit->get_avg_z());
0913           hrz->Fill(hit.Z(), hit.Perp(), this_hit->get_light_yield());
0914           hxy->Fill(hit.X(), hit.Y(), this_hit->get_light_yield());
0915         }
0916     }
0917 
0918   if (_hcalout_abs_hit_container)
0919     {
0920       PHG4HitContainer::ConstRange hcalout_abs_hit_range =
0921           _hcalout_abs_hit_container->getHits();
0922       for (PHG4HitContainer::ConstIterator hit_iter =
0923           hcalout_abs_hit_range.first; hit_iter != hcalout_abs_hit_range.second;
0924           hit_iter++)
0925         {
0926 
0927           PHG4Hit *this_hit = hit_iter->second;
0928           assert(this_hit);
0929 
0930           ea_hcout += this_hit->get_edep();
0931 
0932         }
0933     }
0934 
0935   if (_hcalin_abs_hit_container)
0936     {
0937       PHG4HitContainer::ConstRange hcalin_abs_hit_range =
0938           _hcalin_abs_hit_container->getHits();
0939       for (PHG4HitContainer::ConstIterator hit_iter = hcalin_abs_hit_range.first;
0940           hit_iter != hcalin_abs_hit_range.second; hit_iter++)
0941         {
0942 
0943           PHG4Hit *this_hit = hit_iter->second;
0944           assert(this_hit);
0945 
0946           ea_hcin += this_hit->get_edep();
0947         }
0948     }
0949 
0950   if (_cemc_abs_hit_container)
0951     {
0952       PHG4HitContainer::ConstRange cemc_abs_hit_range =
0953           _cemc_abs_hit_container->getHits();
0954       for (PHG4HitContainer::ConstIterator hit_iter = cemc_abs_hit_range.first;
0955           hit_iter != cemc_abs_hit_range.second; hit_iter++)
0956         {
0957 
0958           PHG4Hit *this_hit = hit_iter->second;
0959           assert(this_hit);
0960 
0961           ea_cemc += this_hit->get_edep();
0962 
0963         }
0964     }
0965 
0966   h = (TH1F*) hm->getHisto("EMCalAna_h_CEMC_SF");
0967   assert(h);
0968   h->Fill(e_cemc / (e_cemc + ea_cemc) + 1e-9);
0969   h = (TH1F*) hm->getHisto("EMCalAna_h_CEMC_VSF");
0970   assert(h);
0971   h->Fill(ev_cemc / (e_cemc + ea_cemc) + 1e-9);
0972 
0973   h = (TH1F*) hm->getHisto("EMCalAna_h_HCALOUT_SF");
0974   assert(h);
0975   h->Fill(e_hcout / (e_hcout + ea_hcout) + 1e-9);
0976   h = (TH1F*) hm->getHisto("EMCalAna_h_HCALOUT_VSF");
0977   assert(h);
0978   h->Fill(ev_hcout / (e_hcout + ea_hcout) + 1e-9);
0979 
0980   h = (TH1F*) hm->getHisto("EMCalAna_h_HCALIN_SF");
0981   assert(h);
0982   h->Fill(e_hcin / (e_hcin + ea_hcin) + 1e-9);
0983   h = (TH1F*) hm->getHisto("EMCalAna_h_HCALIN_VSF");
0984   assert(h);
0985   h->Fill(ev_hcin / (e_hcin + ea_hcin) + 1e-9);
0986 
0987   return Fun4AllReturnCodes::EVENT_OK;
0988 }
0989 
0990 int
0991 EMCalCalib::Init_Tower(PHCompositeNode *topNode)
0992 {
0993 
0994   Fun4AllHistoManager *hm = get_HistoManager();
0995   assert(hm);
0996 
0997   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_1x1", //
0998       "h_CEMC_TOWER_1x1", 5000, 0, 50));
0999   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_1x1_max", //
1000       "EMCalAna_h_CEMC_TOWER_1x1_max", 5000, 0, 50));
1001   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_1x1_max_trigger_ADC", //
1002       "EMCalAna_h_CEMC_TOWER_1x1_max_trigger_ADC", 5000, 0, 50));
1003 
1004   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_2x2", //
1005       "h_CEMC_TOWER_2x2", 5000, 0, 50));
1006   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_2x2_max", //
1007       "EMCalAna_h_CEMC_TOWER_2x2_max", 5000, 0, 50));
1008   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_2x2_max_trigger_ADC", //
1009       "EMCalAna_h_CEMC_TOWER_2x2_max_trigger_ADC", 5000, 0, 50));
1010   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_2x2_slide2_max_trigger_ADC", //
1011       "EMCalAna_h_CEMC_TOWER_2x2_slide2_max_trigger_ADC", 5000, 0, 50));
1012 
1013   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_3x3", //
1014       "h_CEMC_TOWER_3x3", 5000, 0, 50));
1015   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_3x3_max", //
1016       "EMCalAna_h_CEMC_TOWER_3x3_max", 5000, 0, 50));
1017   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_3x3_max_trigger_ADC", //
1018       "EMCalAna_h_CEMC_TOWER_3x3_max_trigger_ADC", 5000, 0, 50));
1019 
1020   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_4x4", //
1021       "h_CEMC_TOWER_4x4", 5000, 0, 50));
1022   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_4x4_max", //
1023       "EMCalAna_h_CEMC_TOWER_4x4_max", 5000, 0, 50));
1024   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_4x4_max_trigger_ADC", //
1025       "EMCalAna_h_CEMC_TOWER_4x4_max_trigger_ADC", 5000, 0, 50));
1026   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_4x4_slide2_max_trigger_ADC", //
1027       "EMCalAna_h_CEMC_TOWER_4x4_slide2_max_trigger_ADC", 5000, 0, 50));
1028 
1029   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_5x5", //
1030       "h_CEMC_TOWER_4x4", 5000, 0, 50));
1031   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_5x5_max", //
1032       "EMCalAna_h_CEMC_TOWER_5x5_max", 5000, 0, 50));
1033   hm->registerHisto(new TH1F("EMCalAna_h_CEMC_TOWER_5x5_max_trigger_ADC", //
1034       "EMCalAna_h_CEMC_TOWER_5x5_max_trigger_ADC", 5000, 0, 50));
1035 
1036   return Fun4AllReturnCodes::EVENT_OK;
1037 }
1038 
1039 int
1040 EMCalCalib::process_event_Tower(PHCompositeNode *topNode)
1041 {
1042   const string detector("CEMC");
1043 
1044   if (verbosity > 2)
1045     cout << "EMCalCalib::process_event_SF() entered" << endl;
1046 
1047   string towernodename = "TOWER_CALIB_" + detector;
1048   // Grab the towers
1049   RawTowerContainer* towers = findNode::getClass<RawTowerContainer>(topNode,
1050       towernodename.c_str());
1051   if (!towers)
1052     {
1053       std::cout << PHWHERE << ": Could not find node " << towernodename.c_str()
1054           << std::endl;
1055       return Fun4AllReturnCodes::DISCARDEVENT;
1056     }
1057   string towergeomnodename = "TOWERGEOM_" + detector;
1058   RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(
1059       topNode, towergeomnodename.c_str());
1060   if (!towergeom)
1061     {
1062       cout << PHWHERE << ": Could not find node " << towergeomnodename.c_str()
1063           << endl;
1064       return Fun4AllReturnCodes::ABORTEVENT;
1065     }
1066 
1067   static const double trigger_ADC_bin = 45. / 256.; //8-bit ADC max to 45 GeV
1068   static const int max_size = 5;
1069   map<int, string> size_label;
1070   size_label[1] = "1x1";
1071   size_label[2] = "2x2";
1072   size_label[3] = "3x3";
1073   size_label[4] = "4x4";
1074   size_label[5] = "5x5";
1075   map<int, double> max_energy;
1076   map<int, double> max_energy_trigger_ADC;
1077   map<int, double> slide2_max_energy_trigger_ADC;
1078   map<int, TH1F*> energy_hist_list;
1079   map<int, TH1F*> max_energy_hist_list;
1080   map<int, TH1F*> max_energy_trigger_ADC_hist_list;
1081   map<int, TH1F*> slide2_max_energy_trigger_ADC_hist_list;
1082 
1083   Fun4AllHistoManager *hm = get_HistoManager();
1084   assert(hm);
1085   for (int size = 1; size <= max_size; ++size)
1086     {
1087       max_energy[size] = 0;
1088       max_energy_trigger_ADC[size] = 0;
1089 
1090       TH1F* h = NULL;
1091 
1092       h = (TH1F*) hm->getHisto("EMCalAna_h_CEMC_TOWER_" + size_label[size]);
1093       assert(h);
1094       energy_hist_list[size] = h;
1095       h = (TH1F*) hm->getHisto(
1096           "EMCalAna_h_CEMC_TOWER_" + size_label[size] + "_max");
1097       assert(h);
1098       max_energy_hist_list[size] = h;
1099 
1100       h = (TH1F*) hm->getHisto(
1101           "EMCalAna_h_CEMC_TOWER_" + size_label[size] + "_max_trigger_ADC");
1102       assert(h);
1103       max_energy_trigger_ADC_hist_list[size] = h;
1104 
1105       if (size == 2 or size == 4)
1106         {
1107           // sliding window made from 2x2 sums
1108           slide2_max_energy_trigger_ADC[size] = 0;
1109 
1110           h = (TH1F*) hm->getHisto(
1111               "EMCalAna_h_CEMC_TOWER_" + size_label[size]
1112                   + "_slide2_max_trigger_ADC");
1113           assert(h);
1114           slide2_max_energy_trigger_ADC_hist_list[size] = h;
1115 
1116         }
1117 
1118     }
1119 
1120   for (int binphi = 0; binphi < towergeom->get_phibins(); ++binphi)
1121     {
1122       for (int bineta = 0; bineta < towergeom->get_etabins(); ++bineta)
1123         {
1124           for (int size = 1; size <= max_size; ++size)
1125             {
1126               double energy = 0;
1127               double energy_trigger_ADC = 0;
1128               double slide2_energy_trigger_ADC = 0;
1129 
1130               for (int iphi = binphi; iphi < binphi + size; ++iphi)
1131                 {
1132                   for (int ieta = bineta; ieta < bineta + size; ++ieta)
1133                     {
1134                       if (ieta > towergeom->get_etabins())
1135                         continue;
1136 
1137                       // wrap around
1138                       int wrapphi = iphi;
1139                       assert(wrapphi >= 0);
1140                       if (wrapphi >= towergeom->get_phibins())
1141                         {
1142                           wrapphi = wrapphi - towergeom->get_phibins();
1143                         }
1144 
1145                       RawTower* tower = towers->getTower(ieta, wrapphi);
1146 
1147                       if (tower)
1148                         {
1149                           const double e_intput = tower->get_energy();
1150 
1151                           const double e_trigger_ADC = round(
1152                               e_intput / trigger_ADC_bin) * trigger_ADC_bin;
1153 
1154                           energy += e_intput;
1155                           energy_trigger_ADC += e_trigger_ADC;
1156 
1157                           if ((size == 2 or size == 4) and (binphi % 2 == 0)
1158                               and (bineta % 2 == 0))
1159                             {
1160                               // sliding window made from 2x2 sums
1161 
1162                               slide2_energy_trigger_ADC += e_trigger_ADC;
1163                             }
1164                         }
1165                     }
1166                 }
1167 
1168               energy_hist_list[size]->Fill(energy);
1169 
1170               if (energy > max_energy[size])
1171                 max_energy[size] = energy;
1172               if (energy_trigger_ADC > max_energy_trigger_ADC[size])
1173                 max_energy_trigger_ADC[size] = energy_trigger_ADC;
1174 
1175               if ((size == 2 or size == 4) and (binphi % 2 == 0)
1176                   and (bineta % 2 == 0))
1177                 {
1178                   // sliding window made from 2x2 sums
1179 
1180                   if (slide2_energy_trigger_ADC
1181                       > slide2_max_energy_trigger_ADC[size])
1182                     slide2_max_energy_trigger_ADC[size] =
1183                         slide2_energy_trigger_ADC;
1184                 }
1185 
1186             } //          for (int size = 1; size <= 4; ++size)
1187         }
1188     }
1189 
1190   for (int size = 1; size <= max_size; ++size)
1191     {
1192       max_energy_hist_list[size]->Fill(max_energy[size]);
1193       max_energy_trigger_ADC_hist_list[size]->Fill(
1194           max_energy_trigger_ADC[size]);
1195 
1196       if (size == 2 or size == 4)
1197         {
1198           // sliding window made from 2x2 sums
1199           slide2_max_energy_trigger_ADC_hist_list[size]->Fill(
1200               slide2_max_energy_trigger_ADC[size]);
1201         }
1202     }
1203   return Fun4AllReturnCodes::EVENT_OK;
1204 }
[ 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!