Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:18:08 

0001 //____________________________________________________________________________..
0002 //
0003 // This is a template for a Fun4All SubsysReco module with all methods from the
0004 // $OFFLINE_MAIN/include/fun4all/SubsysReco.h baseclass
0005 // You do not have to implement all of them, you can just remove unused methods
0006 // here and in tutorial.h.
0007 //
0008 // tutorial(const std::string &name = "tutorial")
0009 // everything is keyed to tutorial, duplicate names do work but it makes
0010 // e.g. finding culprits in logs difficult or getting a pointer to the module
0011 // from the command line
0012 //
0013 // tutorial::~tutorial()
0014 // this is called when the Fun4AllServer is deleted at the end of running. Be
0015 // mindful what you delete - you do loose ownership of object you put on the node tree
0016 //
0017 // int tutorial::Init(PHCompositeNode *topNode)
0018 // This method is called when the module is registered with the Fun4AllServer. You
0019 // can create historgrams here or put objects on the node tree but be aware that
0020 // modules which haven't been registered yet did not put antyhing on the node tree
0021 //
0022 // int tutorial::InitRun(PHCompositeNode *topNode)
0023 // This method is called when the first event is read (or generated). At
0024 // this point the run number is known (which is mainly interesting for raw data
0025 // processing). Also all objects are on the node tree in case your module's action
0026 // depends on what else is around. Last chance to put nodes under the DST Node
0027 // We mix events during readback if branches are added after the first event
0028 //
0029 // int tutorial::process_event(PHCompositeNode *topNode)
0030 // called for every event. Return codes trigger actions, you find them in
0031 // $OFFLINE_MAIN/include/fun4all/Fun4AllReturnCodes.h
0032 //   everything is good:
0033 //     return Fun4AllReturnCodes::EVENT_OK
0034 //   abort event reconstruction, clear everything and process next event:
0035 //     return Fun4AllReturnCodes::ABORT_EVENT; 
0036 //   proceed but do not save this event in output (needs output manager setting):
0037 //     return Fun4AllReturnCodes::DISCARD_EVENT; 
0038 //   abort processing:
0039 //     return Fun4AllReturnCodes::ABORT_RUN
0040 // all other integers will lead to an error and abort of processing
0041 //
0042 // int tutorial::ResetEvent(PHCompositeNode *topNode)
0043 // If you have internal data structures (arrays, stl containers) which needs clearing
0044 // after each event, this is the place to do that. The nodes under the DST node are cleared
0045 // by the framework
0046 //
0047 // int tutorial::EndRun(const int runnumber)
0048 // This method is called at the end of a run when an event from a new run is
0049 // encountered. Useful when analyzing multiple runs (raw data). Also called at
0050 // the end of processing (before the End() method)
0051 //
0052 // int tutorial::End(PHCompositeNode *topNode)
0053 // This is called at the end of processing. It needs to be called by the macro
0054 // by Fun4AllServer::End(), so do not forget this in your macro
0055 //
0056 // int tutorial::Reset(PHCompositeNode *topNode)
0057 // not really used - it is called before the dtor is called
0058 //
0059 // void tutorial::Print(const std::string &what) const
0060 // Called from the command line - useful to print information when you need it
0061 //
0062 //____________________________________________________________________________..
0063 
0064 #include "tutorial.h"
0065 
0066 #include <globalvertex/GlobalVertexMap.h>
0067 #include <globalvertex/SvtxVertexMap.h>
0068 #include <ffarawobjects/Gl1Packet.h>
0069 #include <trackbase_historic/SvtxPHG4ParticleMap_v1.h>
0070 #include <trackbase_historic/SvtxTrack.h>
0071 #include <trackbase_historic/SvtxTrack_v1.h>
0072 #include <trackbase_historic/SvtxTrack_v2.h>
0073 #include <trackbase_historic/SvtxTrackMap.h>
0074 #include <trackbase_historic/SvtxTrackMap_v1.h>
0075 #include <trackbase_historic/SvtxTrackState_v1.h>
0076 #include <trackbase_historic/TrackSeed.h>
0077 
0078 #include <calobase/RawClusterContainer.h>
0079 #include <calobase/RawTowerGeomContainer.h>
0080 #include <calobase/RawCluster.h>
0081 #include <calobase/RawClusterUtility.h>
0082 #include <calobase/RawTowerDefs.h>
0083 #include <calobase/RawTowerGeom.h>
0084 #include <calobase/RawTowerGeomv5.h>
0085 #include <calobase/TowerInfoContainer.h>
0086 #include <calobase/TowerInfo.h>
0087 #include <calobase/TowerInfoDefs.h>
0088 
0089 #include <ffaobjects/EventHeaderv1.h>
0090 
0091 #include <fun4all/Fun4AllReturnCodes.h>
0092 
0093 #include <globalvertex/GlobalVertex.h>
0094 #include <globalvertex/GlobalVertexMap.h>
0095 #include <globalvertex/MbdVertexMap.h>
0096 #include <globalvertex/MbdVertex.h>
0097 #include <globalvertex/SvtxVertexMap.h>
0098 #include <globalvertex/SvtxVertex.h>
0099 
0100 #include <phool/getClass.h>
0101 #include <phool/PHCompositeNode.h>
0102 
0103 #include <Acts/Geometry/GeometryIdentifier.hpp>
0104 #include <Acts/MagneticField/ConstantBField.hpp>
0105 #include <Acts/MagneticField/MagneticFieldProvider.hpp>
0106 #include <Acts/Surfaces/CylinderSurface.hpp>
0107 #include <Acts/Surfaces/PerigeeSurface.hpp>
0108 #include <Acts/Geometry/TrackingGeometry.hpp>
0109 
0110 #include <CLHEP/Vector/ThreeVector.h>
0111 #include <math.h>
0112 #include <vector>
0113 
0114 #include <TFile.h>
0115 #include <TTree.h>
0116 #include <TH1F.h>
0117 #include <TH2F.h>
0118 #include <TLorentzVector.h>
0119 
0120 //____________________________________________________________________________..
0121 tutorial::tutorial(
0122     const std::string & name_in,
0123     std::string output_path_in,
0124     std::string output_rootfile_name_in):
0125     SubsysReco(name_in),
0126     output_path(output_path_in),
0127     output_rootfile_name(output_rootfile_name_in),
0128     file_out(nullptr),
0129     output( nullptr )
0130 {
0131   std::cout << "tutorial::tutorial(const std::string &name) Calling ctor" << std::endl;
0132 }
0133 
0134 //____________________________________________________________________________..
0135 tutorial::~tutorial()
0136 {
0137   std::cout << "tutorial::~tutorial() Calling dtor" << std::endl;
0138 }
0139 
0140 //____________________________________________________________________________..
0141 int tutorial::Init(PHCompositeNode *topNode)
0142 {
0143     std::cout << topNode << std::endl;
0144   std::cout << "tutorial::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0145 
0146   ////////////////////////////////////////////////////////
0147   // Initialization of the member                       //
0148   ////////////////////////////////////////////////////////
0149 //   output = new TFile( output_path.c_str(), "RECREATE" );
0150   
0151     std::cout << "tutorial::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0152 
0153     eventID = 0;
0154     NClus = 0;
0155     clus_system.clear();
0156     clus_layer.clear();
0157     clus_adc.clear();
0158     clus_X.clear();
0159     clus_Y.clear();
0160     clus_Z.clear();
0161     clus_size.clear();
0162     clus_phi_size.clear();
0163     clus_z_size.clear();
0164 
0165     nTowers = 0;
0166     tower_X.clear();
0167     tower_Y.clear();
0168     tower_Z.clear();
0169     tower_R.clear();
0170     tower_Eta.clear();
0171     tower_Phi.clear();
0172     tower_Eta_test.clear();
0173     tower_Phi_test.clear();
0174     tower_Eta_bin.clear();
0175     tower_Phi_bin.clear();
0176     tower_edep.clear();
0177     tower_system.clear();
0178 
0179     tower_int_X.clear();
0180     tower_int_Y.clear();
0181     tower_int_Z.clear();
0182     tower_int_R.clear();
0183 
0184     nCaloClus = 0;
0185     caloClus_X.clear();
0186     caloClus_Y.clear();
0187     caloClus_Z.clear();
0188     caloClus_R.clear();
0189     caloClus_Phi.clear();
0190     caloClus_edep.clear();
0191     caloClus_system.clear();
0192 
0193     caloClus_innr_X.clear();
0194     caloClus_innr_Y.clear();
0195     caloClus_innr_Z.clear();
0196     caloClus_innr_R.clear();
0197     caloClus_innr_Phi.clear();
0198     caloClus_innr_edep.clear();
0199 
0200     // note : Truth primary vertex information
0201     TruthPV_trig_x_ = -999;
0202     TruthPV_trig_y_ = -999;
0203     TruthPV_trig_z_ = -999;
0204     NTruthVtx_ = 0;
0205 
0206     // note : PHG4 information (from all PHG4Particles)
0207     NPrimaryG4P_ = 0;
0208     NPrimaryG4P_promptChargeHadron_ = 0;
0209     PrimaryG4P_Pt_.clear();
0210     PrimaryG4P_Eta_.clear();
0211     PrimaryG4P_Phi_.clear();
0212     PrimaryG4P_E_.clear();
0213     PrimaryG4P_PID_.clear();
0214     PrimaryG4P_ParticleClass_.clear();
0215     PrimaryG4P_isStable_.clear();
0216     PrimaryG4P_Charge_.clear();
0217     PrimaryG4P_isChargeHadron_.clear();
0218 
0219     _CEMC_Hit_Evis.clear();
0220     _CEMC_Hit_Edep.clear();
0221     _CEMC_Hit_ch.clear();
0222     _CEMC_Hit_x.clear();
0223     _CEMC_Hit_y.clear();
0224     _CEMC_Hit_z.clear();
0225 
0226     _CEMC_Pr_Hit_x.clear();
0227     _CEMC_Pr_Hit_y.clear();
0228     _CEMC_Pr_Hit_z.clear();
0229     _CEMC_Pr_Hit_R.clear();
0230     _CEMC_Pr_Hit_deltaT.clear();
0231 
0232     primary_electron_tracks.clear();
0233 
0234     ////////////////////////////////////////////////////////
0235     // Initialization of the member                       //
0236     ////////////////////////////////////////////////////////
0237     file_out = new TFile((output_path + "/" + output_rootfile_name).c_str(), "RECREATE");
0238     tree_out = new TTree("tree", "sPHENIX info.");
0239     
0240     // note : the tracker cluster
0241     tree_out -> Branch("trk_NClus", & NClus);
0242     tree_out -> Branch("trk_system", & clus_system);
0243     tree_out -> Branch("trk_layer", & clus_layer);
0244     tree_out -> Branch("trk_adc", & clus_adc);
0245     tree_out -> Branch("trk_X", & clus_X);
0246     tree_out -> Branch("trk_Y", & clus_Y);
0247     tree_out -> Branch("trk_Z", & clus_Z);
0248     tree_out -> Branch("trk_size", & clus_size);
0249     tree_out -> Branch("trk_phi_size", & clus_phi_size);
0250     tree_out -> Branch("trk_Z_size", & clus_z_size);
0251 
0252 
0253     // note : the calorimeter raw tower information with the calibration, and introduction of noise hits
0254     tree_out -> Branch("nTowers", & nTowers);
0255     tree_out -> Branch("tower_system", & tower_system);
0256     tree_out -> Branch("tower_X", & tower_X);
0257     tree_out -> Branch("tower_Y", & tower_Y);
0258     tree_out -> Branch("tower_Z", & tower_Z);
0259     tree_out -> Branch("tower_R", & tower_R);
0260     tree_out -> Branch("tower_Eta", & tower_Eta);
0261     tree_out -> Branch("tower_Phi", & tower_Phi);
0262     tree_out -> Branch("tower_Eta_test", & tower_Eta_test);
0263     tree_out -> Branch("tower_Phi_test", & tower_Phi_test);
0264     tree_out -> Branch("tower_Eta_bin", & tower_Eta_bin);
0265     tree_out -> Branch("tower_Phi_bin", & tower_Phi_bin);
0266     tree_out -> Branch("tower_edep", & tower_edep);
0267 
0268     tree_out -> Branch("tower_int_X", & tower_int_X);
0269     tree_out -> Branch("tower_int_Y", & tower_int_Y);
0270     tree_out -> Branch("tower_int_Z", & tower_int_Z);
0271     tree_out -> Branch("tower_int_R", & tower_int_R);
0272 
0273 
0274     // note : the calorimeter raw tower information with the calibration, and introduction of noise hits
0275     tree_out -> Branch("nCaloClus", & nCaloClus);
0276     tree_out -> Branch("caloClus_system", & caloClus_system);
0277     tree_out -> Branch("caloClus_X", & caloClus_X);
0278     tree_out -> Branch("caloClus_Y", & caloClus_Y);
0279     tree_out -> Branch("caloClus_Z", & caloClus_Z);
0280     tree_out -> Branch("caloClus_R", & caloClus_R);
0281     tree_out -> Branch("caloClus_Phi", & caloClus_Phi);
0282     tree_out -> Branch("caloClus_edep", & caloClus_edep);
0283 
0284     tree_out -> Branch("caloClus_innr_X", & caloClus_innr_X);
0285     tree_out -> Branch("caloClus_innr_Y", & caloClus_innr_Y);
0286     tree_out -> Branch("caloClus_innr_Z", & caloClus_innr_Z);
0287     tree_out -> Branch("caloClus_innr_R", & caloClus_innr_R);
0288     tree_out -> Branch("caloClus_innr_Phi", & caloClus_innr_Phi);
0289     tree_out -> Branch("caloClus_innr_edep", & caloClus_innr_edep);
0290 
0291     // note : true vertex information from G4Particle 
0292     tree_out->Branch("NTruthVtx", &NTruthVtx_);
0293     tree_out->Branch("TruthPV_trig_x", &TruthPV_trig_x_);
0294     tree_out->Branch("TruthPV_trig_y", &TruthPV_trig_y_);
0295     tree_out->Branch("TruthPV_trig_z", &TruthPV_trig_z_);
0296     // note : Primary PHG4Particle information
0297     tree_out->Branch("NPrimaryG4P", &NPrimaryG4P_);
0298     tree_out->Branch("PrimaryG4P_Pt", &PrimaryG4P_Pt_);
0299     tree_out->Branch("PrimaryG4P_Eta", &PrimaryG4P_Eta_);
0300     tree_out->Branch("PrimaryG4P_Phi", &PrimaryG4P_Phi_);
0301     tree_out->Branch("PrimaryG4P_E", &PrimaryG4P_E_);
0302     tree_out->Branch("PrimaryG4P_PID", &PrimaryG4P_PID_);
0303     tree_out->Branch("PrimaryG4P_isChargeHadron", &PrimaryG4P_isChargeHadron_);
0304 
0305     // truthinfo G4hit info
0306     tree_out->Branch("CEMC_Hit_Evis", &_CEMC_Hit_Evis);
0307     tree_out->Branch("CEMC_Hit_Edep", &_CEMC_Hit_Edep);
0308     tree_out->Branch("CEMC_Hit_ch", &_CEMC_Hit_ch); 
0309     tree_out->Branch("CEMC_Hit_x", &_CEMC_Hit_x);
0310     tree_out->Branch("CEMC_Hit_y", &_CEMC_Hit_y);
0311     tree_out->Branch("CEMC_Hit_z", &_CEMC_Hit_z);
0312 
0313     // Primary particle truth hit on cemc 
0314     tree_out->Branch("CEMC_Pr_Hit_x", &_CEMC_Pr_Hit_x);
0315     tree_out->Branch("CEMC_Pr_Hit_y", &_CEMC_Pr_Hit_y);
0316     tree_out->Branch("CEMC_Pr_Hit_z", &_CEMC_Pr_Hit_z);
0317     tree_out->Branch("CEMC_Pr_Hit_R", &_CEMC_Pr_Hit_R);
0318     tree_out->Branch("_CEMC_Pr_Hit_deltaT", &_CEMC_Pr_Hit_deltaT);
0319 
0320     tree_out->Branch("PrG4_TTPRO_dD", &_PrG4_TTPRO_dD);
0321     tree_out->Branch("PrG4_TTPRO_dR", &_PrG4_TTPRO_dR);
0322     tree_out->Branch("PrG4_TTPRO_dphi", &_PrG4_TTPRO_dphi);
0323 
0324     return Fun4AllReturnCodes::EVENT_OK;
0325 }
0326 
0327 //____________________________________________________________________________..
0328 int tutorial::InitRun(PHCompositeNode * /*topNode*/)
0329 {
0330   std::cout << "tutorial::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0331   return Fun4AllReturnCodes::EVENT_OK;
0332 }
0333 
0334 
0335 
0336 //____________________________________________________________________________..
0337 int tutorial::prepareTracker(PHCompositeNode * topNode) 
0338 {
0339   std::cout << "tutorial::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0340 
0341   //////////////////////////////////////////////////////////////////////////////
0342   // Getting Nodes                                                            //
0343   /////////////////////////////////////////////////////////////////////////////  
0344   // TRKR_CLUSTER node: Information of TrkrCluster
0345   node_cluster_map = findNode::getClass<TrkrClusterContainerv4>(topNode, "TRKR_CLUSTER");
0346 
0347   if (!node_cluster_map)
0348     {
0349       std::cerr << PHWHERE << "TrkrClusterContainer node is missing." << std::endl;
0350       return Fun4AllReturnCodes::ABORTEVENT;
0351     }
0352 
0353   // ActsGeometry node: for the global coordinate
0354   node_acts = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0355   if ( !node_acts )
0356     {
0357       std::cout << PHWHERE << "No ActsGeometry on node tree. Bailing." << std::endl;
0358       return Fun4AllReturnCodes::ABORTEVENT;
0359     }
0360 
0361   
0362     std::vector <TrkrCluster *> mvtx_clusters; mvtx_clusters.clear();
0363     std::vector <TrkrCluster *> intt_clusters; intt_clusters.clear();
0364     std::vector <TrkrCluster *> tpc_clusters; tpc_clusters.clear();
0365 
0366     // note : mvtx
0367     for (unsigned int mvtxlayer = 0; mvtxlayer < 3; mvtxlayer++) {
0368         //      cout << " INTT layer " << mvtxlayer << endl;
0369         //      int layer= ( mvtxlayer < 2 ? 0 : 1 );
0370 
0371         // loop over all hits
0372         for (const auto & hitsetkey: node_cluster_map -> getHitSetKeys(TrkrDefs::TrkrId::mvtxId, mvtxlayer)) {
0373 
0374             // type: std::pair<ConstIterator, ConstIterator> ConstRange
0375             // here, MMap::const_iterator ConstIterator;
0376             auto range = node_cluster_map -> getClusters(hitsetkey);
0377 
0378             // loop over iterators of this cluster
0379             for (auto clusIter = range.first; clusIter != range.second; ++clusIter) {
0380                 const auto cluskey = clusIter -> first;
0381                 const auto cluster = clusIter -> second;
0382                 mvtx_clusters.push_back(cluster);
0383 
0384                 const auto globalPos = node_acts -> getGlobalPosition(cluskey, cluster);
0385 
0386                 // int ladder_z   = InttDefs::getLadderZId(cluskey);
0387                 // int ladder_phi = InttDefs::getLadderPhiId(cluskey);
0388                 int size = cluster -> getSize();
0389 
0390                 cluster -> setPosition(0, globalPos.x());
0391                 cluster -> setPosition(1, globalPos.y());
0392                 cluster -> setPosition(2, globalPos.z());
0393 
0394                 clus_system.push_back(0); // note : the cluster system is 0 for MVTX
0395                 clus_layer.push_back(mvtxlayer);
0396                 clus_adc.push_back(cluster -> getAdc());
0397                 clus_X.push_back(globalPos.x());
0398                 clus_Y.push_back(globalPos.y());
0399                 clus_Z.push_back(globalPos.z());
0400                 clus_size.push_back(size);
0401                 int phisize = cluster -> getPhiSize();
0402                 // if (phisize <= 0) {phisize += 256;}
0403                 clus_phi_size.push_back(phisize);
0404                 clus_z_size.push_back(cluster -> getZSize());
0405 
0406                 // cluster->setPosition(0,  globalPos.x() );
0407                 // cluster->setPosition(1,  globalPos.y() );
0408                 // cluster->setPosition(2,  globalPos.z() );
0409             }
0410         }
0411     }
0412 
0413     // note : intt
0414     for (unsigned int inttlayer = 0; inttlayer < 4; inttlayer++) {
0415         //      cout << " INTT layer " << inttlayer << endl;
0416         //      int layer= ( inttlayer < 2 ? 0 : 1 );
0417 
0418         // loop over all hits
0419         for (const auto & hitsetkey: node_cluster_map -> getHitSetKeys(TrkrDefs::TrkrId::inttId, inttlayer + 3)) {
0420 
0421             // type: std::pair<ConstIterator, ConstIterator> ConstRange
0422             // here, MMap::const_iterator ConstIterator;
0423             auto range = node_cluster_map -> getClusters(hitsetkey);
0424 
0425             // loop over iterators of this cluster
0426             for (auto clusIter = range.first; clusIter != range.second; ++clusIter) {
0427                 const auto cluskey = clusIter -> first;
0428                 const auto cluster = clusIter -> second;
0429                 intt_clusters.push_back(cluster);
0430 
0431                 const auto globalPos = node_acts -> getGlobalPosition(cluskey, cluster);
0432 
0433                 // int ladder_z   = InttDefs::getLadderZId(cluskey);
0434                 // int ladder_phi = InttDefs::getLadderPhiId(cluskey);
0435                 int size = cluster -> getSize();
0436 
0437                 cluster -> setPosition(0, globalPos.x());
0438                 cluster -> setPosition(1, globalPos.y());
0439                 cluster -> setPosition(2, globalPos.z());
0440 
0441                 clus_system.push_back(1); // note : the cluster system is 1 for INTT
0442                 clus_layer.push_back(inttlayer + 3);
0443                 clus_adc.push_back(cluster -> getAdc());
0444                 clus_X.push_back(globalPos.x());
0445                 clus_Y.push_back(globalPos.y());
0446                 clus_Z.push_back(globalPos.z());
0447                 clus_size.push_back(size);
0448                 int phisize = cluster -> getPhiSize();
0449                 if (phisize <= 0) {
0450                     phisize += 256;
0451                 }
0452                 clus_phi_size.push_back(phisize);
0453                 clus_z_size.push_back(cluster -> getZSize());
0454 
0455 
0456                 // cluster->setPosition(0,  globalPos.x() );
0457                 // cluster->setPosition(1,  globalPos.y() );
0458                 // cluster->setPosition(2,  globalPos.z() );
0459             }
0460         }
0461     }
0462 
0463     // note : TPC
0464     for (const auto & hitsetkey: node_cluster_map -> getHitSetKeys(TrkrDefs::TrkrId::tpcId)) {
0465 
0466         // type: std::pair<ConstIterator, ConstIterator> ConstRange
0467         // here, MMap::const_iterator ConstIterator;
0468         auto range = node_cluster_map -> getClusters(hitsetkey);
0469 
0470         // loop over iterators of this cluster
0471         for (auto clusIter = range.first; clusIter != range.second; ++clusIter) {
0472             const auto cluskey = clusIter -> first;
0473             const auto cluster = clusIter -> second;
0474             tpc_clusters.push_back(cluster);
0475 
0476             const auto globalPos = node_acts -> getGlobalPosition(cluskey, cluster);
0477 
0478             // int ladder_z   = InttDefs::getLadderZId(cluskey);
0479             // int ladder_phi = InttDefs::getLadderPhiId(cluskey);
0480             int size = cluster -> getSize();
0481 
0482             cluster -> setPosition(0, globalPos.x());
0483             cluster -> setPosition(1, globalPos.y());
0484             cluster -> setPosition(2, globalPos.z());
0485 
0486             clus_system.push_back(2); // note : the cluster system is 2 for TPC
0487             clus_layer.push_back(-999);
0488             clus_adc.push_back(cluster -> getAdc());
0489             clus_X.push_back(globalPos.x());
0490             clus_Y.push_back(globalPos.y());
0491             clus_Z.push_back(globalPos.z());
0492             clus_size.push_back(size);
0493             int phisize = cluster -> getPhiSize();
0494             // if (phisize <= 0) {phisize += 256;}
0495             clus_phi_size.push_back(phisize);
0496             clus_z_size.push_back(cluster -> getZSize());
0497 
0498             // cluster->setPosition(0,  globalPos.x() );
0499             // cluster->setPosition(1,  globalPos.y() );
0500             // cluster->setPosition(2,  globalPos.z() );
0501         }
0502     }
0503 
0504     NClus = mvtx_clusters.size() + intt_clusters.size() + tpc_clusters.size();
0505 
0506     std::cout << "N mvtx cluster : " << mvtx_clusters.size() << std::endl;
0507     std::cout << "N intt cluster : " << intt_clusters.size() << std::endl;
0508     std::cout << "N tpc cluster : " << tpc_clusters.size() << std::endl;
0509 
0510     return Fun4AllReturnCodes::EVENT_OK;
0511 
0512 }
0513 
0514 
0515 // CEMC (PHCompositeNode)/
0516 //        G4HIT_ABSORBER_CEMC (IO,PHG4HitContainer)
0517 //        G4HIT_CEMC (IO,PHG4HitContainer)
0518 //        G4CELL_CEMC (IO,PHG4CellContainer)
0519 //        TOWER_SIM_CEMC (IO,RawTowerContainer)
0520 //        TOWERINFO_SIM_CEMC (IO,TowerInfoContainerv1)
0521 //        TOWER_RAW_CEMC (IO,RawTowerContainer)
0522 //        TOWERINFO_RAW_CEMC (IO,TowerInfoContainerv1)
0523 //        TOWER_CALIB_CEMC (IO,RawTowerContainer)
0524 //        TOWERINFO_CALIB_CEMC (IO,TowerInfoContainerv1)
0525 //        CLUSTER_CEMC (IO,RawClusterContainer)
0526 //        CLUSTER_POS_COR_CEMC (IO,RawClusterContainer)
0527 
0528 int tutorial::prepareEMCal(PHCompositeNode * topNode) 
0529 {
0530     geomEM = findNode::getClass <RawTowerGeomContainer> (topNode, m_TowerGeomNodeName);
0531 
0532     EMCal_tower_sim_info = findNode::getClass <TowerInfoContainerv1> (topNode, "TOWERINFO_SIM_CEMC");
0533     EMCal_tower_sim = findNode::getClass <RawTowerContainer> (topNode, "TOWER_SIM_CEMC");
0534     EMCal_tower_calib = findNode::getClass <TowerInfoContainerv2> (topNode, emcal_node_name);
0535 
0536     if (!EMCal_tower_calib) 
0537     {
0538         std::cout << "tutorial::process_event Could not find node " << emcal_node_name << std::endl;
0539         exit(1);
0540     }
0541 
0542     // note : EMCal
0543     for (int i = 0; i < EMCal_tower_calib -> size(); ++i) //loop over channels 
0544     {
0545         TowerInfo * tower = EMCal_tower_calib -> get_tower_at_channel(i); //get EMCal tower
0546         int key = EMCal_tower_calib -> encode_key(i);
0547         int etabin = EMCal_tower_calib -> getTowerEtaBin(key);
0548         int phibin = EMCal_tower_calib -> getTowerPhiBin(key);
0549         // float time = EMCal_tower_calib->get_tower_at_channel(i)->get_time_float(); //get time
0550 
0551         if (tower -> get_energy() <= 0.07) {continue;}
0552 
0553         const RawTowerDefs::keytype geomkey = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, etabin, phibin);
0554         RawTowerGeom * tower_geom = geomEM -> get_tower_geometry(geomkey); //encode tower geometry
0555         double EMCal_pos_x   = tower_geom -> get_center_x();
0556         double EMCal_pos_y   = tower_geom -> get_center_y();
0557         double EMCal_pos_z   = tower_geom -> get_center_z();
0558         double EMCal_pos_R   = sqrt(EMCal_pos_x*EMCal_pos_x + EMCal_pos_y*EMCal_pos_y);
0559         double EMCal_pos_eta = tower_geom -> get_eta();
0560         double EMCal_pos_phi = tower_geom -> get_phi();
0561         double EMCal_energy  = tower -> get_energy();
0562 
0563         tower_system.push_back(0);
0564         tower_X.push_back(EMCal_pos_x);
0565         tower_Y.push_back(EMCal_pos_y);
0566         tower_Z.push_back(EMCal_pos_z);
0567         tower_R.push_back(EMCal_pos_R);
0568         tower_Eta.push_back(EMCal_pos_eta);
0569         tower_Phi.push_back(EMCal_pos_phi);
0570         tower_Eta_bin.push_back(etabin);
0571         tower_Phi_bin.push_back(phibin);
0572         tower_edep.push_back(EMCal_energy);
0573     
0574         double EMCal_int_x   = tower_geom -> get_center_int_x();
0575         double EMCal_int_y   = tower_geom -> get_center_int_y();
0576         double EMCal_int_z   = tower_geom -> get_center_int_z();
0577         double EMCal_int_R   = sqrt(EMCal_int_x*EMCal_int_x + EMCal_int_y*EMCal_int_y);
0578 
0579         tower_int_X.push_back(EMCal_int_x);
0580         tower_int_Y.push_back(EMCal_int_y);
0581         tower_int_Z.push_back(EMCal_int_z);
0582         tower_int_R.push_back(EMCal_int_R);
0583     }   
0584 
0585     return Fun4AllReturnCodes::EVENT_OK;
0586 }
0587 
0588 int tutorial::prepareiHCal(PHCompositeNode * topNode) {
0589     // HCALIN (PHCompositeNode)/
0590     //   G4HIT_ABSORBER_HCALIN (IO,PHG4HitContainer)
0591     //   G4HIT_HCALIN (IO,PHG4HitContainer)
0592     //   G4CELL_HCALIN (IO,PHG4CellContainer)
0593     //   TOWER_SIM_HCALIN (IO,RawTowerContainer)
0594     //   TOWERINFO_SIM_HCALIN (IO,TowerInfoContainerv1)
0595     //   TOWER_RAW_HCALIN (IO,RawTowerContainer)
0596     //   TOWERINFO_RAW_HCALIN (IO,TowerInfoContainerv1)
0597     //   TOWER_CALIB_HCALIN (IO,RawTowerContainer)
0598     //   TOWERINFO_CALIB_HCALIN (IO,TowerInfoContainerv1)
0599     //   CLUSTER_HCALIN (IO,RawClusterContainer)
0600 
0601     geomIH = findNode::getClass <RawTowerGeomContainer> (topNode, "TOWERGEOM_HCALIN");
0602     iHCal_tower_sim = findNode::getClass <RawTowerContainer> (topNode, "TOWER_SIM_HCALIN");
0603     iHCal_tower_calib = findNode::getClass <TowerInfoContainerv2> (topNode, ihcal_node_name);
0604 
0605     if (!iHCal_tower_calib) {
0606         std::cout << "tutorial::process_event Could not find node " << ihcal_node_name << std::endl;
0607         exit(1);
0608     }
0609 
0610     // note : iHCal
0611     for (int i = 0; i < iHCal_tower_calib -> size(); ++i) //loop over channels 
0612     {
0613         TowerInfo * tower = iHCal_tower_calib -> get_tower_at_channel(i); //get iHCal tower
0614         int key = iHCal_tower_calib -> encode_key(i);
0615         int etabin = iHCal_tower_calib -> getTowerEtaBin(key);
0616         int phibin = iHCal_tower_calib -> getTowerPhiBin(key);
0617         float time = iHCal_tower_calib->get_tower_at_channel(i)->get_time_float(); //get time
0618         if (tower -> get_energy() <= 0.07) {continue;}
0619 
0620         const RawTowerDefs::keytype geomkey = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, etabin, phibin);
0621         RawTowerGeom * tower_geom = geomIH -> get_tower_geometry(geomkey); //encode tower geometry
0622         double iHCal_pos_x   = tower_geom -> get_center_x();
0623         double iHCal_pos_y   = tower_geom -> get_center_y();
0624         double iHCal_pos_z   = tower_geom -> get_center_z();
0625         double iHCal_pos_eta = tower_geom -> get_eta();
0626         double iHCal_pos_phi = tower_geom -> get_phi();
0627         double iHCal_energy  = tower -> get_energy();
0628 
0629         tower_system.push_back(1);
0630         tower_X.push_back(iHCal_pos_x);
0631         tower_Y.push_back(iHCal_pos_y);
0632         tower_Z.push_back(iHCal_pos_z);
0633         tower_Eta.push_back(iHCal_pos_eta);
0634         tower_Phi.push_back(iHCal_pos_phi);
0635         tower_Eta_bin.push_back(etabin);
0636         tower_Phi_bin.push_back(phibin);
0637         tower_edep.push_back(iHCal_energy);
0638         // std::cout<<"test iHCal 15"<<std::endl;
0639     }
0640 
0641     return Fun4AllReturnCodes::EVENT_OK;
0642 }
0643 
0644 int tutorial::prepareoHCal(PHCompositeNode * topNode) {
0645     //     HCALOUT (PHCompositeNode)/
0646     //        G4HIT_ABSORBER_HCALOUT (IO,PHG4HitContainer)
0647     //        G4HIT_HCALOUT (IO,PHG4HitContainer)
0648     //        G4CELL_HCALOUT (IO,PHG4CellContainer)
0649     //        TOWER_SIM_HCALOUT (IO,RawTowerContainer)
0650     //        TOWERINFO_SIM_HCALOUT (IO,TowerInfoContainerv1)
0651     //        TOWER_RAW_HCALOUT (IO,RawTowerContainer)
0652     //        TOWERINFO_RAW_HCALOUT (IO,TowerInfoContainerv1)
0653     //        TOWER_CALIB_HCALOUT (IO,RawTowerContainer)
0654     //        TOWERINFO_CALIB_HCALOUT (IO,TowerInfoContainerv1)
0655     //        CLUSTER_HCALOUT (IO,RawClusterContainer)
0656 
0657     geomOH = findNode::getClass <RawTowerGeomContainer> (topNode, "TOWERGEOM_HCALOUT");
0658     oHCal_tower_sim = findNode::getClass <RawTowerContainer> (topNode, "TOWER_SIM_HCALOUT");
0659     oHCal_tower_calib = findNode::getClass <TowerInfoContainerv2> (topNode, ohcal_node_name);
0660 
0661     if (!oHCal_tower_calib) {
0662         std::cout << "tutorial::process_event Could not find node " << ohcal_node_name << std::endl;
0663         exit(1);
0664     }
0665 
0666     // note : oHCal
0667     for (int i = 0; i < oHCal_tower_calib -> size(); ++i) //loop over channels 
0668     {
0669         // std::cout<<"test oHCal 1"<<std::endl;
0670         TowerInfo * tower = oHCal_tower_calib -> get_tower_at_channel(i); //get oHCal tower
0671         // std::cout<<"test oHCal 2"<<std::endl;
0672         int key = oHCal_tower_calib -> encode_key(i);
0673         int etabin = oHCal_tower_calib -> getTowerEtaBin(key);
0674         int phibin = oHCal_tower_calib -> getTowerPhiBin(key);
0675 
0676         if (tower -> get_energy() <= 0.07) {continue;}
0677 
0678         const RawTowerDefs::keytype geomkey = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, etabin, phibin);
0679         RawTowerGeom * tower_geom = geomOH -> get_tower_geometry(geomkey); //encode tower geometry
0680         double oHCal_pos_x   = tower_geom -> get_center_x();
0681         double oHCal_pos_y   = tower_geom -> get_center_y();
0682         double oHCal_pos_z   = tower_geom -> get_center_z();
0683         double oHCal_pos_eta = tower_geom -> get_eta();
0684         double oHCal_pos_phi = tower_geom -> get_phi();
0685         double oHCal_energy  = tower -> get_energy();
0686 
0687         tower_system.push_back(2);
0688         tower_X.push_back(oHCal_pos_x);
0689         tower_Y.push_back(oHCal_pos_y);
0690         tower_Z.push_back(oHCal_pos_z);
0691         tower_Eta.push_back(oHCal_pos_eta);
0692         tower_Phi.push_back(oHCal_pos_phi);
0693         tower_Eta_bin.push_back(etabin);
0694         tower_Phi_bin.push_back(phibin);
0695         tower_edep.push_back(oHCal_energy);
0696     }
0697 
0698     return Fun4AllReturnCodes::EVENT_OK;
0699 }
0700 
0701 // === s === ChcKuma add Cal Clus ===============================================
0702 int tutorial::prepareEMCalClus(PHCompositeNode * topNode) {
0703     EMCal_cluster_cont = findNode::getClass <RawClusterContainer> (topNode, "CLUSTER_CEMC");
0704     EMCal_cluster_innr = findNode::getClass <RawClusterContainer> (topNode, "CLUSTERINNER_CEMC");
0705 
0706     if (!EMCal_cluster_cont) {
0707         std::cout << "tutorial::process_event Could not find node " << emcalClus_node_name << std::endl;
0708         exit(1);
0709     }
0710     if (!EMCal_cluster_innr) {
0711         std::cout << "tutorial::process_event Could not find node " << emcalClus_inner_node_name << std::endl;
0712         exit(1);
0713     }
0714 
0715     // note : EMCal clus center
0716     for (int i = 0; i < EMCal_cluster_cont -> size(); ++i) //loop over channels 
0717     {
0718         RawCluster * cluster = EMCal_cluster_cont  -> getCluster(i); //get EMCal tower
0719         if (cluster -> get_energy() <= 0.5) {continue;}
0720 
0721         double EMCal_pos_x   = cluster -> get_x();
0722         double EMCal_pos_y   = cluster -> get_y();
0723         double EMCal_pos_z   = cluster -> get_z();
0724         double EMCal_pos_r   = cluster -> get_r();
0725         double EMCal_pos_phi = cluster -> get_phi();
0726         double EMCal_energy  = cluster -> get_energy();
0727 
0728         caloClus_system.push_back(0);
0729         caloClus_X.push_back(EMCal_pos_x);
0730         caloClus_Y.push_back(EMCal_pos_y);
0731         caloClus_Z.push_back(EMCal_pos_z);
0732         caloClus_R.push_back(EMCal_pos_r);
0733         caloClus_Phi.push_back(EMCal_pos_phi);
0734         caloClus_edep.push_back(EMCal_energy);
0735     }   
0736 
0737     // note : EMCal clus innerface center
0738     for (int i = 0; i < EMCal_cluster_innr -> size(); ++i) //loop over channels 
0739     {
0740         RawCluster * cluster_innr = EMCal_cluster_innr  -> getCluster(i); //get EMCal tower
0741         if (cluster_innr -> get_energy() <= 0.5) {continue;}
0742 
0743         double EMCal_innr_x   = cluster_innr -> get_x();
0744         double EMCal_innr_y   = cluster_innr -> get_y();
0745         double EMCal_innr_z   = cluster_innr -> get_z();
0746         double EMCal_innr_r   = cluster_innr -> get_r();
0747         double EMCal_innr_phi = cluster_innr -> get_phi();
0748         double EMCal_innr_energy  = cluster_innr -> get_energy();
0749 
0750         caloClus_system.push_back(321);
0751         caloClus_innr_X.push_back(EMCal_innr_x);
0752         caloClus_innr_Y.push_back(EMCal_innr_y);
0753         caloClus_innr_Z.push_back(EMCal_innr_z);
0754         caloClus_innr_R.push_back(EMCal_innr_r);
0755         caloClus_innr_Phi.push_back(EMCal_innr_phi);
0756         caloClus_innr_edep.push_back(EMCal_innr_energy);
0757     }
0758 
0759     return Fun4AllReturnCodes::EVENT_OK;
0760 }
0761 
0762 
0763 int tutorial::prepareiHCalClus(PHCompositeNode * topNode) {
0764     iHCal_cluster_cont = findNode::getClass <RawClusterContainer> (topNode, "CLUSTER_HCALIN");
0765 
0766     if (!iHCal_cluster_cont) {
0767         std::cout << "tutorial::process_event Could not find node " << ihcalClus_node_name << std::endl;
0768         exit(1);
0769     }
0770 
0771     // note : iHCal
0772     for (int i = 0; i < iHCal_cluster_cont -> size(); ++i) //loop over channels 
0773     {
0774         RawCluster * cluster = iHCal_cluster_cont  -> getCluster(i); //get EMCal tower
0775         if (cluster -> get_energy() <= 0.5) {continue;}
0776 
0777         double iHCal_pos_x   = cluster -> get_x();
0778         double iHCal_pos_y   = cluster -> get_y();
0779         double iHCal_pos_z   = cluster -> get_z();
0780         double iHCal_pos_r   = cluster -> get_r();
0781         double iHCal_pos_phi = cluster -> get_phi();
0782         double iHCal_energy  = cluster -> get_energy();
0783 
0784         caloClus_system.push_back(1);
0785         caloClus_X.push_back(iHCal_pos_x);
0786         caloClus_Y.push_back(iHCal_pos_y);
0787         caloClus_Z.push_back(iHCal_pos_z);
0788         caloClus_R.push_back(iHCal_pos_r);
0789         caloClus_Phi.push_back(iHCal_pos_phi);
0790         caloClus_edep.push_back(iHCal_energy);
0791     }   
0792     return Fun4AllReturnCodes::EVENT_OK;
0793 }
0794 
0795 
0796 int tutorial::prepareoHCalClus(PHCompositeNode * topNode) {
0797     oHCal_cluster_cont = findNode::getClass <RawClusterContainer> (topNode, "CLUSTER_HCALIN");
0798 
0799     if (!oHCal_cluster_cont) {
0800         std::cout << "tutorial::process_event Could not find node " << ohcalClus_node_name << std::endl;
0801         exit(1);
0802     }
0803 
0804     // note : oHCal
0805     for (int i = 0; i < oHCal_cluster_cont -> size(); ++i) //loop over channels 
0806     {
0807         RawCluster * cluster = oHCal_cluster_cont  -> getCluster(i); //get EMCal tower
0808         if (cluster -> get_energy() <= 0.5) {continue;}
0809 
0810         double oHCal_pos_x   = cluster -> get_x();
0811         double oHCal_pos_y   = cluster -> get_y();
0812         double oHCal_pos_z   = cluster -> get_z();
0813         double oHCal_pos_r   = cluster -> get_r();
0814         double oHCal_pos_phi = cluster -> get_phi();
0815         double oHCal_energy  = cluster -> get_energy();
0816 
0817         caloClus_system.push_back(2);
0818         caloClus_X.push_back(oHCal_pos_x);
0819         caloClus_Y.push_back(oHCal_pos_y);
0820         caloClus_Z.push_back(oHCal_pos_z);
0821         caloClus_R.push_back(oHCal_pos_r);
0822         caloClus_Phi.push_back(oHCal_pos_phi);
0823         caloClus_edep.push_back(oHCal_energy);
0824     }   
0825     return Fun4AllReturnCodes::EVENT_OK;
0826 }
0827 
0828 // === e === ChcKuma add Cal Clus ===============================================
0829 
0830 int tutorial::prepareG4Turth(PHCompositeNode * topNode){
0831     std::cout << "Get PHG4 info.: truth primary vertex" << std::endl;
0832     m_truth_info = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");   
0833 
0834     if (!m_truth_info)
0835     {
0836         std::cout << PHWHERE << "Error, can't find G4TruthInfo" << std::endl;
0837         exit(1);
0838     }
0839 
0840     // note : Truth vertex
0841     auto vrange = m_truth_info->GetPrimaryVtxRange();
0842     int NTruthPV = 0, NTruthPV_Embeded0 = 0;
0843     for (auto iter = vrange.first; iter != vrange.second; ++iter) // process all primary vertices
0844     {
0845         const int point_id = iter->first;
0846         PHG4VtxPoint *point = iter->second;
0847         if (point)
0848         {
0849             if (m_truth_info->isEmbededVtx(point_id) == 0)
0850             {
0851                 TruthPV_trig_x_ = point->get_x();
0852                 TruthPV_trig_y_ = point->get_y();
0853                 TruthPV_trig_z_ = point->get_z();
0854                 
0855                 NTruthPV_Embeded0++;
0856             }
0857             NTruthPV++;
0858         }
0859     }
0860 
0861     NTruthVtx_ = NTruthPV;
0862 
0863     // note : PHG4Particle
0864     std::vector<int> tmpv_chargehadron;
0865     tmpv_chargehadron.clear();
0866     std::cout << "Get PHG4 info.: truth primary G4Particle" << std::endl;
0867     const auto prange = m_truth_info->GetPrimaryParticleRange();
0868     // const auto prange = m_truth_info->GetParticleRange();
0869     for (auto iter = prange.first; iter != prange.second; ++iter)
0870     {
0871         PHG4Particle *ptcl = iter->second;
0872         // particle->identify();
0873         if (ptcl)
0874         {
0875             PrimaryG4P_PID_.push_back(ptcl->get_pid());
0876             TLorentzVector p;
0877             p.SetPxPyPzE(ptcl->get_px(), ptcl->get_py(), ptcl->get_pz(), ptcl->get_e());
0878             PrimaryG4P_E_.push_back(ptcl->get_e());
0879             PrimaryG4P_Pt_.push_back(p.Pt());
0880             PrimaryG4P_Eta_.push_back(p.Eta());
0881             PrimaryG4P_Phi_.push_back(p.Phi());
0882 
0883             TString particleclass = TString(TDatabasePDG::Instance()->GetParticle(ptcl->get_pid())->ParticleClass());
0884             bool isStable = (TDatabasePDG::Instance()->GetParticle(ptcl->get_pid())->Stable() == 1) ? true : false;
0885             double charge = TDatabasePDG::Instance()->GetParticle(ptcl->get_pid())->Charge();
0886             bool isHadron = (particleclass.Contains("Baryon") || particleclass.Contains("Meson"));
0887             bool isChargeHadron = (isStable && (charge != 0) && isHadron);
0888             if (isChargeHadron)
0889                 tmpv_chargehadron.push_back(ptcl->get_pid());
0890 
0891             PrimaryG4P_ParticleClass_.push_back(particleclass);
0892             PrimaryG4P_isStable_.push_back(isStable);
0893             PrimaryG4P_Charge_.push_back(charge);
0894             PrimaryG4P_isChargeHadron_.push_back(isChargeHadron);
0895 
0896             // 找到pr electron track id
0897             if (abs(ptcl->get_pid()) == 11)  // 电子的 PDG ID 是 ±11
0898             {
0899                 primary_electron_tracks.insert(ptcl->get_track_id());
0900             }
0901         }
0902     }
0903     NPrimaryG4P_ = PrimaryG4P_PID_.size();
0904     NPrimaryG4P_promptChargeHadron_ = tmpv_chargehadron.size();
0905 
0906     return Fun4AllReturnCodes::EVENT_OK;
0907 }
0908 
0909 // -------------------------------------------------------------------------------------------------------------------------- 
0910 int tutorial::prepareG4HIT(PHCompositeNode * topNode)
0911 {
0912     if (primary_electron_tracks.empty())
0913     {
0914         std::cout << "No primary electrons found in this event!" << std::endl;
0915         return 0;
0916     }
0917 
0918     hits_CEMC = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_CEMC");
0919     if(!hits_CEMC)
0920     {
0921       std::cout << "PhotonEMC::process_event: G4HIT_CEMC not found!!!" << std::endl;
0922     }
0923 
0924     int prhit0 = 0;
0925     double t_earlist = 0.;
0926     PHG4HitContainer::ConstRange hit_range = hits_CEMC->getHits();
0927     for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0928     {
0929         PHG4Hit *this_hit = hit_iter->second;
0930         float light_yield = hit_iter->second->get_light_yield();
0931         float edep = hit_iter->second->get_edep();
0932         int ch = hit_iter->second->get_layer();
0933         float x = hit_iter->second->get_x(0);
0934         float y = hit_iter->second->get_y(0);
0935         float z = hit_iter->second->get_z(0);
0936         
0937         // int trkid = hit_iter->second->get_trkid();
0938         // PHG4Particle *part = truthinfo->GetParticle(trkid);
0939         // v4.SetPxPyPzE(part->get_px(), part->get_py(), part->get_pz(), part->get_e());
0940         // int pid = part->get_pid();
0941         // _CEMC_Hit_pid.push_back(pid);
0942 
0943         // int vtxid = part->get_vtx_id();
0944         // PHG4VtxPoint *vtx = truthinfo->GetVtx(vtxid);
0945 
0946         // // add trkid to a set
0947         // _CEMC_Hit_Evis.push_back(light_yield);
0948         // _CEMC_Hit_Edep.push_back(edep);
0949         // _CEMC_Hit_ch.push_back(ch);
0950         // _CEMC_Hit_x.push_back(x);
0951         // _CEMC_Hit_y.push_back(y);
0952         // _CEMC_Hit_z.push_back(z);
0953 
0954         // _CEMC_Hit_particle_x.push_back(vtx->get_x());
0955         // _CEMC_Hit_particle_y.push_back(vtx->get_y());
0956         // _CEMC_Hit_particle_z.push_back(vtx->get_z());
0957 
0958         // get primary electron hits on emcal
0959         int track_id = this_hit->get_trkid();
0960         // if (primary_electron_tracks.find(track_id) == primary_electron_tracks.end()) continue;
0961 
0962         // 只选择初级电子的 G4Hit，忽略由 shower 产生的次级粒子
0963         if (primary_electron_tracks.find(track_id) != primary_electron_tracks.end())
0964         {
0965             prhit0 += 1;
0966 
0967             double x0 = this_hit->get_x(0);  // 入口 x
0968             double y0 = this_hit->get_y(0);  // 入口 y
0969             double z0 = this_hit->get_z(0);  // 入口 z
0970             double t0 = this_hit->get_t(0);  // 时间
0971 
0972             double r0 = sqrt(x0*x0 + y0*y0);
0973 
0974             // std::cout << "Primary electron hit0 EMCal at: (x0=" << x0
0975             //           << ", y0=" << y0 << ", z0=" << z0 << ", t0=" << t0 << ", r0="<< r0 << " )" << std::endl;
0976 
0977             double x1 = this_hit->get_x(1);  // 出口 x
0978             double y1 = this_hit->get_y(1);  // 出口 y
0979             double z1 = this_hit->get_z(1);  // 出口 z
0980             double t1 = this_hit->get_t(1);  // 时间
0981           
0982             // std::cout << "Primary electron hit1 EMCal at: (x1=" << x1
0983             //           << ", y1=" << y1 << ", z1=" << z1 << ", t1=" << t1 << ")" << std::endl;
0984 
0985             double delta_t = t0 - t_earlist;
0986 
0987             if (prhit0 == 1) 
0988             {
0989                 t_earlist = t0;
0990 
0991                 _CEMC_Pr_Hit_x.push_back(x0);
0992                 _CEMC_Pr_Hit_y.push_back(y0);
0993                 _CEMC_Pr_Hit_z.push_back(z0);
0994                 _CEMC_Pr_Hit_R.push_back(r0);
0995 
0996                 double PtG_x = x0; 
0997                 double PtG_y = y0; 
0998                 double PtG_r = r0;  
0999                 double PtG_phi = atan2(PtG_y, PtG_x); 
1000                 double PtG_z = z0; 
1001 
1002                 // std::cout << "PrG at EMCal surface: (x=" << PtG_x
1003                 //           << ", y=" << PtG_y << ", r=" << PtG_r << ", phi=" << PtG_phi
1004                 //           << ", z=" << PtG_z << ")" << std::endl;
1005             }
1006             _CEMC_Pr_Hit_deltaT.push_back(delta_t);
1007 
1008         }
1009     }
1010 
1011     return Fun4AllReturnCodes::EVENT_OK;
1012 }
1013 
1014 int tutorial::createTracksFromTruth(PHCompositeNode* topNode)
1015 {
1016     // 获取或创建 SvtxTrackMap 节点
1017     SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
1018     if (!trackmap)
1019     {
1020         trackmap = new SvtxTrackMap_v1();
1021         PHNodeIterator iter(topNode);
1022         PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
1023 
1024         PHIODataNode<PHObject> *truthtracknode = new PHIODataNode<PHObject>(trackmap, "SvtxTrackMap", "PHObject");
1025         dstNode->addNode(truthtracknode);
1026     }
1027 
1028     const auto prange = m_truth_info->GetPrimaryParticleRange();
1029     for (auto iter = prange.first; iter != prange.second; ++iter)
1030     {
1031         PHG4Particle* ptcl = iter->second;
1032         if (!ptcl) continue;
1033 
1034         TLorentzVector p;
1035         p.SetPxPyPzE(ptcl->get_px(), ptcl->get_py(), ptcl->get_pz(), ptcl->get_e());
1036 
1037         SvtxTrack_v2* track = new SvtxTrack_v2();
1038         track->set_id(trackmap->size());
1039         track->set_px(ptcl->get_px());
1040         track->set_py(ptcl->get_py());
1041         track->set_pz(ptcl->get_pz());
1042         track->set_charge(static_cast<int>(TDatabasePDG::Instance()->GetParticle(ptcl->get_pid())->Charge()));
1043         track->set_chisq(1.0);
1044         track->set_ndf(1);
1045         track->set_vertex_id(0);  // 可改为对应 vtx id
1046 
1047         trackmap->insert(track);
1048     }
1049 
1050     return Fun4AllReturnCodes::EVENT_OK;
1051 }
1052 
1053 int tutorial::prepareTruthTrack(PHCompositeNode * topNode) 
1054 {
1055     trackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
1056     if(!trackMap)
1057     {
1058       std::cout << "TrackCaloMatch::process_event not found! Aborting!" << std::endl;
1059       return Fun4AllReturnCodes::ABORTEVENT;
1060     }
1061 
1062     for (auto &iter : *trackMap)
1063     {
1064         truth_track = iter.second;
1065         
1066         double truth_tk_pt = truth_track->get_pt();
1067         std::cout<< "Truth track pt: " << truth_tk_pt << std::endl;
1068 
1069         thisState = truth_track->get_state(99);
1070         double em_x = thisState->get_x();
1071         double em_y = thisState->get_y();
1072         double em_r = sqrt(em_x*em_x + em_y*em_y);
1073         double em_phi = thisState->get_phi();
1074         double em_z = thisState->get_z();
1075         // std::cout << "Truth track state at EMCal: (x=" << em_x
1076         //           << ", y=" << em_y << ", r=" << em_r << ", phi=" << em_phi
1077         //           << ", z=" << em_z << ")" << std::endl;
1078 
1079         _PrG4_TTPRO_dD.push_back(sqrt((em_x-_CEMC_Pr_Hit_x[0])*(em_x-_CEMC_Pr_Hit_x[0])+(em_y-_CEMC_Pr_Hit_y[0])*(em_y-_CEMC_Pr_Hit_y[0]))); // 计算与初级电子的距离差
1080         _PrG4_TTPRO_dR.push_back(em_r-_CEMC_Pr_Hit_R[0]);
1081         _PrG4_TTPRO_dphi.push_back(em_phi-atan2(_CEMC_Pr_Hit_y[0], _CEMC_Pr_Hit_x[0])); // 计算与初级电子的角度差
1082     }
1083 
1084     return Fun4AllReturnCodes::EVENT_OK;
1085 }
1086 
1087 
1088 //____________________________________________________________________________..
1089 int tutorial::process_event(PHCompositeNode * topNode) {
1090 
1091     prepareTracker(topNode);
1092     prepareEMCal(topNode);
1093     prepareiHCal(topNode);
1094     prepareoHCal(topNode);
1095     prepareEMCalClus(topNode);
1096     prepareiHCalClus(topNode);
1097     prepareoHCalClus(topNode);
1098     
1099     prepareG4Turth(topNode);
1100     prepareG4HIT(topNode);
1101 
1102     prepareTruthTrack(topNode);
1103 
1104     // createTracksFromTruth(topNode);
1105 
1106     nTowers = tower_system.size();
1107     
1108     tree_out -> Fill();
1109 
1110     return Fun4AllReturnCodes::EVENT_OK;
1111 }
1112 
1113 
1114 //____________________________________________________________________________..
1115 int tutorial::ResetEvent(PHCompositeNode *topNode)
1116 {
1117     std::cout << "tutorial::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
1118 
1119     NClus = 0;
1120     clus_system.clear();
1121     clus_layer.clear();
1122     clus_adc.clear();
1123     clus_X.clear();
1124     clus_Y.clear();
1125     clus_Z.clear();
1126     clus_size.clear();
1127     clus_phi_size.clear();
1128     clus_z_size.clear();
1129 
1130     nTowers = 0;
1131     tower_system.clear();
1132     tower_X.clear();
1133     tower_Y.clear();
1134     tower_Z.clear();
1135     tower_R.clear();
1136     tower_Eta.clear();
1137     tower_Phi.clear();
1138     tower_Eta_test.clear();
1139     tower_Phi_test.clear();
1140     tower_Eta_bin.clear();
1141     tower_Phi_bin.clear();
1142     tower_edep.clear();
1143 
1144     tower_int_X.clear();
1145     tower_int_Y.clear();
1146     tower_int_Z.clear();
1147     tower_int_R.clear();
1148 
1149     nCaloClus = 0;
1150     caloClus_system.clear();
1151     caloClus_X.clear();
1152     caloClus_Y.clear();
1153     caloClus_Z.clear();
1154     caloClus_R.clear();
1155     caloClus_Phi.clear();
1156     caloClus_edep.clear();
1157 
1158     caloClus_innr_X.clear();
1159     caloClus_innr_Y.clear();
1160     caloClus_innr_Z.clear();
1161     caloClus_innr_R.clear();
1162     caloClus_innr_Phi.clear();
1163     caloClus_innr_edep.clear();
1164 
1165     // note : Truth primary vertex information
1166     TruthPV_trig_x_ = -999;
1167     TruthPV_trig_y_ = -999;
1168     TruthPV_trig_z_ = -999;
1169     NTruthVtx_ = 0;
1170     // note : PHG4 information (from all PHG4Particles)
1171     NPrimaryG4P_ = 0;
1172     NPrimaryG4P_promptChargeHadron_ = 0;
1173     PrimaryG4P_Pt_.clear();
1174     PrimaryG4P_Eta_.clear();
1175     PrimaryG4P_Phi_.clear();
1176     PrimaryG4P_E_.clear();
1177     PrimaryG4P_PID_.clear();
1178     PrimaryG4P_ParticleClass_.clear();
1179     PrimaryG4P_isStable_.clear();
1180     PrimaryG4P_Charge_.clear();
1181     PrimaryG4P_isChargeHadron_.clear();
1182 
1183     _CEMC_Hit_Evis.clear();
1184     _CEMC_Hit_Edep.clear();
1185     _CEMC_Hit_ch.clear();
1186     _CEMC_Hit_x.clear();
1187     _CEMC_Hit_y.clear();
1188     _CEMC_Hit_z.clear();
1189 
1190     _CEMC_Pr_Hit_x.clear();
1191     _CEMC_Pr_Hit_y.clear();
1192     _CEMC_Pr_Hit_z.clear();
1193     _CEMC_Pr_Hit_R.clear();
1194     _CEMC_Pr_Hit_deltaT.clear();
1195 
1196     primary_electron_tracks.clear();
1197 
1198     eventID += 1;
1199 
1200 
1201     // 2. 再对所有 vector 进行 shrink_to_fit 操作，释放多余内存
1202     clus_system.shrink_to_fit();
1203     clus_layer.shrink_to_fit();
1204     clus_adc.shrink_to_fit();
1205     clus_X.shrink_to_fit();
1206     clus_Y.shrink_to_fit();
1207     clus_Z.shrink_to_fit();
1208     clus_size.shrink_to_fit();
1209     clus_phi_size.shrink_to_fit();
1210     clus_z_size.shrink_to_fit();
1211 
1212     tower_system.shrink_to_fit();
1213     tower_X.shrink_to_fit();
1214     tower_Y.shrink_to_fit();
1215     tower_Z.shrink_to_fit();
1216     tower_R.shrink_to_fit();
1217     tower_Eta.shrink_to_fit();
1218     tower_Phi.shrink_to_fit();
1219     tower_Eta_test.shrink_to_fit();
1220     tower_Phi_test.shrink_to_fit();
1221     tower_Eta_bin.shrink_to_fit();
1222     tower_Phi_bin.shrink_to_fit();
1223     tower_edep.shrink_to_fit();
1224 
1225     tower_int_X.shrink_to_fit();
1226     tower_int_Y.shrink_to_fit();
1227     tower_int_Z.shrink_to_fit();
1228     tower_int_R.shrink_to_fit();
1229 
1230     caloClus_system.shrink_to_fit();
1231     caloClus_X.shrink_to_fit();
1232     caloClus_Y.shrink_to_fit();
1233     caloClus_Z.shrink_to_fit();
1234     caloClus_R.shrink_to_fit();
1235     caloClus_Phi.shrink_to_fit();
1236     caloClus_edep.shrink_to_fit();
1237 
1238     caloClus_innr_X.shrink_to_fit();
1239     caloClus_innr_Y.shrink_to_fit();
1240     caloClus_innr_Z.shrink_to_fit();
1241     caloClus_innr_R.shrink_to_fit();
1242     caloClus_innr_Phi.shrink_to_fit();
1243     caloClus_innr_edep.shrink_to_fit();
1244 
1245     PrimaryG4P_Pt_.shrink_to_fit();
1246     PrimaryG4P_Eta_.shrink_to_fit();
1247     PrimaryG4P_Phi_.shrink_to_fit();
1248     PrimaryG4P_E_.shrink_to_fit();
1249     PrimaryG4P_PID_.shrink_to_fit();
1250     PrimaryG4P_ParticleClass_.shrink_to_fit();
1251     PrimaryG4P_isStable_.shrink_to_fit();
1252     PrimaryG4P_Charge_.shrink_to_fit();
1253     PrimaryG4P_isChargeHadron_.shrink_to_fit();
1254 
1255     _CEMC_Hit_Evis.shrink_to_fit();
1256     _CEMC_Hit_Edep.shrink_to_fit();
1257     _CEMC_Hit_ch.shrink_to_fit();
1258     _CEMC_Hit_x.shrink_to_fit();
1259     _CEMC_Hit_y.shrink_to_fit();
1260     _CEMC_Hit_z.shrink_to_fit();
1261 
1262     _CEMC_Pr_Hit_x.shrink_to_fit();
1263     _CEMC_Pr_Hit_y.shrink_to_fit();
1264     _CEMC_Pr_Hit_z.shrink_to_fit();
1265     _CEMC_Pr_Hit_R.shrink_to_fit();
1266     _CEMC_Pr_Hit_deltaT.shrink_to_fit();
1267 
1268     return Fun4AllReturnCodes::EVENT_OK;
1269 }
1270 
1271 //____________________________________________________________________________..
1272 int tutorial::EndRun(const int runnumber)
1273 {
1274     std::cout << "tutorial::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
1275 
1276     return Fun4AllReturnCodes::EVENT_OK;
1277 }
1278 
1279 //____________________________________________________________________________..
1280 int tutorial::End(PHCompositeNode *topNode)
1281 {
1282   std::cout << "tutorial::End(PHCompositeNode *topNode) This is the End..." << std::endl;
1283 
1284   ////////////////////////////////////////////////////////
1285   // Writing objects to the output file                 //
1286   ////////////////////////////////////////////////////////
1287   file_out -> cd();
1288   tree_out -> Write();
1289   file_out -> Close();
1290 
1291   return Fun4AllReturnCodes::EVENT_OK;
1292 }
1293 
1294 //____________________________________________________________________________..
1295 int tutorial::Reset(PHCompositeNode *topNode)
1296 {
1297  std::cout << "tutorial::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
1298   return Fun4AllReturnCodes::EVENT_OK;
1299 }
1300 
1301 //____________________________________________________________________________..
1302 void tutorial::Print(const std::string &what) const
1303 {
1304   std::cout << "tutorial::Print(const std::string &what) const Printing info for " << what << std::endl;
1305 }
1306 
1307 
1308 
1309 //____________________________________________________________________________..
1310 bool tutorial::checkTrack(SvtxTrack* track)
1311 {
1312     return true;
1313 }
[ 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!