Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-03 08:12:42 

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