File indexing completed on 2025-08-06 08:13:14
0001 #include "CaloEvaluator.h"
0002
0003 #include "CaloEvalStack.h"
0004 #include "CaloRawClusterEval.h"
0005 #include "CaloRawTowerEval.h"
0006 #include "CaloTruthEval.h"
0007
0008 #include <g4main/PHG4Particle.h>
0009 #include <g4main/PHG4Shower.h>
0010 #include <g4main/PHG4TruthInfoContainer.h>
0011 #include <g4main/PHG4VtxPoint.h>
0012
0013 #include <g4vertex/GlobalVertex.h>
0014 #include <g4vertex/GlobalVertexMap.h>
0015
0016 #include <calobase/RawCluster.h>
0017 #include <calobase/RawClusterContainer.h>
0018 #include <calobase/RawClusterUtility.h>
0019 #include <calobase/RawTower.h>
0020 #include <calobase/RawTowerContainer.h>
0021 #include <calobase/RawTowerGeom.h>
0022 #include <calobase/RawTowerGeomContainer.h>
0023
0024 #include <fun4all/Fun4AllReturnCodes.h>
0025 #include <fun4all/SubsysReco.h>
0026
0027 #include <phool/getClass.h>
0028 #include <phool/phool.h>
0029
0030 #include <TFile.h>
0031 #include <TNtuple.h>
0032 #include <TTree.h>
0033
0034 #include <CLHEP/Vector/ThreeVector.h>
0035
0036 #include <cmath>
0037 #include <cstdlib>
0038 #include <iostream>
0039 #include <set>
0040 #include <utility>
0041
0042 using namespace std;
0043
0044 CaloEvaluator::CaloEvaluator(const string& name, const string& caloname, const string& filename)
0045 : SubsysReco(name)
0046 , _caloname(caloname)
0047 , _ievent(0)
0048 , _towerID_debug(0)
0049 , _ieta_debug(0)
0050 , _iphi_debug(0)
0051 , _eta_debug(0)
0052 , _phi_debug(0)
0053 , _e_debug(0)
0054 , _x_debug(0)
0055 , _y_debug(0)
0056 , _z_debug(0)
0057 , _truth_trace_embed_flags()
0058 , _truth_e_threshold(0.0)
0059 ,
0060 _reco_e_threshold(0.0)
0061 ,
0062 _caloevalstack(nullptr)
0063 , _strict(false)
0064 , _do_gpoint_eval(false)
0065 , _do_gshower_eval(true)
0066 , _do_tower_eval(false)
0067 , _do_cluster_eval(true)
0068 , _ntp_gpoint(nullptr)
0069 , _ntp_gshower(nullptr)
0070 , _ntp_tower(nullptr)
0071 , _tower_debug(nullptr)
0072 , _ntp_cluster(nullptr)
0073 , _filename(filename)
0074 , _tfile(nullptr)
0075 {
0076 }
0077
0078 int CaloEvaluator::Init(PHCompositeNode* )
0079 {
0080 _ievent = 0;
0081
0082 _tfile = new TFile(_filename.c_str(), "RECREATE");
0083
0084 if (_do_gpoint_eval) _ntp_gpoint = new TNtuple("ntp_gpoint", "primary vertex => best (first) vertex",
0085 "event:gvx:gvy:gvz:"
0086 "vx:vy:vz");
0087
0088 if (_do_gshower_eval) _ntp_gshower = new TNtuple("ntp_gshower", "truth shower => best cluster",
0089 "event:gparticleID:gflavor:gnhits:"
0090 "geta:gphi:ge:gpt:gvx:gvy:gvz:gembed:gedep:"
0091 "clusterID:ntowers:eta:x:y:z:phi:e:efromtruth");
0092
0093
0094 if (_do_tower_eval)
0095 {
0096 _ntp_tower = new TNtuple("ntp_tower", "tower => max truth primary",
0097 "event:towerID:ieta:iphi:eta:phi:e:x:y:z:"
0098 "gparticleID:gflavor:gnhits:"
0099 "geta:gphi:ge:gpt:gvx:gvy:gvz:"
0100 "gembed:gedep:"
0101 "efromtruth");
0102
0103
0104 _tower_debug = new TTree("tower_debug", "tower => max truth primary");
0105
0106 _tower_debug->Branch("event", &_ievent, "event/I");
0107 _tower_debug->Branch("towerID", &_towerID_debug, "towerID/I");
0108 _tower_debug->Branch("ieta", &_ieta_debug, "ieta/I");
0109 _tower_debug->Branch("iphi", &_iphi_debug, "iphi/I");
0110 _tower_debug->Branch("eta", &_eta_debug, "eta/F");
0111 _tower_debug->Branch("phi", &_phi_debug, "phi/F");
0112 _tower_debug->Branch("e", &_e_debug, "e/F");
0113 _tower_debug->Branch("x", &_x_debug, "x/F");
0114 _tower_debug->Branch("y", &_y_debug, "y/F");
0115 _tower_debug->Branch("z", &_z_debug, "z/F");
0116 }
0117
0118 if (_do_cluster_eval) {
0119 _ntp_cluster = new TNtuple("ntp_cluster", "cluster => max truth primary",
0120 "event:clusterID:ntowers:eta_detector:eta:x:y:z:phi:e:ecore:"
0121 "gparticleID:gflavor:gnhits:"
0122 "geta:gphi:ge:gpt:gvx:gvy:gvz:"
0123 "gembed:gedep:"
0124 "efromtruth");
0125
0126 _cluster_tower_info = new TTree("cluster_tower_info", "Cluster Tower Info");
0127 _cluster_tower_info->Branch("towerEtas", &_towerEtas);
0128 _cluster_tower_info->Branch("towerPhis", &_towerPhis);
0129 _cluster_tower_info->Branch("towerEnergies", &_towerEnergies);
0130 }
0131
0132
0133 set_reco_tracing_energy_threshold(1);
0134
0135 return Fun4AllReturnCodes::EVENT_OK;
0136 }
0137
0138 int CaloEvaluator::process_event(PHCompositeNode* topNode)
0139 {
0140 if (!_caloevalstack)
0141 {
0142 _caloevalstack = new CaloEvalStack(topNode, _caloname);
0143 _caloevalstack->set_strict(_strict);
0144 _caloevalstack->set_verbosity(Verbosity() + 1);
0145 }
0146 else
0147 {
0148 _caloevalstack->next_event(topNode);
0149 }
0150
0151
0152
0153
0154
0155 printInputInfo(topNode);
0156
0157
0158
0159
0160
0161 fillOutputNtuples(topNode);
0162
0163
0164
0165
0166
0167 printOutputInfo(topNode);
0168
0169 ++_ievent;
0170
0171 return Fun4AllReturnCodes::EVENT_OK;
0172 }
0173
0174 int CaloEvaluator::End(PHCompositeNode* )
0175 {
0176 _tfile->cd();
0177
0178 if (_do_gpoint_eval) _ntp_gpoint->Write();
0179 if (_do_gshower_eval) _ntp_gshower->Write();
0180 if (_do_tower_eval)
0181 {
0182 _ntp_tower->Write();
0183 _tower_debug->Write();
0184 }
0185 if (_do_cluster_eval) {
0186 _ntp_cluster->Write();
0187 _cluster_tower_info->Write();
0188 }
0189
0190 _tfile->Close();
0191
0192 delete _tfile;
0193
0194 if (Verbosity() > 0)
0195 {
0196 cout << "========================= " << Name() << "::End() ============================" << endl;
0197 cout << " " << _ievent << " events of output written to: " << _filename << endl;
0198 cout << "===========================================================================" << endl;
0199 }
0200
0201 if (_caloevalstack) delete _caloevalstack;
0202
0203 return Fun4AllReturnCodes::EVENT_OK;
0204 }
0205
0206 void CaloEvaluator::printInputInfo(PHCompositeNode* topNode)
0207 {
0208 if (Verbosity() > 2) cout << "CaloEvaluator::printInputInfo() entered" << endl;
0209
0210
0211
0212 if (Verbosity() > 1)
0213 {
0214 cout << endl;
0215 cout << PHWHERE << " NEW INPUT FOR EVENT " << _ievent << endl;
0216 cout << endl;
0217
0218
0219 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0220 if (!truthinfo)
0221 {
0222 cerr << PHWHERE << " ERROR: Can't find G4TruthInfo" << endl;
0223 exit(-1);
0224 }
0225
0226 cout << Name() << ": PHG4TruthInfoContainer contents: " << endl;
0227
0228 PHG4TruthInfoContainer::Range truthrange = truthinfo->GetParticleRange();
0229 for (PHG4TruthInfoContainer::Iterator truthiter = truthrange.first;
0230 truthiter != truthrange.second;
0231 ++truthiter)
0232 {
0233 PHG4Particle* particle = truthiter->second;
0234
0235 cout << truthiter->first << " => pid: " << particle->get_pid()
0236 << " pt: " << sqrt(pow(particle->get_px(), 2) + pow(particle->get_py(), 2)) << endl;
0237 }
0238 }
0239
0240 return;
0241 }
0242
0243 void CaloEvaluator::printOutputInfo(PHCompositeNode* topNode)
0244 {
0245 if (Verbosity() > 2) cout << "CaloEvaluator::printOutputInfo() entered" << endl;
0246
0247 CaloRawClusterEval* clustereval = _caloevalstack->get_rawcluster_eval();
0248 CaloTruthEval* trutheval = _caloevalstack->get_truth_eval();
0249
0250
0251
0252
0253
0254 if (Verbosity() > 1)
0255 {
0256
0257 cout << endl;
0258 cout << PHWHERE << " NEW OUTPUT FOR EVENT " << _ievent << endl;
0259 cout << endl;
0260
0261
0262 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0263 if (!truthinfo)
0264 {
0265 cerr << PHWHERE << " ERROR: Can't find G4TruthInfo" << endl;
0266 exit(-1);
0267 }
0268
0269
0270 GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0271
0272 PHG4VtxPoint* gvertex = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
0273 float gvx = gvertex->get_x();
0274 float gvy = gvertex->get_y();
0275 float gvz = gvertex->get_z();
0276
0277 float vx = NAN;
0278 float vy = NAN;
0279 float vz = NAN;
0280 if (vertexmap)
0281 {
0282 if (!vertexmap->empty())
0283 {
0284 GlobalVertex* vertex = (vertexmap->begin()->second);
0285
0286 vx = vertex->get_x();
0287 vy = vertex->get_y();
0288 vz = vertex->get_z();
0289 }
0290 }
0291
0292 cout << "vtrue = (" << gvx << "," << gvy << "," << gvz << ") => vreco = (" << vx << "," << vy << "," << vz << ")" << endl;
0293
0294 PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
0295 for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0296 iter != range.second;
0297 ++iter)
0298 {
0299 PHG4Particle* primary = iter->second;
0300
0301 cout << endl;
0302
0303 cout << "===Primary PHG4Particle=========================================" << endl;
0304 cout << " particle id = " << primary->get_track_id() << endl;
0305 cout << " flavor = " << primary->get_pid() << endl;
0306 cout << " (px,py,pz,e) = (";
0307
0308 float gpx = primary->get_px();
0309 float gpy = primary->get_py();
0310 float gpz = primary->get_pz();
0311 float ge = primary->get_e();
0312
0313 cout.width(5);
0314 cout << gpx;
0315 cout << ",";
0316 cout.width(5);
0317 cout << gpy;
0318 cout << ",";
0319 cout.width(5);
0320 cout << gpz;
0321 cout << ",";
0322 cout.width(5);
0323 cout << ge;
0324 cout << ")" << endl;
0325
0326 float gpt = sqrt(gpx * gpx + gpy * gpy);
0327 float geta = NAN;
0328 if (gpt != 0.0) geta = asinh(gpz / gpt);
0329 float gphi = atan2(gpy, gpx);
0330
0331 cout << "(eta,phi,e,pt) = (";
0332 cout.width(5);
0333 cout << geta;
0334 cout << ",";
0335 cout.width(5);
0336 cout << gphi;
0337 cout << ",";
0338 cout.width(5);
0339 cout << ge;
0340 cout << ",";
0341 cout.width(5);
0342 cout << gpt;
0343 cout << ")" << endl;
0344
0345 PHG4VtxPoint* vtx = trutheval->get_vertex(primary);
0346 gvx = vtx->get_x();
0347 gvy = vtx->get_y();
0348 gvz = vtx->get_z();
0349
0350 cout << " vtrue = (";
0351 cout.width(5);
0352 cout << gvx;
0353 cout << ",";
0354 cout.width(5);
0355 cout << gvy;
0356 cout << ",";
0357 cout.width(5);
0358 cout << gvz;
0359 cout << ")" << endl;
0360
0361 cout << " embed = " << trutheval->get_embed(primary) << endl;
0362 cout << " edep = " << trutheval->get_shower_energy_deposit(primary) << endl;
0363
0364 std::set<RawCluster*> clusters = clustereval->all_clusters_from(primary);
0365 for (std::set<RawCluster*>::iterator clusiter = clusters.begin();
0366 clusiter != clusters.end();
0367 ++clusiter)
0368 {
0369 RawCluster* cluster = (*clusiter);
0370
0371 float ntowers = cluster->getNTowers();
0372 float x = cluster->get_x();
0373 float y = cluster->get_y();
0374 float z = cluster->get_z();
0375 float phi = cluster->get_phi();
0376 float e = cluster->get_energy();
0377
0378 float efromtruth = clustereval->get_energy_contribution(cluster, primary);
0379
0380 cout << " => #" << cluster->get_id() << " (x,y,z,phi,e) = (";
0381 cout.width(5);
0382 cout << x;
0383 cout << ",";
0384 cout.width(5);
0385 cout << y;
0386 cout << ",";
0387 cout.width(5);
0388 cout << z;
0389 cout << ",";
0390 cout.width(5);
0391 cout << phi;
0392 cout << ",";
0393 cout.width(5);
0394 cout << e;
0395 cout << "), ntowers = " << ntowers << ", efromtruth = " << efromtruth << endl;
0396 }
0397 }
0398 cout << endl;
0399 }
0400
0401 return;
0402 }
0403
0404 void CaloEvaluator::fillOutputNtuples(PHCompositeNode* topNode)
0405 {
0406 if (Verbosity() > 2) cout << "CaloEvaluator::fillOutputNtuples() entered" << endl;
0407
0408 CaloRawClusterEval* clustereval = _caloevalstack->get_rawcluster_eval();
0409 CaloRawTowerEval* towereval = _caloevalstack->get_rawtower_eval();
0410 CaloTruthEval* trutheval = _caloevalstack->get_truth_eval();
0411
0412
0413
0414
0415
0416 if (_do_gpoint_eval)
0417 {
0418
0419 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0420 if (!truthinfo)
0421 {
0422 cerr << PHWHERE << " ERROR: Can't find G4TruthInfo" << endl;
0423 exit(-1);
0424 }
0425
0426
0427 GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0428
0429 PHG4VtxPoint* gvertex = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
0430 float gvx = gvertex->get_x();
0431 float gvy = gvertex->get_y();
0432 float gvz = gvertex->get_z();
0433
0434 float vx = NAN;
0435 float vy = NAN;
0436 float vz = NAN;
0437 if (vertexmap)
0438 {
0439 if (!vertexmap->empty())
0440 {
0441 GlobalVertex* vertex = (vertexmap->begin()->second);
0442
0443 vx = vertex->get_x();
0444 vy = vertex->get_y();
0445 vz = vertex->get_z();
0446 }
0447 }
0448
0449 float gpoint_data[7] = {(float) _ievent,
0450 gvx,
0451 gvy,
0452 gvz,
0453 vx,
0454 vy,
0455 vz};
0456
0457 _ntp_gpoint->Fill(gpoint_data);
0458 }
0459
0460
0461
0462
0463
0464 if (_ntp_gshower)
0465 {
0466 if (Verbosity() > 1) cout << Name() << " CaloEvaluator::filling gshower ntuple..." << endl;
0467
0468 GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0469
0470 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0471 if (!truthinfo)
0472 {
0473 cerr << PHWHERE << " ERROR: Can't find G4TruthInfo" << endl;
0474 exit(-1);
0475 }
0476
0477 PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
0478 for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0479 iter != range.second;
0480 ++iter)
0481 {
0482 PHG4Particle* primary = iter->second;
0483
0484 if (primary->get_e() < _truth_e_threshold) continue;
0485
0486 if (!_truth_trace_embed_flags.empty())
0487 {
0488 if (_truth_trace_embed_flags.find(trutheval->get_embed(primary)) ==
0489 _truth_trace_embed_flags.end()) continue;
0490 }
0491
0492 float gparticleID = primary->get_track_id();
0493 float gflavor = primary->get_pid();
0494
0495 PHG4Shower* shower = trutheval->get_primary_shower(primary);
0496 float gnhits = NAN;
0497 if (shower)
0498 gnhits = shower->get_nhits(trutheval->get_caloid());
0499 else
0500 gnhits = 0.0;
0501 float gpx = primary->get_px();
0502 float gpy = primary->get_py();
0503 float gpz = primary->get_pz();
0504 float ge = primary->get_e();
0505
0506 float gpt = sqrt(gpx * gpx + gpy * gpy);
0507 float geta = NAN;
0508 if (gpt != 0.0) geta = asinh(gpz / gpt);
0509 float gphi = atan2(gpy, gpx);
0510
0511 PHG4VtxPoint* vtx = trutheval->get_vertex(primary);
0512 float gvx = vtx->get_x();
0513 float gvy = vtx->get_y();
0514 float gvz = vtx->get_z();
0515
0516 float gembed = trutheval->get_embed(primary);
0517 float gedep = trutheval->get_shower_energy_deposit(primary);
0518
0519 RawCluster* cluster = clustereval->best_cluster_from(primary);
0520
0521 float clusterID = NAN;
0522 float ntowers = NAN;
0523 float eta = NAN;
0524 float x = NAN;
0525 float y = NAN;
0526 float z = NAN;
0527 float phi = NAN;
0528 float e = NAN;
0529
0530 float efromtruth = NAN;
0531
0532 if (cluster)
0533 {
0534 clusterID = cluster->get_id();
0535 ntowers = cluster->getNTowers();
0536 x = cluster->get_x();
0537 y = cluster->get_y();
0538 z = cluster->get_z();
0539 phi = cluster->get_phi();
0540 e = cluster->get_energy();
0541
0542 efromtruth = clustereval->get_energy_contribution(cluster, primary);
0543
0544
0545 if (vertexmap)
0546 {
0547 if (!vertexmap->empty())
0548 {
0549 GlobalVertex* vertex = (vertexmap->begin()->second);
0550
0551 eta =
0552 RawClusterUtility::GetPseudorapidity(
0553 *cluster,
0554 CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
0555 }
0556 }
0557 }
0558
0559 float shower_data[] = {(float) _ievent,
0560 gparticleID,
0561 gflavor,
0562 gnhits,
0563 geta,
0564 gphi,
0565 ge,
0566 gpt,
0567 gvx,
0568 gvy,
0569 gvz,
0570 gembed,
0571 gedep,
0572 clusterID,
0573 ntowers,
0574 eta,
0575 x,
0576 y,
0577 z,
0578 phi,
0579 e,
0580 efromtruth};
0581
0582 _ntp_gshower->Fill(shower_data);
0583 }
0584 }
0585
0586
0587
0588
0589
0590 if (_do_tower_eval)
0591 {
0592 if (Verbosity() > 1) cout << "CaloEvaluator::filling tower ntuple..." << endl;
0593
0594 string towernode = "TOWER_CALIB_" + _caloname;
0595 RawTowerContainer* towers = findNode::getClass<RawTowerContainer>(topNode, towernode.c_str());
0596 if (!towers)
0597 {
0598 cerr << PHWHERE << " ERROR: Can't find " << towernode << endl;
0599 exit(-1);
0600 }
0601
0602 string towergeomnode = "TOWERGEOM_" + _caloname;
0603 RawTowerGeomContainer* towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnode.c_str());
0604 if (!towergeom)
0605 {
0606 cerr << PHWHERE << " ERROR: Can't find " << towergeomnode << endl;
0607 exit(-1);
0608 }
0609
0610 RawTowerContainer::ConstRange begin_end = towers->getTowers();
0611 RawTowerContainer::ConstIterator rtiter;
0612 for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0613 {
0614 RawTower* tower = rtiter->second;
0615
0616 if (tower->get_energy() < _reco_e_threshold) continue;
0617
0618 RawTowerGeom* tower_geom = towergeom->get_tower_geometry(tower->get_id());
0619 if (!tower_geom)
0620 {
0621 cerr << PHWHERE << " ERROR: Can't find tower geometry for this tower hit: ";
0622 tower->identify();
0623 exit(-1);
0624 }
0625
0626
0627 const float towerid = tower->get_id();
0628 const float ieta = tower->get_bineta();
0629 const float iphi = tower->get_binphi();
0630 const float eta = tower_geom->get_eta();
0631 const float phi = tower_geom->get_phi();
0632 const float e = tower->get_energy();
0633 const float x = tower_geom->get_center_x();
0634 const float y = tower_geom->get_center_y();
0635 const float z = tower_geom->get_center_z();
0636
0637
0638 _towerID_debug = tower->get_id();
0639 _ieta_debug = tower->get_bineta();
0640 _iphi_debug = tower->get_binphi();
0641 _eta_debug = tower_geom->get_eta();
0642 _phi_debug = tower_geom->get_phi();
0643 _e_debug = tower->get_energy();
0644 _x_debug = tower_geom->get_center_x();
0645 _y_debug = tower_geom->get_center_y();
0646 _z_debug = tower_geom->get_center_z();
0647
0648 PHG4Particle* primary = towereval->max_truth_primary_particle_by_energy(tower);
0649
0650 float gparticleID = NAN;
0651 float gflavor = NAN;
0652 float gnhits = NAN;
0653 float gpx = NAN;
0654 float gpy = NAN;
0655 float gpz = NAN;
0656 float ge = NAN;
0657
0658 float gpt = NAN;
0659 float geta = NAN;
0660 float gphi = NAN;
0661
0662 float gvx = NAN;
0663 float gvy = NAN;
0664 float gvz = NAN;
0665
0666 float gembed = NAN;
0667 float gedep = NAN;
0668
0669 float efromtruth = NAN;
0670
0671 if (primary)
0672 {
0673 gparticleID = primary->get_track_id();
0674 gflavor = primary->get_pid();
0675
0676 PHG4Shower* shower = trutheval->get_primary_shower(primary);
0677 if (shower)
0678 gnhits = shower->get_nhits(trutheval->get_caloid());
0679 else
0680 gnhits = 0.0;
0681 gpx = primary->get_px();
0682 gpy = primary->get_py();
0683 gpz = primary->get_pz();
0684 ge = primary->get_e();
0685
0686 gpt = sqrt(gpx * gpx + gpy * gpy);
0687 if (gpt != 0.0) geta = asinh(gpz / gpt);
0688 gphi = atan2(gpy, gpx);
0689
0690 PHG4VtxPoint* vtx = trutheval->get_vertex(primary);
0691
0692 if (vtx)
0693 {
0694 gvx = vtx->get_x();
0695 gvy = vtx->get_y();
0696 gvz = vtx->get_z();
0697 }
0698
0699 gembed = trutheval->get_embed(primary);
0700 gedep = trutheval->get_shower_energy_deposit(primary);
0701
0702 efromtruth = towereval->get_energy_contribution(tower, primary);
0703 }
0704
0705 float tower_data[] = {(float) _ievent,
0706 towerid,
0707 ieta,
0708 iphi,
0709 eta,
0710 phi,
0711 e,
0712 x,
0713 y,
0714 z,
0715 gparticleID,
0716 gflavor,
0717 gnhits,
0718 geta,
0719 gphi,
0720 ge,
0721 gpt,
0722 gvx,
0723 gvy,
0724 gvz,
0725 gembed,
0726 gedep,
0727 efromtruth};
0728
0729 _ntp_tower->Fill(tower_data);
0730 _tower_debug->Fill();
0731 }
0732 }
0733
0734
0735
0736
0737
0738 if (_do_cluster_eval)
0739 {
0740 if (Verbosity() > 1) cout << "CaloEvaluator::filling gcluster ntuple..." << endl;
0741
0742 GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0743
0744 string clusternode = "CLUSTER_" + _caloname;
0745 RawClusterContainer* clusters = findNode::getClass<RawClusterContainer>(topNode, clusternode.c_str());
0746 if (!clusters)
0747 {
0748 cerr << PHWHERE << " ERROR: Can't find " << clusternode << endl;
0749 exit(-1);
0750 }
0751
0752 string towernode = "TOWER_CALIB_" + _caloname;
0753 RawTowerContainer* towers = findNode::getClass<RawTowerContainer>(topNode, towernode.c_str());
0754 if (!towers) {
0755 cerr << PHWHERE << " ERROR: Can't find towers " << towernode << endl;
0756 exit(-1);
0757 }
0758
0759
0760
0761 for (const auto& iterator : clusters->getClustersMap())
0762 {
0763 RawCluster* cluster = iterator.second;
0764
0765
0766
0767
0768
0769 if (cluster->get_energy() < _reco_e_threshold) continue;
0770
0771 float clusterID = cluster->get_id();
0772 float ntowers = cluster->getNTowers();
0773 float x = cluster->get_x();
0774 float y = cluster->get_y();
0775 float z = cluster->get_z();
0776 float eta = NAN;
0777 float eta_detector = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(0,0,0));
0778 float phi = cluster->get_phi();
0779 float e = cluster->get_energy();
0780 float ecore = cluster->get_ecore();
0781
0782
0783 if (vertexmap)
0784 {
0785 if (!vertexmap->empty())
0786 {
0787 GlobalVertex* vertex = (vertexmap->begin()->second);
0788
0789 eta = RawClusterUtility::GetPseudorapidity(*cluster, CLHEP::Hep3Vector(vertex->get_x(), vertex->get_y(), vertex->get_z()));
0790 }
0791 }
0792
0793 PHG4Particle* primary = clustereval->max_truth_primary_particle_by_energy(cluster);
0794
0795 float gparticleID = NAN;
0796 float gflavor = NAN;
0797
0798 float gnhits = NAN;
0799 float gpx = NAN;
0800 float gpy = NAN;
0801 float gpz = NAN;
0802 float ge = NAN;
0803
0804 float gpt = NAN;
0805 float geta = NAN;
0806 float gphi = NAN;
0807
0808 float gvx = NAN;
0809 float gvy = NAN;
0810 float gvz = NAN;
0811
0812 float gembed = NAN;
0813 float gedep = NAN;
0814
0815 float efromtruth = NAN;
0816
0817 if (primary)
0818 {
0819 gparticleID = primary->get_track_id();
0820 gflavor = primary->get_pid();
0821
0822 PHG4Shower* shower = trutheval->get_primary_shower(primary);
0823 if (shower)
0824 gnhits = shower->get_nhits(trutheval->get_caloid());
0825 else
0826 gnhits = 0.0;
0827 gpx = primary->get_px();
0828 gpy = primary->get_py();
0829 gpz = primary->get_pz();
0830 ge = primary->get_e();
0831
0832 gpt = sqrt(gpx * gpx + gpy * gpy);
0833 if (gpt != 0.0) geta = asinh(gpz / gpt);
0834 gphi = atan2(gpy, gpx);
0835
0836 PHG4VtxPoint* vtx = trutheval->get_vertex(primary);
0837
0838 if (vtx)
0839 {
0840 gvx = vtx->get_x();
0841 gvy = vtx->get_y();
0842 gvz = vtx->get_z();
0843 }
0844
0845 gembed = trutheval->get_embed(primary);
0846 gedep = trutheval->get_shower_energy_deposit(primary);
0847
0848 efromtruth = clustereval->get_energy_contribution(cluster,
0849 primary);
0850 }
0851
0852 float cluster_data[] = {(float) _ievent,
0853 clusterID,
0854 ntowers,
0855 eta_detector,
0856 eta,
0857 x,
0858 y,
0859 z,
0860 phi,
0861 e,
0862 ecore,
0863 gparticleID,
0864 gflavor,
0865 gnhits,
0866 geta,
0867 gphi,
0868 ge,
0869 gpt,
0870 gvx,
0871 gvy,
0872 gvz,
0873 gembed,
0874 gedep,
0875 efromtruth};
0876
0877
0878 RawCluster::TowerConstRange towersConstRange = cluster->get_towers();
0879 RawCluster::TowerConstIterator toweriter;
0880
0881 _towerEtas.clear();
0882 _towerPhis.clear();
0883 _towerEnergies.clear();
0884
0885 for (toweriter = towersConstRange.first; toweriter != towersConstRange.second; ++toweriter) {
0886 RawTower* tower = towers->getTower(toweriter->first);
0887
0888 unsigned char towerEta = tower->get_bineta();
0889 unsigned char towerPhi = tower->get_binphi();
0890 float towerEnergy = tower->get_energy();
0891
0892
0893 _towerEtas.push_back(towerEta);
0894 _towerPhis.push_back(towerPhi);
0895 _towerEnergies.push_back(towerEnergy);
0896 }
0897 _cluster_tower_info->Fill();
0898 _ntp_cluster->Fill(cluster_data);
0899 }
0900 }
0901 return;
0902 }
