ParticleGen-Jingyu/physiTuto/CaloOnly.cc

0001 #include "CaloOnly.h"
0002
0003 #include <calobase/RawTowerContainer.h>
0004 #include <calobase/RawTower.h>
0005 #include <calobase/RawClusterContainer.h>
0006 #include <calobase/RawTowerGeomContainer.h>
0007 #include <calobase/RawCluster.h>
0008 #include <calobase/RawClusterUtility.h>
0009 #include <calobase/RawTowerDefs.h>
0010 #include <calobase/RawTowerGeom.h>
0011 #include <calobase/TowerInfoContainer.h>
0012 #include <calobase/TowerInfo.h>
0013 #include <calobase/TowerInfoDefs.h>
0014
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016
0017 #include <globalvertex/GlobalVertex.h>
0018 #include <globalvertex/GlobalVertexMap.h>
0019 #include <globalvertex/SvtxVertexMap.h>
0020 #include <globalvertex/SvtxVertex.h>
0021
0022 #include <phool/getClass.h>
0023 #include <phool/PHCompositeNode.h>
0024
0025 #include <trackbase/TrkrCluster.h>
0026 #include <trackbase/TrkrClusterContainer.h>
0027 #include <trackbase/TrkrDefs.h>
0028 #include <trackbase_historic/SvtxTrackMap.h>
0029 #include <trackbase_historic/SvtxTrackState_v1.h>
0030 #include <trackbase_historic/TrackSeedContainer.h>
0031 #include <trackbase_historic/TrackSeed.h>
0032 #include <trackreco/ActsPropagator.h>
0033
0034 #include <Acts/Geometry/GeometryIdentifier.hpp>
0035 #include <Acts/MagneticField/ConstantBField.hpp>
0036 #include <Acts/MagneticField/MagneticFieldProvider.hpp>
0037 #include <Acts/Surfaces/CylinderSurface.hpp>
0038 #include <Acts/Surfaces/PerigeeSurface.hpp>
0039 #include <Acts/Geometry/TrackingGeometry.hpp>
0040
0041 #include <CLHEP/Vector/ThreeVector.h>
0042 #include <math.h>
0043 #include <vector>
0044
0045 #include <TFile.h>
0046 #include <TTree.h>
0047 #include <TH1F.h>
0048 #include <TH2F.h>
0049
0050 //____________________________________________________________________________..
0051 CaloOnly::CaloOnly(const std::string &name, const std::string &file):
0052     SubsysReco(name),
0053     _outfilename(file),
0054     _outfile(nullptr),
0055     _tree(nullptr)
0056 {
0057     std::cout << "CaloOnly::CaloOnly(const std::string &name) Calling ctor" << std::endl;
0058 }
0059
0060 //____________________________________________________________________________..
0061 CaloOnly::~CaloOnly()
0062 {
0063     std::cout << "CaloOnly::~CaloOnly() Calling dtor" << std::endl;
0064 }
0065
0066 //____________________________________________________________________________..
0067 int CaloOnly::Init(PHCompositeNode *topNode)
0068 {
0069     std::cout << topNode << std::endl;
0070     std::cout << "CaloOnly::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0071     _outfile = new TFile(_outfilename.c_str(), "RECREATE");
0072     delete _tree;
0073
0074     _tree = new TTree("tree", "A tree with track/calo info");
0075     _tree->Branch("_run_test", &_run_test);
0076     // EMCal variables
0077     _tree->Branch("_emcalgeo_id", &_emcalgeo_id);
0078     _tree->Branch("_emcalgeo_phibin", &_emcalgeo_phibin);
0079     _tree->Branch("_emcalgeo_etabin", &_emcalgeo_etabin);
0080
0081     _tree->Branch("_emcal_e", &_emcal_e);
0082     _tree->Branch("_emcal_phi", &_emcal_phi);
0083     _tree->Branch("_emcal_eta", &_emcal_eta);
0084     _tree->Branch("_emcal_iphi", &_emcal_iphi);
0085     _tree->Branch("_emcal_ieta", &_emcal_ieta);
0086     _tree->Branch("_emcal_time", &_emcal_time);
0087     _tree->Branch("_emcal_chi2", &_emcal_chi2);
0088     _tree->Branch("_emcal_pedestal", &_emcal_pedestal);
0089     // IHCAL variables
0090     _tree->Branch("_ihcal_e", &_ihcal_e);
0091     _tree->Branch("_ihcal_phi", &_ihcal_phi);
0092     _tree->Branch("_ihcal_eta", &_ihcal_eta);
0093     _tree->Branch("_ihcal_iphi", &_ihcal_iphi);
0094     _tree->Branch("_ihcal_ieta", &_ihcal_ieta);
0095     _tree->Branch("_ihcal_time", &_ihcal_time);
0096     _tree->Branch("_ihcal_chi2", &_ihcal_chi2);
0097     _tree->Branch("_ihcal_pedestal", &_ihcal_pedestal);
0098     // OHCAL variables
0099     _tree->Branch("_ohcal_e", &_ohcal_e);
0100     _tree->Branch("_ohcal_phi", &_ohcal_phi);
0101     _tree->Branch("_ohcal_eta", &_ohcal_eta);
0102     _tree->Branch("_ohcal_iphi", &_ohcal_iphi);
0103     _tree->Branch("_ohcal_ieta", &_ohcal_ieta);
0104     _tree->Branch("_ohcal_time", &_ohcal_time);
0105     _tree->Branch("_ohcal_chi2", &_ohcal_chi2);
0106     _tree->Branch("_ohcal_pedestal", &_ohcal_pedestal);
0107
0108     // EMCal cluster information
0109     _tree->Branch("_emcal_cluster_id", &_emcal_cluster_id);
0110     _tree->Branch("_emcal_cluster_e", &_emcal_cluster_e);
0111     _tree->Branch("_emcal_cluster_phi", &_emcal_cluster_phi);
0112     _tree->Branch("_emcal_cluster_eta", &_emcal_cluster_eta);
0113     _tree->Branch("_emcal_cluster_x", &_emcal_cluster_x);
0114     _tree->Branch("_emcal_cluster_y", &_emcal_cluster_y);
0115     _tree->Branch("_emcal_cluster_z", &_emcal_cluster_z);
0116     _tree->Branch("_emcal_cluster_R", &_emcal_cluster_R);
0117     _tree->Branch("_emcal_cluster_ecore", &_emcal_cluster_ecore);
0118     _tree->Branch("_emcal_cluster_chi2", &_emcal_cluster_chi2);
0119     _tree->Branch("_emcal_cluster_prob", &_emcal_cluster_prob);
0120     // EMCal full cluster corrections
0121     _tree->Branch("_emcal_clusfull_e", &_emcal_clusfull_e);
0122     _tree->Branch("_emcal_clusfull_eta", &_emcal_clusfull_eta);
0123     _tree->Branch("_emcal_clusfull_phi", &_emcal_clusfull_phi);
0124     _tree->Branch("_emcal_clusfull_x", &_emcal_clusfull_x);
0125     _tree->Branch("_emcal_clusfull_y", &_emcal_clusfull_y);
0126     _tree->Branch("_emcal_clusfull_z", &_emcal_clusfull_z);
0127     _tree->Branch("_emcal_clusfull_R", &_emcal_clusfull_R);
0128     _tree->Branch("_emcal_clusfull_pt", &_emcal_clusfull_pt);
0129     // EMCal core cluster corrections
0130     _tree->Branch("_emcal_cluscore_e", &_emcal_cluscore_e);
0131     _tree->Branch("_emcal_cluscore_eta", &_emcal_cluscore_eta);
0132     _tree->Branch("_emcal_cluscore_phi", &_emcal_cluscore_phi);
0133     _tree->Branch("_emcal_cluscore_x", &_emcal_cluscore_x);
0134     _tree->Branch("_emcal_cluscore_y", &_emcal_cluscore_y);
0135     _tree->Branch("_emcal_cluscore_z", &_emcal_cluscore_z);
0136     _tree->Branch("_emcal_cluscore_R", &_emcal_cluscore_R);
0137     _tree->Branch("_emcal_cluscore_pt", &_emcal_cluscore_pt);
0138
0139     _tree->Branch("_mbd_z", &_mbd_z);
0140
0141     return Fun4AllReturnCodes::EVENT_OK;
0142 }
0143
0144 //____________________________________________________________________________..
0145 int CaloOnly::process_event(PHCompositeNode *topNode)
0146 {
0147     bool hasMBDvertex = true;
0148
0149     GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0150     if (!vertexmap)
0151     {
0152         std::cout << "CaloOnly::process_event - Fatal Error - GlobalVertexMap node is missing. Please turn on the do_global flag in the main  macro in order to reconstruct the global vertex." << std::endl;
0153         //assert(vertexmap);  // force quit
0154
0155         hasMBDvertex = false;
0156     }
0157
0158     if (vertexmap->empty())
0159     {
0160         std::cout << "CaloOnly::process_event - Fatal Error - GlobalVertexMap node is empty. Please turn on the do_global flag in the main macro  in order to reconstruct the global vertex." << std::endl;
0161
0162         hasMBDvertex = false;
0163     }
0164
0165     GlobalVertex *vtx = vertexmap->begin()->second;
0166     if (vtx == nullptr)
0167     {
0168         hasMBDvertex = false;
0169     }
0170
0171     if(hasMBDvertex == true)
0172     {
0173         _mbd_z.push_back(vtx->get_z());
0174     }
0175
0176     // ---------------------------- CALO info ---------------------------------------------------------------------
0177     // ---------------------------- CALO info ---------------------------------------------------------------------
0178     // ---------------------------- CALO info ---------------------------------------------------------------------
0179
0180     // 获取 RawTowerGeomContainer
0181     towerGeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0182     if (!towerGeom)
0183     {
0184         std::cout << "CaloOnly::process_event: Cannot find CEMC RawTowerGeomContainer!" << std::endl;
0185     }
0186
0187     // calo cluster
0188     clustersEM = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_CEMC");
0189     if (!clustersEM)
0190     {
0191         std::cout << "CaloOnly::process_event: cannot find cluster container " << "CLUSTERINFO_CEMC" << std::endl;
0192     }
0193
0194     // clustersIH = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_HCALIN");
0195     // if (!clustersIH)
0196     // {
0197     //     std::cout << "CaloOnly::process_event: cannot find cluster container " << "CLUSTERINFO_HCALIN" << std::endl;
0198     // }
0199
0200     // clustersOH = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_HCALOUT");
0201     // if (!clustersOH)
0202     // {
0203     //     std::cout << "CaloOnly::process_event: cannot find cluster container " << "CLUSTERINFO_HCALOUT" << std::endl;
0204     // }
0205
0206
0207     // calo tower
0208     EMCAL_Container = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0209     if(!EMCAL_Container)
0210     {
0211         std::cout << "CaloOnly::process_event: TOWERINFO_CALIB_CEMC not found!!!" << std::endl;
0212     }
0213
0214     IHCAL_Container = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0215     if(!IHCAL_Container)
0216     {
0217         std::cout << "CaloOnly::process_event: TOWERINFO_CALIB_HCALIN not found! Aborting!" << std::endl;
0218         return Fun4AllReturnCodes::ABORTEVENT;
0219     }
0220
0221     OHCAL_Container = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0222     if(!OHCAL_Container)
0223     {
0224         std::cout << "OHCAL_Container not found! Aborting!" << std::endl;
0225         return Fun4AllReturnCodes::ABORTEVENT;
0226     }
0227
0228
0229     // calo geometry
0230     EMCalGeo = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0231     if(!EMCalGeo)
0232     {
0233         std::cout << "CaloOnly::process_event: TOWERGEOM_CEMC not found!!!" << std::endl;
0234     }
0235
0236     IHCalGeo = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0237     if(!IHCalGeo)
0238     {
0239         std::cout << "CaloOnly::process_event: TOWERGEOM_HCALIN not found!!!" << std::endl;
0240     }
0241
0242     OHCalGeo = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0243     if(!OHCalGeo)
0244     {
0245         std::cout << "CaloOnly::process_event: TOWERGEOM_HCALOUT not found!!!" << std::endl;
0246     }
0247
0248     // // Rawtowercontainer
0249     // EMCal_RawTowerContainer = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
0250     // if(!EMCal_RawTowerContainer)
0251     // {
0252     //     std::cout << "CaloOnly::process_event: TOWER_CALIB_CEMC not found!!!" << std::endl;
0253     // }
0254
0255     ResetTreeVectors();
0256
0257     _run_test.push_back(1);
0258     FillTree();
0259
0260     // _tree->Fill();
0261
0262     return Fun4AllReturnCodes::EVENT_OK;
0263 }
0264
0265 //____________________________________________________________________________..
0266 int CaloOnly::End(PHCompositeNode *topNode)
0267 {
0268     std::cout << topNode << std::endl;
0269     _outfile->cd();
0270     _outfile->Write();
0271     _outfile->Close();
0272     return Fun4AllReturnCodes::EVENT_OK;
0273 }
0274
0275 void CaloOnly::ResetTreeVectors()
0276 {
0277     _run_test.clear();
0278     // EMCal vectors
0279     _emcalgeo_id.clear();
0280     _emcalgeo_phibin.clear();
0281     _emcalgeo_etabin.clear();
0282
0283     _emcal_e.clear();
0284     _emcal_phi.clear();
0285     _emcal_eta.clear();
0286     _emcal_iphi.clear();
0287     _emcal_ieta.clear();
0288     _emcal_time.clear();
0289     _emcal_chi2.clear();
0290     _emcal_pedestal.clear();
0291
0292     // IHCAL vectors
0293     _ihcal_e.clear();
0294     _ihcal_phi.clear();
0295     _ihcal_eta.clear();
0296     _ihcal_iphi.clear();
0297     _ihcal_ieta.clear();
0298     _ihcal_time.clear();
0299     _ihcal_chi2.clear();
0300     _ihcal_pedestal.clear();
0301
0302     // OHCAL vectors
0303     _ohcal_e.clear();
0304     _ohcal_phi.clear();
0305     _ohcal_eta.clear();
0306     _ohcal_iphi.clear();
0307     _ohcal_ieta.clear();
0308     _ohcal_time.clear();
0309     _ohcal_chi2.clear();
0310     _ohcal_pedestal.clear();
0311
0312     // EMCal cluster information
0313     _emcal_cluster_id.clear();
0314     _emcal_cluster_e.clear();
0315     _emcal_cluster_phi.clear();
0316     _emcal_cluster_eta.clear();
0317     _emcal_cluster_x.clear();
0318     _emcal_cluster_y.clear();
0319     _emcal_cluster_z.clear();
0320     _emcal_cluster_ecore.clear();
0321     _emcal_cluster_chi2.clear();
0322     _emcal_cluster_prob.clear();
0323
0324     // EMCal full cluster corrections
0325     _emcal_clusfull_e.clear();
0326     _emcal_clusfull_eta.clear();
0327     _emcal_clusfull_phi.clear();
0328     _emcal_clusfull_x.clear();
0329     _emcal_clusfull_y.clear();
0330     _emcal_clusfull_z.clear();
0331     _emcal_clusfull_pt.clear();
0332
0333     // EMCal core cluster corrections
0334     _emcal_cluscore_e.clear();
0335     _emcal_cluscore_eta.clear();
0336     _emcal_cluscore_phi.clear();
0337     _emcal_cluscore_x.clear();
0338     _emcal_cluscore_y.clear();
0339     _emcal_cluscore_z.clear();
0340     _emcal_cluscore_pt.clear();
0341 }
0342
0343
0344 // //____________________________________________________________________________..
0345 void CaloOnly::FillTree()
0346 {
0347     if (!clustersEM || !EMCAL_Container || !EMCalGeo || !IHCalGeo || !OHCalGeo || !IHCAL_Container || !OHCAL_Container )
0348     {
0349         std::cout << PHWHERE << "missing node trees, can't continue with track calo matching"
0350               << std::endl;
0351         return;
0352     }
0353
0354     // 遍历所有塔
0355     for (unsigned int i = 0; i < EMCalGeo->size(); i++)
0356     {
0357         RawTowerGeom* geom = EMCalGeo->get_tower_geometry(i);
0358         if (!geom) continue;
0359         std::cout<<"rawtower id is: "<<geom->get_id()<<std::endl;
0360
0361         std::cout << "Tower ID: " << i  // 这里的 first 是 tower id
0362                   << std::endl;
0363     }
0364
0365     // loop over all calo tower
0366     TowerInfo *tInfo = nullptr;
0367
0368     for(unsigned int iem = 0; iem < EMCAL_Container->size(); iem++)
0369     {
0370         tInfo = EMCAL_Container->get_tower_at_channel(iem);
0371
0372         if(!tInfo)
0373         {
0374             continue;
0375         }
0376
0377         if(!tInfo->get_isGood())
0378         {
0379             continue;
0380         }
0381
0382         //if(tInfo->get_energy() < 0.2) continue;
0383
0384         unsigned int towerinfo_key = EMCAL_Container->encode_key(iem);
0385         int ti_ieta = EMCAL_Container->getTowerEtaBin(towerinfo_key);
0386         int ti_iphi = EMCAL_Container->getTowerPhiBin(towerinfo_key);
0387
0388         const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ti_ieta, ti_iphi);
0389         RawTowerGeom *tower_geom = EMCalGeo->get_tower_geometry(key);
0390
0391         // _emcalgeo_id.push_back(tower_geom->get_id());
0392         // _emcalgeo_phibin.push_back(ti_iphi);
0393         // _emcalgeo_etabin.push_back(ti_ieta);
0394         // std::cout<<"rawtower id is: "<<tower_geom->get_id()<<std::endl;
0395
0396         _emcal_e.push_back(tInfo->get_energy());
0397         _emcal_phi.push_back(tower_geom->get_phi());
0398         _emcal_eta.push_back(tower_geom->get_eta());
0399         _emcal_iphi.push_back(ti_iphi);
0400         _emcal_ieta.push_back(ti_ieta);
0401         _emcal_time.push_back(tInfo->get_time());
0402         _emcal_chi2.push_back(tInfo->get_chi2());
0403         _emcal_pedestal.push_back(tInfo->get_pedestal());
0404     }
0405
0406     for(unsigned int ihcal = 0; ihcal < IHCAL_Container->size(); ihcal++)
0407     {
0408         tInfo = IHCAL_Container->get_tower_at_channel(ihcal);
0409
0410         if(!tInfo)
0411         {
0412             continue;
0413         }
0414
0415         if(!tInfo->get_isGood())
0416         {
0417             continue;
0418         }
0419
0420         //if(tInfo->get_energy() < 0.2) continue;
0421
0422         unsigned int towerinfo_key = IHCAL_Container->encode_key(ihcal);
0423         int ti_ieta = IHCAL_Container->getTowerEtaBin(towerinfo_key);
0424         int ti_iphi = IHCAL_Container->getTowerPhiBin(towerinfo_key);
0425
0426         const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ti_ieta, ti_iphi);
0427         RawTowerGeom *tower_geom = IHCalGeo->get_tower_geometry(key);
0428
0429         _ihcal_e.push_back(tInfo->get_energy());
0430         _ihcal_phi.push_back(tower_geom->get_phi());
0431         _ihcal_eta.push_back(tower_geom->get_eta());
0432         _ihcal_iphi.push_back(ti_iphi);
0433         _ihcal_ieta.push_back(ti_ieta);
0434         _ihcal_time.push_back(tInfo->get_time());
0435         _ihcal_chi2.push_back(tInfo->get_chi2());
0436         _ihcal_pedestal.push_back(tInfo->get_pedestal());
0437     }
0438
0439     for(unsigned int ohcal = 0; ohcal < OHCAL_Container->size(); ohcal++)
0440     {
0441         tInfo = OHCAL_Container->get_tower_at_channel(ohcal);
0442
0443         if(!tInfo)
0444         {
0445             continue;
0446         }
0447
0448         if(!tInfo->get_isGood())
0449         {
0450             continue;
0451         }
0452
0453         //if(tInfo->get_energy() < 0.2) continue;
0454
0455         unsigned int towerinfo_key = OHCAL_Container->encode_key(ohcal);
0456         int ti_ieta = OHCAL_Container->getTowerEtaBin(towerinfo_key);
0457         int ti_iphi = OHCAL_Container->getTowerPhiBin(towerinfo_key);
0458
0459         const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ti_ieta, ti_iphi);
0460         RawTowerGeom *tower_geom = OHCalGeo->get_tower_geometry(key);
0461
0462         _ohcal_e.push_back(tInfo->get_energy());
0463         _ohcal_phi.push_back(tower_geom->get_phi());
0464         _ohcal_eta.push_back(tower_geom->get_eta());
0465         _ohcal_iphi.push_back(ti_iphi);
0466         _ohcal_ieta.push_back(ti_ieta);
0467         _ohcal_time.push_back(tInfo->get_time());
0468         _ohcal_chi2.push_back(tInfo->get_chi2());
0469         _ohcal_pedestal.push_back(tInfo->get_pedestal());
0470     }
0471
0472     // Loop over the EMCal clusters
0473     CLHEP::Hep3Vector vertex(0., 0., 0.);
0474
0475     RawCluster *cluster = nullptr;
0476
0477     RawClusterContainer::Range begin_end_EMC = clustersEM->getClusters();
0478     RawClusterContainer::Iterator clusIter_EMC;
0479     for (clusIter_EMC = begin_end_EMC.first; clusIter_EMC != begin_end_EMC.second; ++clusIter_EMC)
0480     {
0481         cluster = clusIter_EMC->second;
0482         if(cluster->get_energy() < m_emcal_e_low_cut) continue;
0483
0484         // cluster info
0485         _emcal_cluster_id.push_back(clusIter_EMC->first);
0486         _emcal_cluster_e.push_back(cluster->get_energy());
0487         _emcal_cluster_phi.push_back(RawClusterUtility::GetAzimuthAngle(*cluster, vertex));
0488         _emcal_cluster_eta.push_back(RawClusterUtility::GetPseudorapidity(*cluster, vertex));
0489         _emcal_cluster_x.push_back(cluster->get_x());
0490         _emcal_cluster_y.push_back(cluster->get_y());
0491         _emcal_cluster_z.push_back(cluster->get_z());
0492         _emcal_cluster_ecore.push_back(cluster->get_ecore());
0493         _emcal_cluster_chi2.push_back(cluster->get_chi2());
0494         _emcal_cluster_prob.push_back(cluster->get_prob());
0495
0496         // after correction full cluster total energy info
0497         CLHEP::Hep3Vector E_vec_cluster_Full = RawClusterUtility::GetEVec(*cluster, vertex);
0498         float clusE = E_vec_cluster_Full.mag();
0499         float clus_eta = E_vec_cluster_Full.pseudoRapidity();
0500         float clus_phi = E_vec_cluster_Full.phi();
0501         float clus_x = E_vec_cluster_Full.x();
0502         float clus_y = E_vec_cluster_Full.y();
0503         float clus_z = E_vec_cluster_Full.z();
0504         float clus_pt = E_vec_cluster_Full.perp();
0505
0506         _emcal_clusfull_e.push_back(clusE);
0507         _emcal_clusfull_eta.push_back(clus_eta);
0508         _emcal_clusfull_phi.push_back(clus_phi);
0509         _emcal_clusfull_x.push_back(clus_x);
0510         _emcal_clusfull_y.push_back(clus_y);
0511         _emcal_clusfull_z.push_back(clus_z);
0512         _emcal_clusfull_pt.push_back(clus_pt);
0513
0514         // double phi_emcal_corr = clus_phi;
0515
0516         // after correction cluster core energy info
0517         CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*cluster, vertex);
0518         float clusEcore = E_vec_cluster.mag();
0519         float clus_coreeta = E_vec_cluster.pseudoRapidity();
0520         float clus_corephi = E_vec_cluster.phi();
0521         float clus_corex = E_vec_cluster.x();
0522         float clus_corey = E_vec_cluster.y();
0523         float clus_corez = E_vec_cluster.z();
0524         float clus_corept = E_vec_cluster.perp();
0525
0526         _emcal_cluscore_e.push_back(clusEcore);
0527         _emcal_cluscore_eta.push_back(clus_coreeta);
0528         _emcal_cluscore_phi.push_back(clus_corephi);
0529         _emcal_cluscore_x.push_back(clus_corex);
0530         _emcal_cluscore_y.push_back(clus_corey);
0531         _emcal_cluscore_z.push_back(clus_corez);
0532         _emcal_cluscore_pt.push_back(clus_corept);
0533
0534         // // zi ji suan ge cluster
0535         // RawCluster::TowerConstRange towers = cluster->get_towers();
0536         // RawCluster::TowerConstIterator toweriter;
0537
0538         // TowerInfo *towerInfo = nullptr;
0539
0540         // for (toweriter = towers.first; toweriter != towers.second; ++toweriter)
0541         // {
0542         //     _emcal_tower_cluster_id.push_back(clusIter_EMC->first);
0543         //     RawTowerGeom *tower_geom = EMCalGeo->get_tower_geometry(toweriter->first);
0544         //     _emcal_tower_phi.push_back(tower_geom->get_phi());
0545         //     _emcal_tower_eta.push_back(tower_geom->get_eta());
0546         //     _emcal_tower_e.push_back(toweriter->second);
0547         //     unsigned int key = TowerInfoDefs::encode_emcal(tower_geom->get_bineta(), tower_geom->get_binphi());
0548         //     towerInfo = EMCAL_Container->get_tower_at_key(key);
0549         //     _emcal_tower_status.push_back(towerInfo->get_status());
0550         // }
0551     }
0552
0553     // // loop over truth primary particles
0554     // PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
0555     // for (PHG4TruthInfoContainer::ConstIterator truth_itr = range.first; truth_itr != range.second; ++truth_itr)
0556     // {
0557     //     PHG4Particle *truth = truth_itr->second;
0558     //     _truth_px.push_back(truth->get_px());
0559     //     _truth_py.push_back(truth->get_py());
0560     //     _truth_pz.push_back(truth->get_pz());
0561     //     _truth_e.push_back(truth->get_e());
0562     //     _truth_pt.push_back(sqrt(truth->get_px() * truth->get_px() + truth->get_py() * truth->get_py()));
0563     //     _truth_eta.push_back(atanh(truth->get_pz() / sqrt(truth->get_px() * truth->get_px() + truth->get_py() * truth->get_py() + truth->get_pz() * truth->get_pz())));
0564     //     _truth_phi.push_back(atan2(truth->get_py(),truth->get_px()));
0565     //     _truth_pid.push_back(truth->get_pid());
0566     //     //std::cout<<"pid = "<<truth->get_pid()<<std::endl;
0567     // }
0568
0569     // // map between track id and charge deposit
0570     // std::map<int, float> trkmap;
0571
0572     // hits_CEMC = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_CEMC");
0573     // if(!hits_CEMC)
0574     // {
0575     //   std::cout << "CaloOnly::process_event: G4HIT_CEMC not found!!!" << std::endl;
0576     // }
0577
0578     // TLorentzVector v4;
0579     // PHG4HitContainer::ConstRange hit_range = hits_CEMC->getHits();
0580     // for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0581     // {
0582     //     PHG4Hit *this_hit = hit_iter->second;
0583     //     float light_yield = hit_iter->second->get_light_yield();
0584     //     float edep = hit_iter->second->get_edep();
0585     //     int ch = hit_iter->second->get_layer();
0586     //     float x = hit_iter->second->get_x(0);
0587     //     float y = hit_iter->second->get_y(0);
0588     //     float z = hit_iter->second->get_z(0);
0589     //     float t = hit_iter->second->get_t(0);
0590     //     int trkid = hit_iter->second->get_trkid();
0591     //     int truthtrkid = trkid;
0592     //     PHG4Particle *part = truthinfo->GetParticle(trkid);
0593
0594     //     v4.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e());
0595     //     int pid = part->get_pid();
0596     //     _CEMC_Hit_pid.push_back(pid);
0597
0598     //     int vtxid = part->get_vtx_id();
0599     //     PHG4VtxPoint *vtx = truthinfo->GetVtx(vtxid);
0600
0601     //     // add trkid to a set
0602     //     _CEMC_Hit_Evis.push_back(light_yield);
0603     //     _CEMC_Hit_Edep.push_back(edep);
0604     //     _CEMC_Hit_ch.push_back(ch);
0605     //     _CEMC_Hit_x.push_back(x);
0606     //     _CEMC_Hit_y.push_back(y);
0607     //     _CEMC_Hit_z.push_back(z);
0608     //     _CEMC_Hit_t.push_back(t);
0609     //     _CEMC_Hit_particle_x.push_back(vtx->get_x());
0610     //     _CEMC_Hit_particle_y.push_back(vtx->get_y());
0611     //     _CEMC_Hit_particle_z.push_back(vtx->get_z());
0612     // }
0613
0614     _tree->Fill();
0615 }
0616
0617
0618 // -----------------------------------------------------------------------------------------------------------------------
0619 // -----------------------------------------------------------------------------------------------------------------------
0620 // -----------------------------------------------------------------------------------------------------------------------
0621 // -----------------------------------------------------------------------------------------------------------------------
0622 // -----------------------------------------------------------------------------------------------------------------------
0623 // -----------------------------------------------------------------------------------------------------------------------
0624 // -----------------------------------------------------------------------------------------------------------------------
0625 // -----------------------------------------------------------------------------------------------------------------------
0626 // -----------------------------------------------------------------------------------------------------------------------
0627 // get some calo info:
0628 // RawTowerContainer - RawTower : Tower of calorimeters, records a tower ID# -> GeV of calibrated energy
0629 // findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
0630 // findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_EEMC");
0631 // findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
0632 // findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
0633 // findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
0634
0635 // findNode::getClass<RawTowerContainer>(topNode, "TOWER_RAW_SPILL_WARBLER");
0636 // findNode::getClass<RawTowerContainer>(topNode, "TOWER_RAW_CEMC");
0637 // findNode::getClass<RawTowerContainer>(topNode, "TOWER_RAW_LG_CEMC");
0638
0639 // findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER");
0640
0641
0642 // RawClusterContainer - RawCluster : Cluster as groups of calorimeter towers
0643 // findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
0644 // findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALIN");
0645 // findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALOUT");
0646
0647 // findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC_MOD");
0648 // findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALIN_MOD");
0649 // findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_HCALOUT_MOD");
0650
0651 // RawClusterContainer *EMCAL_RawClusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_POS_COR_CEMC");
0652 // RawClusterContainer *clustersEM = findNode::getClass<RawClusterContainer>(topNode, "TOPOCLUSTER_EMCAL");
0653 // RawClusterContainer *clustersHAD = findNode::getClass<RawClusterContainer>(topNode, "TOPOCLUSTER_HCAL");
0654
0655 // findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_" + m_detector);