FullJetFinder/src/FullJetFinder.cc

0001 /**
0002  * @file FullJetFinder.cc
0003  *
0004  * @brief module for producing a TTree with full jets (tracks + calos) studies
0005  *        Based on JetValidation macro and HF group macros
0006  *
0007  * @ingroup FullJetFinder
0008  *
0009  * @author Jakub Kvapil
0010  * Contact: jakub.kvapil@cern.ch
0011  *
0012  */
0013
0014 #include "FullJetFinder.h"
0015
0016 #include <fun4all/Fun4AllReturnCodes.h>
0017 #include <fun4all/PHTFileServer.h>
0018
0019 #include <phool/PHCompositeNode.h>
0020 #include <phool/getClass.h>
0021
0022 #include <jetbase/JetContainer.h>
0023 #include <jetbase/JetMap.h>
0024 #include <jetbase/Jet.h>
0025
0026 #include <centrality/CentralityInfo.h>
0027
0028 #include <globalvertex/GlobalVertex.h>
0029 #include <globalvertex/GlobalVertexMap.h>
0030
0031 #include <particleflowreco/ParticleFlowElementv1.h>
0032 #include <particleflowreco/ParticleFlowElementContainer.h>
0033
0034 #define HomogeneousField
0035 #include <KFParticle.h>
0036
0037 #include <g4main/PHG4TruthInfoContainer.h>
0038 #include <g4main/PHG4Particle.h>
0039
0040 #include <trackbase_historic/SvtxTrackMap.h>
0041 #include <trackbase_historic/TrackAnalysisUtils.h>
0042
0043 #include <trackbase/MvtxDefs.h>
0044
0045 /// HEPMC truth includes
0046 #pragma GCC diagnostic push
0047 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0048 #include <HepMC/GenEvent.h>
0049 #include <HepMC/GenVertex.h>
0050 #pragma GCC diagnostic pop
0051
0052 #include <phhepmc/PHHepMCGenEvent.h>
0053 #include <phhepmc/PHHepMCGenEventMap.h>
0054
0055 #include <TDatabasePDG.h>
0056
0057 #include <TTree.h>
0058 #include <TH1.h>
0059
0060 #include <cmath>
0061
0062
0063
0064 //____________________________________________________________________________..
0065 FullJetFinder::FullJetFinder(const std::string& outputfilename, FullJetFinder::TYPE jet_type):
0066  SubsysReco("FullJetFinder")
0067  , m_recoJetName()
0068  , m_truthJetName()
0069  , m_outputFileName(outputfilename)
0070  , m_ptRangeReco(0, 100)
0071  , m_ptRangeTruth(0, 100)
0072  , m_doTruthJets(0)
0073  , m_T()
0074  , m_jet_type(jet_type)
0075 {
0076   std::cout << "FullJetFinder::FullJetFinder(const std::string &name) Calling ctor" << std::endl;
0077 }
0078
0079 //____________________________________________________________________________..
0080 FullJetFinder::~FullJetFinder()
0081 {
0082   std::cout << "FullJetFinder::~FullJetFinder() Calling dtor" << std::endl;
0083 }
0084
0085 //_________________________________________________________
0086 void FullJetFinder::Container::Reset()
0087 {
0088   reco_jet_n = -1;
0089   truth_jet_n = -1;
0090   centrality = -1;
0091   impactparam = -1;
0092   recojets.clear();
0093   truthjets.clear();
0094   //primaryVertex.clear();
0095 }
0096
0097 //____________________________________________________________________________..
0098 int FullJetFinder::Init(PHCompositeNode *topNode)
0099 {
0100   std::cout << "FullJetFinder::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0101   PHTFileServer::get().open(m_outputFileName, "RECREATE");
0102   std::cout << "FullJetFinder::Init - Output to " << m_outputFileName << std::endl;
0103
0104   if (m_inputs == 0 || m_inputs > 10){
0105     std::cout << "MyJetAnalysis::Init - Error number of inputs outside (0,10> " << std::endl;
0106     exit(-1);
0107   }
0108
0109   // configure Tree
0110   for(int i = 0 ; i < m_inputs; i++){
0111     //m_T[i] = new TTree(m_TreeNameCollection.at(i).data(), m_TreeNameCollection.at(i).data());
0112     m_T[i] = new TTree("Data", "Data");
0113     m_container[i] = new Container;
0114     m_T[i]->Branch( "JetContainer", &m_container[i] );
0115
0116     TString tmp = "";
0117     m_chi2ndf[i] = new TH1D(tmp +"chi2ndf_" + m_TreeNameCollection.at(i).data(),tmp +"chi2ndf_" + m_TreeNameCollection.at(i).data(),1000,0,100);
0118     m_mvtxcl[i] = new TH1I(tmp +"mvtxcl_" + m_TreeNameCollection.at(i).data(),tmp +"mvtxcl_" + m_TreeNameCollection.at(i).data(),6,0,6);
0119     m_inttcl[i] = new TH1I(tmp +"inttcl_" + m_TreeNameCollection.at(i).data(),tmp +"inttcl_" + m_TreeNameCollection.at(i).data(),6,0,6);
0120     m_mtpccl[i]= new TH1I(tmp +"tpccl_" + m_TreeNameCollection.at(i).data(),tmp +"tpccl_" + m_TreeNameCollection.at(i).data(),100,0,100);
0121   }
0122
0123   m_stat = new TH1I("event_stat","event_stat",10,-0.5,9.5);
0124   m_stat->GetXaxis()->SetBinLabel(1,"n_events");
0125   m_stat->GetXaxis()->SetBinLabel(2,"GV_exists");
0126   m_stat->GetXaxis()->SetBinLabel(3,"GV_notempty");
0127   m_stat->GetXaxis()->SetBinLabel(4,"GV_SVTX_vtx");
0128   m_stat->GetXaxis()->SetBinLabel(5,"GV_in_10cm");
0129
0130
0131
0132   return Fun4AllReturnCodes::EVENT_OK;
0133 }
0134
0135 //____________________________________________________________________________..
0136 int FullJetFinder::InitRun(PHCompositeNode *topNode)
0137 {
0138   std::cout << "FullJetFinder::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0139   return Fun4AllReturnCodes::EVENT_OK;
0140 }
0141
0142 //____________________________________________________________________________..
0143 int FullJetFinder::process_event(PHCompositeNode *topNode)
0144 {
0145     //std::cout<<"NEW EVENT"<<std::endl;
0146   //std::cout<<"NEW Fun4All EVENT"<<std::endl<<std::endl;
0147   m_stat->Fill(0);
0148   //centrality
0149   CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0150   if (!cent_node){
0151     std::cout << "MyJetAnalysis::process_event - Error can not find centrality node " << std::endl;
0152     exit(-1);
0153   }
0154
0155   // interface to reco jets
0156   if ((m_recoJetName.size() != m_truthJetName.size()) || m_recoJetName.empty() || m_truthJetName.empty()){
0157     std::cout << "MyJetAnalysis::process_event - Error in m_recoJetName size !=  m_truthJetName.size() or empty" << std::endl;
0158     exit(-1);
0159   }
0160
0161   if (m_recoJetName.size() != static_cast<long unsigned int>(m_inputs)){
0162     std::cout << "MyJetAnalysis::process_event - Error number of AddInput() calls does not match with number of inputs specified in constructor" << std::endl;
0163     exit(-1);
0164   }
0165
0166
0167
0168   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0169   if (!vertexmap){
0170     std::cout << "MyJetAnalysis::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;
0171     assert(vertexmap);  // force quit
0172     exit(-1);
0173   }
0174
0175   m_stat->Fill(1);
0176
0177   if (vertexmap->empty()){
0178     std::cout << "MyJetAnalysis::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;
0179     return Fun4AllReturnCodes::ABORTEVENT;
0180   }
0181
0182   m_stat->Fill(2);
0183
0184   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0185   if (!truthinfo)
0186   {
0187     std::cout << PHWHERE << " ERROR: Can't find G4TruthInfo" << std::endl;
0188      exit(-1);
0189   }
0190
0191   //loop over multiple JetReco->add_algo()
0192   for(long unsigned int input = 0; input < m_recoJetName.size(); input++){
0193     //reset per event container
0194     m_container[input]->Reset();
0195
0196     //get the event centrality/impact parameter from HIJING
0197     m_container[input]->centrality =  cent_node->get_centile(CentralityInfo::PROP::bimp);
0198     m_container[input]->impactparam =  cent_node->get_quantity(CentralityInfo::PROP::bimp);
0199
0200
0201
0202     // interface to reco jets
0203     JetContainer* jets = findNode::getClass<JetContainer>(topNode, m_recoJetName.at(input));
0204     if (!jets){
0205       std::cout << "MyJetAnalysis::process_event - Error can not find DST Reco JetContainer node " << m_recoJetName.at(input) << std::endl;
0206       exit(-1);
0207     }
0208
0209     //interface to truth jets
0210     //JetMap* jetsMC = findNode::getClass<JetMap>(topNode, m_truthJetName.at(input));
0211     JetContainer* jetsMC = findNode::getClass<JetContainer>(topNode, m_truthJetName.at(input));
0212     if (!jetsMC && m_doTruthJets){
0213       std::cout << "MyJetAnalysis::process_event - Error can not find DST Truth JetContainer node " << m_truthJetName.at(input) << std::endl;
0214       exit(-1);
0215     }
0216
0217     //get associated hepMCGenEvent
0218     HepMC::GenEvent *hepMCGenEvent = nullptr;
0219     if(m_doTruthJets){
0220       PHHepMCGenEventMap *hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0221
0222       /// If the node was not properly put on the tree, return
0223       if (!hepmceventmap){
0224         std::cout << PHWHERE
0225             << "HEPMC event map node is missing, can't collected HEPMC truth particles"
0226             << std::endl;
0227         return 0;
0228       }
0229
0230       PHHepMCGenEvent *hepmcevent = hepmceventmap->get(1);
0231
0232       if (!hepmcevent)
0233       {
0234         std::cout << PHWHERE
0235             << "PHHepMCGenEvent node is missing, can't collected HEPMC truth particles"
0236             << std::endl;
0237         return 0;
0238       }
0239
0240       hepMCGenEvent = hepmcevent->getEvent();
0241
0242       if(!hepMCGenEvent) return Fun4AllReturnCodes::ABORTEVENT;
0243     }
0244
0245     //get Particle Flow container
0246     ParticleFlowElementContainer *pflowContainer = findNode::getClass<ParticleFlowElementContainer>(topNode, "ParticleFlowElements");
0247
0248     if(!pflowContainer && m_jet_type==TYPE::FULLJET){
0249       std::cout << PHWHERE
0250           << "ParticleFlowElements node is missing, can't collect particles"
0251           << std::endl;
0252       return -1;
0253     }
0254
0255
0256     SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0257   if (!trackmap)
0258   {
0259     std::cout << PHWHERE
0260           << "SvtxTrackMap node is missing, can't collect particles"
0261           << std::endl;
0262       return -1;
0263   }
0264
0265   GlobalVertex *prim_vtx = nullptr;
0266
0267     for(auto vertex : *vertexmap){
0268     //std::cout<<"map entry"<<std::endl;
0269     PrimaryVertex primary;
0270     std::vector<GlobalVertex::VTXTYPE> source;
0271     source.clear();
0272     primary.id = vertex.second->get_id();
0273     primary.x = vertex.second->get_x();
0274     primary.y = vertex.second->get_y();
0275     primary.z = vertex.second->get_z();
0276     primary.x_unc = sqrt(vertex.second->get_error(0,0));
0277     primary.y_unc = sqrt(vertex.second->get_error(1,1));
0278     primary.z_unc = sqrt(vertex.second->get_error(2,2));
0279     primary.chisq = vertex.second->get_chisq();
0280     primary.ndf = vertex.second->get_ndof();
0281
0282      //std::cout<<std::endl<<"A "<<vertex.second->get_x()<<" "<<vertex.second->get_y()<<" "<<vertex.second->get_z()<<" "<<vertex.second->get_id()<<" "<<vertex.first<<std::endl;
0283      //for(auto vx :vertex.second->begin_vertexes()->second){
0284      //   std::cout<<"vx "<<vx->get_x()<<" "<<vx->get_y()<<" "<<vx->get_z()<<" "<<vx->get_chisq()<<vx->get_chisq()<<std::endl;
0285      //}
0286
0287     for (auto iter = vertex.second->begin_vertexes();   iter != vertex.second->end_vertexes();   ++iter){
0288       source.push_back(iter->first);
0289       GlobalVertex::VertexVector vtx = iter->second;
0290       //std::cout<<"vertex source "<<iter->first<<std::endl;
0291       if(iter->first == 400){
0292
0293         prim_vtx = vertex.second;
0294       }
0295
0296
0297       //for(auto vx :vtx){
0298       //  std::cout<<"vx "<<vx->get_x()<<" "<<vx->get_y()<<" "<<vx->get_z()<<" "<<vx->get_chisq()<<std::endl;
0299
0300
0301       //}
0302     }
0303
0304   if(std::find(source.begin(), source.end(), GlobalVertex::VTXTYPE::SVTX) == source.end()) return Fun4AllReturnCodes::ABORTEVENT;
0305   m_stat->Fill(3);
0306
0307     //if((std::find(source.begin(), source.end(), GlobalVertex::VTXTYPE::SVTX) != source.end()) && (std::find(source.begin(), source.end(), GlobalVertex::VTXTYPE::MBD) != source.end())) primary.vtxtype = GlobalVertex::VTXTYPE::SVTX_MBD;
0308     if(std::find(source.begin(), source.end(), GlobalVertex::VTXTYPE::SVTX) != source.end()) primary.vtxtype = GlobalVertex::VTXTYPE::SVTX;
0309     else if(std::find(source.begin(), source.end(), GlobalVertex::VTXTYPE::MBD) != source.end()) primary.vtxtype = GlobalVertex::VTXTYPE::MBD;
0310     m_container[input]->primaryVertex = primary;
0311   }
0312
0313   if(prim_vtx->get_z() > 10 || prim_vtx->get_z() < -10) return Fun4AllReturnCodes::ABORTEVENT;
0314  m_stat->Fill(4);
0315
0316
0317
0318     int nrecojet = -1;
0319
0320     Jet::PROPERTY recojet_area_index = jets->property_index(Jet::PROPERTY::prop_area);
0321
0322     for (Jet* jet : *jets){
0323       if(jet->get_pt() < m_ptRangeReco.first || jet->get_pt() > m_ptRangeReco.second) continue;
0324       if (not (std::abs(jet->get_eta()) <= 1.1 - (doFiducial?jetR.at(input):0))) continue;
0325
0326       //std::cout<<"jet pt "<<jet->get_pt()<<std::endl;
0327       nrecojet++;
0328
0329       int nChtracks = 0;
0330
0331       RecoJets recojet;
0332
0333       //loop over jet constituents
0334       for (const auto& comp : jet->get_comp_vec()){
0335         ParticleFlowElement *pflow = nullptr;
0336
0337         SvtxTrack *trk = nullptr;
0338
0339         if(m_jet_type==TYPE::FULLJET){
0340           pflow = pflowContainer->getParticleFlowElement(comp.second);
0341           trk = pflow->get_track();
0342         }
0343         else if(m_jet_type==TYPE::CHARGEJET){
0344           trk = trackmap->get(comp.second);
0345         }
0346
0347          //for (SvtxTrackMap::ConstIter iter = trackmap->begin(); iter != trackmap->end(); ++iter){
0348             //const SvtxTrack *track = iter->second;
0349
0350             //std::cout<<iter->first<<" "<<track->get_id()<<std::endl;
0351           //}
0352
0353         //if charged track
0354         if(trk){
0355           chConstituent track_properties;
0356
0357              int n_mvtx_hits = 0;
0358           int n_intt_hits = 0;
0359           int n_tpc_hits = 0;
0360 //std::cout<<"C"<<std::endl;
0361           for (const auto& ckey : TrackAnalysisUtils::get_cluster_keys(trk)){
0362             switch (TrkrDefs::getTrkrId(ckey)){
0363               case TrkrDefs::mvtxId:
0364                 n_mvtx_hits++;
0365                 break;
0366               case TrkrDefs::inttId:
0367                 n_intt_hits++;
0368                 break;
0369               case TrkrDefs::tpcId:
0370                 n_tpc_hits++;
0371                 break;
0372             }
0373           }
0374
0375           //std::cout << "x:"<<trk->get_x()<< " y:"<<trk->get_y()<< " z:"<<trk->get_z()<< " eta:"<<trk->get_eta()<< " pt:"<<trk->get_pt()<< " chi2:"<<trk->get_chisq()/trk->get_ndf()<<" "<<n_mvtx_hits<<" "<<n_intt_hits<<" "<<n_tpc_hits<<"charge: "<<trk->get_charge();
0376
0377           //get vertex associated to tarck
0378           int id = trk->get_vertex_id();
0379           //std::cout<<id<<std::endl;
0380           GlobalVertex *vtx = prim_vtx;//vertexmap->get(id);
0381           //skip if track without vertex
0382
0383
0384           float DCA_xy = -999;
0385           float DCA_xy_unc = -999;
0386           float DCA_3d = -999;
0387           float chi2_3d = 1;
0388           double sign = -999;
0389           double sign_3d= -999;
0390
0391           if (vtx == nullptr){
0392             std::cout << "MyJetAnalysis::process_event - Track does not have assigned vertex" << std::endl;
0393             //std::cout << "x:"<<pflow->get_track()->get_x()<< " y:"<<pflow->get_track()->get_y()<< " z:"<<pflow->get_track()->get_z()<< " eta:"<<pflow->get_track()->get_eta()<< " pt:"<<pflow->get_track()->get_pt()<<std::endl;
0394             //continue;
0395           } else if(n_mvtx_hits > 0 && n_intt_hits > 0){
0396               GetDistanceFromVertex(trk,vtx,DCA_xy,DCA_xy_unc,DCA_3d,chi2_3d);
0397               double dot = (trk->get_x() - vtx->get_x()) * jet->get_px() + (trk->get_y() - vtx->get_y()) * jet->get_py();
0398               sign = int(dot/std::abs(dot));
0399
0400               double dot_3d = (trk->get_x() - vtx->get_x()) * jet->get_px() + (trk->get_y() - vtx->get_y()) * jet->get_py() + (trk->get_z() - vtx->get_z()) * jet->get_pz();
0401               sign_3d = int(dot_3d/std::abs(dot_3d));
0402 //std::cout<<"DCA3d "<<DCA_3d<<" chi23d "<<chi2_3d<<" dot_3d "<<dot_3d<<" sign "<<sign_3d<<std::endl;
0403               //std::cout<<" signdca:"<< sign*std::abs(DCA_xy/DCA_xy_unc)<<std::endl;
0404           }
0405           nChtracks++;
0406 //std::cout<<"B"<<std::endl;
0407
0408            //rerutn DCA_XY vector, val ||DCA|| + unc DCA
0409
0410
0411           //tarcking team way - no uncertainty of vertex in DCA calculation
0412           /*Acts::Vector3 vtxActs(vtx->get_x(), vtx->get_y(), vtx->get_z());
0413           std::pair<std::pair<float, float>, std::pair<float, float>> DCApair;
0414           DCApair = TrackAnalysisUtils::get_dca(pflow->get_track(),vtxActs);
0415           double dot2 = (pflow->get_track()->get_x() - vtx->get_x()) * jet->get_px() + (pflow->get_track()->get_y() - vtx->get_y()) * jet->get_py();
0416           double sign2 = int(dot2/std::abs(dot2));*/
0417
0418           //get some track quality values
0419
0420 //std::cout<<"D"<<std::endl;
0421
0422           /*m_chi2ndf[input]->Fill(pflow->get_track()->get_chisq()/pflow->get_track()->get_ndf());
0423           m_mvtxcl[input]->Fill(m_nmaps);
0424           m_inttcl[input]->Fill(m_nintt);
0425           m_mtpccl[input]->Fill(m_ntpc);*/
0426           if(m_jet_type==TYPE::FULLJET)track_properties.pflowtype = pflow->get_type();
0427           track_properties.vtx_id = id;
0428           if(m_jet_type==TYPE::FULLJET)track_properties.e = pflow->get_e();
0429           track_properties.eta = trk->get_eta();
0430           track_properties.phi = trk->get_phi();
0431           track_properties.pt = trk->get_pt();
0432           track_properties.charge = trk->get_charge();
0433           track_properties.DCA_xy = DCA_xy;
0434           track_properties.DCA_xy_unc = DCA_xy_unc;
0435           track_properties.sDCA_xy = sign*std::abs(DCA_xy/DCA_xy_unc);
0436           track_properties.DCA3d = DCA_3d;
0437 //std::cout<<"final "<<sign_3d<<" "<<chi2_3d<<" "<<sign_3d*std::sqrt(chi2_3d)<<std::endl;
0438           track_properties.sDCA3d = sign_3d*std::sqrt(chi2_3d);
0439           track_properties.n_mvtx = n_mvtx_hits;
0440           track_properties.n_intt = n_intt_hits;
0441           track_properties.n_tpc = n_tpc_hits;
0442           track_properties.chisq = trk->get_chisq();
0443           track_properties.ndf = trk->get_ndf();
0444 //std::cout<<"E"<<std::endl;
0445           recojet.chConstituents.push_back(track_properties);
0446         } // end if(pflow->get_track())
0447         else if(m_jet_type==TYPE::FULLJET){ //neutral tracl
0448           neConstituent neutral_properties;
0449           neutral_properties.pflowtype = pflow->get_type();
0450           neutral_properties.e = pflow->get_e();
0451           neutral_properties.eta = pflow->get_eta();
0452           neutral_properties.phi = pflow->get_phi();
0453           recojet.neConstituents.push_back(neutral_properties);
0454         }
0455       } // end for (const auto& comp : jet->get_comp_vec())
0456 //std::cout<<"F"<<std::endl;
0457       //recojet.id = jet->get_id();
0458       recojet.id = nrecojet;
0459       recojet.area = jet->get_property(recojet_area_index);
0460       recojet.num_Constituents = static_cast<int>(jet->get_comp_vec().size());
0461       recojet.num_ChConstituents = nChtracks;
0462       recojet.px = jet->get_px();
0463       recojet.py = jet->get_py();
0464       recojet.pz = jet->get_pz();
0465       recojet.pt = jet->get_pt();
0466       recojet.eta = jet->get_eta();
0467       recojet.phi = jet->get_phi();
0468       recojet.m = jet->get_mass();
0469       recojet.e = jet->get_e();
0470 //std::cout<<"G"<<std::endl;
0471       m_container[input]->recojets.push_back(recojet);
0472     } // end for (Jet* jet : *jets)
0473     m_container[input]->reco_jet_n = static_cast<int>(nrecojet+1);
0474
0475 //std::cout<<"H"<<std::endl;
0476
0477     //get truth jets
0478     if(m_doTruthJets){
0479       int ntruthjet = -1;
0480       TruthJets mtruthjet;
0481       Jet::PROPERTY truthjet_area_index = jetsMC->property_index(Jet::PROPERTY::prop_area);
0482
0483       //Jet::PROPERTY recojet_area_index = jetsMC->property_index(Jet::PROPERTY::prop_area);
0484
0485       //truth jet loop
0486       //for (JetMap::Iter iter = jetsMC->begin(); iter != jetsMC->end(); ++iter){
0487
0488       for (Jet* truthjet : *jetsMC){
0489
0490         //Jet* truthjet = iter->second;
0491         if(truthjet->get_pt() < m_ptRangeTruth.first || truthjet->get_pt() > m_ptRangeTruth.second) continue;
0492         if (not (std::abs(truthjet->get_eta()) <= 1.1 - (doFiducial?jetR.at(input):0))) continue;
0493         if(truthjet->get_pt() < 10.0) continue;
0494         ntruthjet++;
0495         //std::cout<<"new Truth Jet "<<std::endl;
0496
0497
0498         //mtruthjet.id = truthjet->get_id();
0499         mtruthjet.id = ntruthjet;
0500         mtruthjet.area = truthjet->get_property(truthjet_area_index);
0501         mtruthjet.num_Constituents = truthjet->size_comp();
0502         mtruthjet.px = truthjet->get_px();
0503         mtruthjet.py = truthjet->get_py();
0504         mtruthjet.pz = truthjet->get_pz();
0505         mtruthjet.pt = truthjet->get_pt();
0506         mtruthjet.eta = truthjet->get_eta();
0507         mtruthjet.phi = truthjet->get_phi();
0508         mtruthjet.m = truthjet->get_mass();
0509         mtruthjet.e = truthjet->get_e();
0510
0511
0512         //jet-tarck PDG identification
0513         TString taggedPDGIDs[] = {"B-Meson","CharmedMeson","CharmedBaryon", "B-Baryon"};
0514         int forbiddenPDGIDs[] = {1,2,3,4,5,6,7,8,21,22};
0515         int nChTrackstruth = 0;
0516
0517         //loop over jet constituents
0518         //for(auto citer = truthjet->begin_comp(); citer != truthjet->end_comp(); ++citer){
0519
0520         for (const auto& comp : truthjet->get_comp_vec()){
0521           //std::cout<<"constituent: "<<std::endl;
0522           //PHG4Particle *g4part = truthinfo->GetPrimaryParticle(citer->second);
0523           PHG4Particle *g4part = truthinfo->GetPrimaryParticle(comp.second);
0524             HepMC::GenParticle* constituent = hepMCGenEvent->barcode_to_particle(g4part->get_barcode());
0525           if(std::abs((TDatabasePDG::Instance()->GetParticle((constituent)->pdg_id()))->Charge()) > 10e-4) nChTrackstruth++;
0526
0527           //mother track tagged, stop search
0528           if (std::find(std::begin(taggedPDGIDs), std::end(taggedPDGIDs), TString((TDatabasePDG::Instance()->GetParticle((constituent)->pdg_id()))->ParticleClass())) != std::end(taggedPDGIDs)){
0529             //std::cout<<"mother tagged, pushing pdg id: "<<(constituent)->pdg_id()<<std::endl;
0530               mtruthjet.constituents_PDG_ID.push_back((constituent)->pdg_id());
0531               //continue; keep searching
0532           }
0533           //std::cout<<"recursive search: "<<std::endl;
0534           int i = -1;
0535           //not yet tagged, begin history search
0536             while (!constituent->is_beam()){
0537             i++;
0538             //std::cout<<"i: "<<i<<std::endl;
0539             bool breakOut = false;
0540             bool taggedHF = false;
0541             for (HepMC::GenVertex::particle_iterator mother = constituent->production_vertex()->particles_begin(HepMC::parents); mother != constituent->production_vertex()->particles_end(HepMC::parents); ++mother){
0542               //std::cout<<"inside for"<<std::endl;
0543               //found HF in parent tree
0544               if (std::find(std::begin(taggedPDGIDs), std::end(taggedPDGIDs), TString((TDatabasePDG::Instance()->GetParticle((*mother)->pdg_id()))->ParticleClass())) != std::end(taggedPDGIDs)){
0545                 //taggedHF = true;
0546                 //std::cout<<"taggedHF, pushing pdg id: "<<(*mother)->pdg_id()<<std::endl;
0547                 mtruthjet.constituents_PDG_ID.push_back((*mother)->pdg_id());
0548                 //break;
0549               }
0550               //reached partonic level, break search
0551               if (std::find(std::begin(forbiddenPDGIDs), std::end(forbiddenPDGIDs), abs((*mother)->pdg_id())) != std::end(forbiddenPDGIDs)){
0552                 //std::cout<<"partonic reached, pdg id: "<<(*mother)->pdg_id()<<std::endl;
0553                 //(*mother)->production_vertex()->print();
0554                 /*for (HepMC::GenVertex::particle_iterator mother2 = (*mother)->production_vertex()->particles_begin(HepMC::parents); mother2 != constituent->production_vertex()->particles_end(HepMC::parents); ++mother2){
0555                     (*mother2)->production_vertex()->print();
0556                 }*/
0557                 if(abs((*mother)->pdg_id()) == 4 || abs((*mother)->pdg_id()) == 5){
0558                   quark HFquark;
0559                   HFquark.vtx_barcode = (*mother)->production_vertex()->barcode();
0560                   HFquark.pdgid = (*mother)->pdg_id();
0561                   HFquark.px = (*mother)->momentum().px();
0562                   HFquark.py = (*mother)->momentum().py();
0563                   HFquark.pz = (*mother)->momentum().pz();
0564                   HFquark.e = (*mother)->momentum().e();
0565                   mtruthjet.constituents_origin_quark.push_back(HFquark);
0566                 //(*mother)->production_vertex()->print();
0567                 //(*mother)->print();
0568                 //std::cout<<"v "<<(*mother)->production_vertex()->barcode()<<" "<<(*mother)->momentum().e()<<" "<<(*mother)->momentum().px()<<" "<<(*mother)->momentum().py()<<" "<<(*mother)->momentum().pz()<<std::endl;
0569                 }
0570                 breakOut = true;
0571                 break;
0572               }
0573               //if(constituent->is_beam()) break;
0574               //std::cout<<"no HF or partonic, pdgid: "<<(*mother)->pdg_id()<<" continue search in history"<<std::endl;
0575
0576               constituent = *mother;
0577             }
0578             //if (taggedHF)std::cout<<"TAGGED HF"<<std::endl;
0579             if (breakOut || taggedHF){
0580               //std::cout<<"breaking breakOut: "<<(breakOut?"true":"false")<<" taggedHF "<<(taggedHF?"true":"false")<<std::endl;
0581               break;
0582             }
0583           }
0584         }// end loop over constituents
0585
0586         mtruthjet.num_ChConstituents = nChTrackstruth ;
0587         m_container[input]->truthjets.push_back(mtruthjet);
0588       }// end jet loop
0589       m_container[input]->truth_jet_n = static_cast<int>(ntruthjet+1);
0590     }//end if do truth
0591     //fill tree
0592     m_T[input]->Fill();
0593   }//end input loop
0594
0595   return Fun4AllReturnCodes::EVENT_OK;
0596 }
0597
0598 //____________________________________________________________________________..
0599 int FullJetFinder::ResetEvent(PHCompositeNode *topNode)
0600 {
0601   //std::cout << "FullJetFinder::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0602   return Fun4AllReturnCodes::EVENT_OK;
0603 }
0604
0605 //____________________________________________________________________________..
0606 int FullJetFinder::EndRun(const int runnumber)
0607 {
0608   std::cout << "FullJetFinder::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0609   return Fun4AllReturnCodes::EVENT_OK;
0610 }
0611
0612 //____________________________________________________________________________..
0613 int FullJetFinder::End(PHCompositeNode *topNode)
0614 {
0615   std::cout << "FullJetFinder::End - Output to " << m_outputFileName << std::endl;
0616
0617    /*CutSelection cutselection;
0618     cutselection.jetptmin = 5;
0619     cutselection.jetptmin = 30;*/
0620
0621   if(PHTFileServer::get().cd(m_outputFileName)){
0622     m_stat->Write();
0623   }
0624
0625   for(int i = 0 ; i < m_inputs; i++){
0626     if(m_T[i] ){
0627       if(PHTFileServer::get().cd(m_outputFileName)){
0628         TDirectory *cdtof = gDirectory->mkdir(m_TreeNameCollection.at(i).data());
0629         cdtof->cd();
0630         //cdtof-><CutSelection>WriteObject(cutselection, "Selection");
0631         //cutselection.Write("Selection");
0632         m_T[i]->Write();
0633       }
0634     }
0635   }
0636   std::cout << "FullJetFinder::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0637   return Fun4AllReturnCodes::EVENT_OK;
0638 }
0639
0640 //____________________________________________________________________________..
0641 int FullJetFinder::Reset(PHCompositeNode *topNode)
0642 {
0643  std::cout << "FullJetFinder::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0644   return Fun4AllReturnCodes::EVENT_OK;
0645 }
0646
0647 //____________________________________________________________________________..
0648 void FullJetFinder::Print(const std::string &what) const
0649 {
0650   std::cout << "FullJetFinder::Print(const std::string &what) const Printing info for " << what << std::endl;
0651 }
0652
0653 void FullJetFinder::GetDistanceFromVertex(SvtxTrack *m_dst_track, GlobalVertex *m_dst_vertex,float &val_xy, float &err_xy, float &val_3d, float &chi2_3d){
0654
0655   KFParticle::SetField(1.4e1);
0656
0657   KFParticle kfp_particle;
0658
0659   float f_trackParameters[6] = {m_dst_track->get_x(),
0660                                 m_dst_track->get_y(),
0661                                 m_dst_track->get_z(),
0662                                 m_dst_track->get_px(),
0663                                 m_dst_track->get_py(),
0664                                 m_dst_track->get_pz()};
0665
0666   float f_trackCovariance[21];
0667   unsigned int iterate = 0;
0668   for (unsigned int i = 0; i < 6; ++i)
0669     for (unsigned int j = 0; j <= i; ++j)
0670     {
0671       f_trackCovariance[iterate] = m_dst_track->get_error(i, j);
0672       ++iterate;
0673     }
0674
0675   kfp_particle.Create(f_trackParameters, f_trackCovariance, (Int_t) m_dst_track->get_charge(), -1);
0676   kfp_particle.NDF() = m_dst_track->get_ndf();
0677   kfp_particle.Chi2() = m_dst_track->get_chisq();
0678   kfp_particle.SetId(m_dst_track->get_id());
0679
0680   float f_vertexParameters[6] = {m_dst_vertex->get_x(),
0681                                  m_dst_vertex->get_y(),
0682                                  m_dst_vertex->get_z(), 0, 0, 0};
0683
0684   float f_vertexCovariance[21];
0685   unsigned int iterate2 = 0;
0686   for (unsigned int i = 0; i < 3; ++i)
0687     for (unsigned int j = 0; j <= i; ++j)
0688     {
0689       f_vertexCovariance[iterate2] = m_dst_vertex->get_error(i, j);
0690       ++iterate2;
0691     }
0692
0693   KFParticle kfp_vertex;
0694   kfp_vertex.Create(f_vertexParameters, f_vertexCovariance, 0, -1);
0695   kfp_vertex.NDF() = m_dst_vertex->get_ndof();
0696   kfp_vertex.Chi2() = m_dst_vertex->get_chisq();
0697   kfp_particle.GetDistanceFromVertexXY(kfp_vertex,val_xy,err_xy);
0698   val_3d = kfp_particle.GetDistanceFromVertex(kfp_vertex);
0699   chi2_3d = kfp_particle.GetDeviationFromVertex(kfp_vertex);
0700 }
0701