Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:12:12 

0001 #include "ElectronID.h"
0002 
0003 #include <cstdlib>
0004 #include <cstdio>
0005 #include <iomanip>
0006 #include <fstream>
0007 
0008 #include <phool/PHCompositeNode.h>
0009 #include <phool/recoConsts.h>
0010 #include <phool/PHIODataNode.h>                         // for PHIODataNode
0011 #include <phool/PHNode.h>                               // for PHNode
0012 #include <phool/PHNodeIterator.h>
0013 #include <phool/PHObject.h>                             // for PHObject
0014 #include <phool/getClass.h>
0015 #include <phool/phool.h>  // for PHWHERE
0016 #include <phool/PHRandomSeed.h>
0017 #include <phool/getClass.h>
0018 
0019 #include <fun4all/Fun4AllReturnCodes.h>
0020 
0021 #include <trackbase_historic/SvtxTrackMap.h>
0022 #include <trackbase_historic/SvtxTrack.h>
0023 #include <trackbase_historic/SvtxVertex.h>
0024 #include <trackbase_historic/SvtxVertexMap.h>
0025 #include <trackbase/TrkrDefs.h>
0026 
0027 #include <calobase/RawClusterContainer.h>
0028 #include <calobase/RawCluster.h>
0029 #include <calobase/RawClusterv1.h>
0030 
0031 #include <globalvertex/GlobalVertexMap.h>
0032 #include <globalvertex/GlobalVertex.h>
0033 
0034 #include <g4main/PHG4TruthInfoContainer.h>
0035 #include <g4main/PHG4Particle.h>
0036 #include <g4main/PHG4VtxPoint.h>
0037 
0038 #include "trackpidassoc/TrackPidAssoc.h"
0039 
0040 //TMVA class
0041 #include <vector>
0042 #include <iostream>
0043 #include <map>
0044 #include <string>
0045 #include "TMVA/Tools.h"
0046 #include "TMVA/Reader.h"
0047 #include "TMVA/MethodCuts.h"
0048 #include "TMVA/Factory.h"
0049 #include "TMVA/DataLoader.h"
0050 
0051 
0052 // gsl
0053 #include <gsl/gsl_randist.h>
0054 #include <gsl/gsl_rng.h>
0055 
0056 #include <TVector3.h>
0057 #include <TMatrixFfwd.h>    // for TMatrixF
0058 
0059 #include <iostream>                            // for operator<<, basic_ostream
0060 #include <set>                                 // for _Rb_tree_iterator, set
0061 #include <utility>                             // for pair
0062 
0063 class PHCompositeNode;
0064 
0065 using namespace std;
0066 
0067 ElectronID::ElectronID(const std::string& name, const std::string &filename) : SubsysReco(name)
0068 {
0069   OutputNtupleFile=nullptr;
0070   OutputFileName=filename;
0071   EventNumber=0;
0072   output_ntuple = true;
0073 
0074   /// default limits 
0075   EMOP_lowerlimit = 0.7;
0076   EMOP_higherlimit = 100.0;
0077   HOP_lowerlimit = 0.0;
0078   HinOEM_higherlimit = 100.0;
0079   Pt_lowerlimit = 0.0;
0080   Pt_higherlimit = 100.0;
0081   Nmvtx_lowerlimit = 2;
0082   Nintt_lowerlimit = 0;
0083   Ntpc_lowerlimit = 0;
0084   Nquality_higherlimit = 100;
0085   Ntpc_lowerlimit = 20;
0086   Nquality_higherlimit = 5.;
0087   PROB_cut = 0.;
0088 
0089   /// MVA
0090   BDT_cut_p = 0.0;
0091   BDT_cut_n = 0.0;
0092   ISUSE_BDT_p =0;
0093   ISUSE_BDT_n =0;
0094 
0095  // unsigned int _nlayers_maps = 3;
0096  // unsigned int _nlayers_intt = 4;
0097  // unsigned int _nlayers_tpc = 48;
0098 }
0099 
0100 ElectronID::~ElectronID() 
0101 {
0102 }
0103 
0104 int ElectronID::Init(PHCompositeNode *topNode)
0105 {
0106 
0107   if(output_ntuple) {
0108 
0109     OutputNtupleFile = new TFile(OutputFileName.c_str(),"RECREATE");
0110     std::cout << "PairMaker::Init: output file " << OutputFileName.c_str() << " opened." << endl;
0111         ntpBDTresponse = new TNtuple("ntpbeforecut","","select_p:select_n");
0112 
0113     ntpbeforecut = new TNtuple("ntpbeforecut","","p:pt:cemce3x3:hcaline3x3:hcaloute3x3:cemce3x3overp:hcaline3x3overcemce3x3:hcale3x3overp:charge:pid:quality:e_cluster:EventNumber:z:vtxid:nmvtx:nintt:ntpc:cemc_prob:cemc_ecore:cemc_chi2");
0114         ntpcutEMOP = new TNtuple("ntpcutEMOP","","p:pt:cemce3x3:hcaline3x3:hcaloute3x3:charge:pid:quality:e_cluster:EventNumber:z:vtxid:cemc_prob:cemc_ecore");
0115     ntpcutHOP = new TNtuple("ntpcutHOP","","p:pt:cemce3x3:hcaline3x3:hcaloute3x3:charge:pid:quality:e_cluster:EventNumber:z:vtxid:cemc_prob:cemc_ecore");
0116     ntpcutEMOP_HinOEM = new TNtuple("ntpcutEMOP_HinOEM","","p:pt:cemce3x3:hcaline3x3:hcaloute3x3:charge:pid:quality:e_cluster:EventNumber:z:vtxid:cemc_prob:cemc_ecore");
0117     ntpcutEMOP_HinOEM_Pt = new TNtuple("ntpcutEMOP_HinOEM_Pt","","p:pt:cemce3x3:hcaline3x3:hcaloute3x3:charge:pid:quality:e_cluster:EventNumber:z:vtxid:cemc_prob:cemc_ecore");
0118         ntpcutEMOP_HinOEM_Pt_read = new TNtuple("ntpcutEMOP_HinOEM_Pt_read","","p:pt:cemce3x3:hcaline3x3:hcaloute3x3:charge:pid:quality:e_cluster:EventNumber:z:vtxid:trackid:cemc_prob:cemc_ecore");
0119         ntpcutBDT_read = new TNtuple("ntpcutBDT_read","","p:pt:cemce3x3:hcaline3x3:hcaloute3x3:charge:pid:quality:e_cluster:EventNumber:z:vtxid:trackid:cemc_prob:cemc_ecore:EOP:HOM:cemc_chi2");
0120   }
0121   else {
0122     PHNodeIterator iter(topNode);
0123     PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0124     if (!dstNode)
0125     {
0126       cerr << PHWHERE << " ERROR: Can not find DST node." << endl;
0127           return Fun4AllReturnCodes::ABORTEVENT;
0128     }
0129   }
0130   
0131   return Fun4AllReturnCodes::EVENT_OK;
0132   
0133 }
0134 
0135 int ElectronID::InitRun(PHCompositeNode* topNode)
0136 {
0137   int ret = GetNodes(topNode);
0138   if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
0139 
0140   return Fun4AllReturnCodes::EVENT_OK;
0141 }
0142 
0143 int ElectronID::process_event(PHCompositeNode* topNode)
0144 {
0145   EventNumber++;
0146   float ntp[30];
0147 
0148   SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0149   if(!trackmap) {
0150     cerr << PHWHERE << " ERROR: Can not find SvtxTrackMap node." << endl;
0151     return Fun4AllReturnCodes::ABORTEVENT;
0152   }
0153 
0154   cout<<"EventNumber ===================== " << EventNumber-1 << endl;
0155   if(EventNumber==1) topNode->print();
0156 
0157   SvtxVertexMap *vtxmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0158   if(!vtxmap) {
0159       cout << "SvtxVertexMap node not found!" << endl;
0160       return Fun4AllReturnCodes::ABORTEVENT;
0161   }
0162   //cout << "Number of SvtxVertexMap entries = " << vtxmap->size() << endl;
0163 
0164   RawClusterContainer* cemc_cluster_container = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
0165   if(!cemc_cluster_container) {
0166     cerr << PHWHERE << " ERROR: Can not find CLUSTER_CEMC node." << endl;
0167     return Fun4AllReturnCodes::ABORTEVENT;
0168   }
0169 
0170 
0171 // Truth info
0172   PHG4TruthInfoContainer* truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0173   if(!truth_container) {
0174     cerr << PHWHERE << " ERROR: Can not find G4TruthInfo node." << endl;
0175     return Fun4AllReturnCodes::ABORTEVENT;
0176   }
0177   PHG4TruthInfoContainer::ConstRange range = truth_container->GetPrimaryParticleRange();
0178   cout << "number of MC particles = " << truth_container->size() << " " << truth_container->GetNumPrimaryVertexParticles() << endl;
0179 
0180     int mycount = 0;
0181     for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0182     {
0183       PHG4Particle* g4particle = iter->second;
0184       mycount++;
0185       int gflavor  = g4particle->get_pid();
0186       double gpx = g4particle->get_px();
0187       double gpy = g4particle->get_py();
0188       double gpz= g4particle->get_pz();
0189       double gpt = sqrt(gpx*gpx+gpy*gpy);
0190       double phi = atan2(gpy,gpx);
0191       double eta = asinh(gpz/gpt);
0192       int primid =  g4particle->get_primary_id();
0193       int parentid = g4particle->get_parent_id();
0194       int trackid = g4particle->get_track_id();
0195       if(trackid>truth_container->GetNumPrimaryVertexParticles()-50) cout << trackid << " " << parentid << " " << primid << " " << gflavor << " " << gpt << " " << phi << " " << eta << endl;
0196     }
0197     cout << "mycount = " << mycount << endl;
0198 // end Truth
0199 
0200 
0201 //MVA method setup**********************
0202 
0203   // This loads the library
0204    TMVA::Tools::Instance();
0205 
0206    // Default MVA methods with trained + tested
0207    std::map<std::string,int> Use_p;
0208    std::map<std::string,int> Use_n;
0209    // Boosted Decision Trees
0210    Use_p["BDT"] = ISUSE_BDT_p; // uses Adaptive Boost; for positive particles
0211    Use_n["BDT"] = ISUSE_BDT_n; // uses Adaptive Boost; for negative particles
0212 
0213   // Create the Reader object
0214 
0215    TMVA::Reader *reader_positive = new TMVA::Reader( "!Color:!Silent" );
0216    TMVA::Reader *reader_negative = new TMVA::Reader( "!Color:!Silent" );
0217 
0218    // Create a set of variables and declare them to the reader
0219    // - the variable names MUST corresponds in name and type to those given in the weight file(s) used
0220    Float_t var1_p, var2_p, var3_p;
0221    Float_t var1_n, var2_n, var3_n;
0222    reader_positive->AddVariable( "var1_p", &var1_p );
0223    reader_positive->AddVariable( "var2_p", &var2_p );
0224    reader_positive->AddVariable( "var3_p", &var3_p );
0225 
0226    reader_negative->AddVariable( "var1_n", &var1_n );
0227    reader_negative->AddVariable( "var2_n", &var2_n );
0228    reader_negative->AddVariable( "var3_n", &var3_n );
0229 
0230   // Book the MVA methods
0231    TString dir_p, dir_n;
0232    dir_p = "dataset/Weights_positive/"; // weights for positive particles
0233    dir_n = "dataset/Weights_negative/"; // weights for negative particles
0234    TString prefix = "TMVAClassification";
0235 
0236    // Book method(s)
0237    for (std::map<std::string,int>::iterator it = Use_p.begin(); it != Use_p.end(); it++) {
0238       if (it->second) {
0239          TString methodName_p = TString(it->first) + TString(" method");
0240          TString weightfile_p = dir_p + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
0241          reader_positive->BookMVA( methodName_p, weightfile_p );
0242       }
0243    }
0244    
0245    for (std::map<std::string,int>::iterator it = Use_n.begin(); it != Use_n.end(); it++) {
0246       if (it->second) {
0247          TString methodName_n = TString(it->first) + TString(" method");
0248          TString weightfile_n = dir_n + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
0249          reader_negative->BookMVA( methodName_n, weightfile_n );
0250       }
0251    }
0252 
0253 //end of MVA method setup*******************
0254 
0255   int nmvtx = 0;
0256   int nintt = 0;
0257   int ntpc = 0;
0258   cout << _nlayers_maps << " " << _nlayers_intt << " " << _nlayers_tpc<<endl;
0259   // get the tracks
0260   for(SvtxTrackMap::Iter it = _track_map->begin(); it != _track_map->end(); ++it)
0261     {
0262       SvtxTrack *track = it->second;
0263 
0264       PHG4Particle* g4particle = findMCmatch(track, truth_container);
0265       int gflavor = g4particle->get_pid();
0266       if(gflavor==0) continue;
0267 
0268       nmvtx = 0;
0269       nintt = 0;
0270       ntpc = 0;
0271 
0272       for (SvtxTrack::ConstClusterKeyIter iter = track->begin_cluster_keys(); iter != track->end_cluster_keys(); ++iter)
0273       {
0274         TrkrDefs::cluskey cluser_key = *iter;
0275         int trackerid = TrkrDefs::getTrkrId(cluser_key);
0276       //  cout << "trackerid= " << trackerid << endl; 
0277         if(trackerid==0) nmvtx++;
0278         if(trackerid==1) nintt++;
0279         if(trackerid==2) ntpc++;
0280         //unsigned int layer = TrkrDefs::getLayer(cluser_key);
0281         //if (_nlayers_maps > 0 && layer < _nlayers_maps) nmvtx++;
0282         //if (_nlayers_intt > 0 && layer >= _nlayers_maps && layer < _nlayers_maps + _nlayers_intt) nintt++;
0283         //if (_nlayers_tpc > 0 && layer >= (_nlayers_maps + _nlayers_intt) && layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc)) ntpc++;
0284       }
0285  
0286 
0287       double px = track->get_px();
0288       double py = track->get_py();
0289 
0290       double z = track->get_z();
0291 
0292       unsigned int vtxid = track->get_vertex_id();
0293 
0294       double mom = track->get_p();
0295       double pt = sqrt(px*px + py*py);
0296       int charge = track->get_charge();
0297       int pid = it->first;
0298       float quality = track->get_quality();
0299 
0300       double e_cluster = track->get_cal_cluster_e(SvtxTrack::CAL_LAYER::CEMC);
0301 
0302       double e_cemc_3x3 = track->get_cal_energy_3x3(SvtxTrack::CAL_LAYER::CEMC);
0303       double e_hcal_in_3x3 = track->get_cal_energy_3x3(SvtxTrack::CAL_LAYER::HCALIN);
0304       double e_hcal_out_3x3 = track->get_cal_energy_3x3(SvtxTrack::CAL_LAYER::HCALOUT);
0305 
0306       unsigned int cemc_clusid = track->get_cal_cluster_id(SvtxTrack::CAL_LAYER::CEMC);
0307       double cemc_prob = 0.;
0308       double cemc_ecore = 0.;
0309       double cemc_chi2 = 99999.;
0310       if(cemc_clusid<99999) {
0311         RawCluster* cemc_cluster = cemc_cluster_container->getCluster(cemc_clusid);
0312         cemc_prob = cemc_cluster->get_prob();
0313         cemc_ecore = cemc_cluster->get_ecore();
0314         cemc_chi2 = cemc_cluster->get_chi2();
0315       }
0316 
0317       // CEMC E/p cut
0318       double cemceoverp = e_cemc_3x3 / mom;
0319       //double cemceoverp = cemc_ecore / mom;
0320       // HCaline/CEMCe cut
0321       double hcalineovercemce = e_hcal_in_3x3 / e_cemc_3x3;
0322       // HCal E/p cut
0323       double hcaleoverp = (e_hcal_in_3x3 + e_hcal_out_3x3) / mom;
0324 
0325     //MVA valuables
0326      if(cemceoverp>0.0 && cemceoverp<10.0 && hcalineovercemce>0.0 && hcalineovercemce<10.0){
0327          if(charge>0){
0328              var1_p = cemceoverp;
0329              var2_p = hcalineovercemce;
0330              var3_p = cemc_chi2;
0331           }
0332          if(charge<0){
0333              var1_n = cemceoverp;
0334              var2_n = hcalineovercemce;
0335              var3_n = cemc_chi2;
0336           }
0337       }
0338     
0339     if (Use_p["BDT"] && Use_n["BDT"] && quality < Nquality_higherlimit && nmvtx >= Nmvtx_lowerlimit && nintt >= Nintt_lowerlimit && ntpc >= Ntpc_lowerlimit && pt > Pt_lowerlimit && pt < Pt_higherlimit &&cemceoverp>0.0 && cemceoverp<10.0 && hcalineovercemce>0.0 && hcalineovercemce<10.0) {
0340           float select_p=reader_positive->EvaluateMVA("BDT method");
0341       float select_n=reader_negative->EvaluateMVA("BDT method");
0342           ntp[0] = select_p;
0343           ntp[1] = select_n;
0344       if(output_ntuple) { ntpBDTresponse -> Fill(ntp); }
0345           if(select_p>BDT_cut_p && select_n>BDT_cut_n){
0346               // add to the association map
0347           _track_pid_assoc->addAssoc(TrackPidAssoc::electron, it->second->get_id());
0348           }
0349 
0350     }// end of BDT cut
0351     else{ // for traditional cuts
0352 
0353       ntp[0] = mom;
0354       ntp[1] = pt;
0355       ntp[2] = e_cemc_3x3;
0356       ntp[3] = e_hcal_in_3x3;
0357       ntp[4] = e_hcal_out_3x3;
0358       ntp[5] = cemceoverp;
0359       ntp[6] = hcalineovercemce;
0360       ntp[7] = hcaleoverp;
0361       ntp[8] = charge;
0362       ntp[9] = pid;
0363       ntp[10] = quality;
0364       ntp[11] = e_cluster;
0365       ntp[12] = EventNumber;
0366       ntp[13] = z;
0367       ntp[14] = vtxid;
0368       ntp[15] = nmvtx;
0369       ntp[16] = nintt;
0370       ntp[17] = ntpc;
0371       ntp[18] = cemc_prob;
0372       ntp[19] = cemc_ecore;
0373       ntp[20] = cemc_chi2;
0374       if(output_ntuple) { ntpbeforecut -> Fill(ntp); }
0375 
0376     //std::cout << " Pt_lowerlimit " << Pt_lowerlimit << " Pt_higherlimit " << Pt_higherlimit << " HOP_lowerlimit " << HOP_lowerlimit <<std::endl;
0377         //std::cout << " EMOP_lowerlimit " << EMOP_lowerlimit << " EMOP_higherlimit " << EMOP_higherlimit << " HinOEM_higherlimit " << HinOEM_higherlimit <<std::endl;
0378 
0379      ///////////////////////////////////////////////////////////////electrons
0380       if(cemceoverp > EMOP_lowerlimit && cemceoverp < EMOP_higherlimit && quality < Nquality_higherlimit && nmvtx >= Nmvtx_lowerlimit && nintt >= Nintt_lowerlimit && ntpc >= Ntpc_lowerlimit)
0381     {
0382     
0383       ntp[0] = mom;
0384           ntp[1] = pt;
0385           ntp[2] = e_cemc_3x3;
0386           ntp[3] = e_hcal_in_3x3;
0387           ntp[4] = e_hcal_out_3x3;
0388           ntp[5] = charge;
0389           ntp[6] = pid;
0390       ntp[7] = quality;
0391       ntp[8] = e_cluster;
0392       ntp[9] = EventNumber;
0393       ntp[10] = z;
0394       ntp[11] = vtxid;
0395       ntp[12] = cemc_prob;
0396       ntp[13] = cemc_ecore;
0397       if(output_ntuple) { ntpcutEMOP -> Fill(ntp); }
0398 
0399       if(hcalineovercemce < HinOEM_higherlimit)
0400          {
0401 
0402           ntp[0] = mom;
0403               ntp[1] = pt;
0404               ntp[2] = e_cemc_3x3;
0405               ntp[3] = e_hcal_in_3x3;
0406               ntp[4] = e_hcal_out_3x3;
0407               ntp[5] = charge;
0408               ntp[6] = pid;
0409           ntp[7] = quality;
0410           ntp[8] = e_cluster;
0411           ntp[9] = EventNumber;
0412           ntp[10] = z;
0413           ntp[11] = vtxid;
0414           ntp[12] = cemc_prob;
0415           ntp[13] = cemc_ecore;
0416           if(output_ntuple) { ntpcutEMOP_HinOEM -> Fill(ntp); }
0417 
0418           if( pt > Pt_lowerlimit && pt < Pt_higherlimit)
0419             {
0420 
0421               ntp[0] = mom;
0422                   ntp[1] = pt;
0423                   ntp[2] = e_cemc_3x3;
0424                   ntp[3] = e_hcal_in_3x3;
0425                   ntp[4] = e_hcal_out_3x3;
0426                   ntp[5] = charge;
0427                   ntp[6] = pid;
0428               ntp[7] = quality;
0429               ntp[8] = e_cluster;
0430               ntp[9] = EventNumber;
0431               ntp[10] = z;
0432               ntp[11] = vtxid;
0433               ntp[12] = cemc_prob;
0434               ntp[13] = cemc_ecore;
0435               if(output_ntuple) { ntpcutEMOP_HinOEM_Pt -> Fill(ntp); }
0436         
0437               if(Verbosity() > 0) {
0438                 std::cout << " Track " << it->first  << " identified as electron " << "    mom " << mom << " e_cemc_3x3 " << e_cemc_3x3  << " cemceoverp " << cemceoverp << " e_hcal_in_3x3 " << e_hcal_in_3x3 << " e_hcal_out_3x3 " << e_hcal_out_3x3 << std::endl; 
0439                }
0440           
0441               // add to the association map
0442               _track_pid_assoc->addAssoc(TrackPidAssoc::electron, it->second->get_id());
0443            }
0444          }
0445     }
0446      }//end of BDT else
0447       
0448       /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////hadrons
0449      
0450      // if(hcaleoverp > HOP_lowerlimit && quality < 10)
0451       if(hcaleoverp > HOP_lowerlimit && quality < 5)
0452     {
0453       
0454       ntp[0] = mom;
0455           ntp[1] = pt;
0456           ntp[2] = e_cemc_3x3;
0457           ntp[3] = e_hcal_in_3x3;
0458           ntp[4] = e_hcal_out_3x3;
0459           ntp[5] = charge;
0460           ntp[6] = pid;
0461       ntp[7] = quality;
0462           ntp[8] = e_cluster;
0463       ntp[9] = EventNumber;
0464       ntp[10] = z;
0465       ntp[11] = vtxid;
0466       ntp[12] = cemc_prob;
0467       ntp[13] = cemc_ecore;
0468       if(output_ntuple) { ntpcutHOP -> Fill(ntp); }
0469 
0470       if(Verbosity() > 0) {
0471         std::cout << " Track " << it->first  << " identified as hadron " << "    mom " << mom << " e_cemc_3x3 " << e_cemc_3x3  << " hcaleoverp " << hcaleoverp << " e_hcal_in_3x3 " << e_hcal_in_3x3 << " e_hcal_out_3x3 " << e_hcal_out_3x3 << std::endl; 
0472         }
0473 
0474       // add to the association map
0475       _track_pid_assoc->addAssoc(TrackPidAssoc::hadron, it->second->get_id());
0476         }
0477       
0478 
0479     }//end of evet loop.
0480 
0481   delete reader_positive;
0482   delete reader_negative;
0483   
0484   // Read back the association map
0485   if(Verbosity() > 1)
0486     std::cout << "Read back the association map electron entries" << std::endl;
0487   auto electrons = _track_pid_assoc->getTracks(TrackPidAssoc::electron);
0488   for(auto it = electrons.first; it != electrons.second; ++it)
0489     {
0490       SvtxTrack *tr = _track_map->get(it->second);
0491 
0492       double z = tr->get_z();
0493       unsigned int vtxid = tr->get_vertex_id();
0494       double mom = tr->get_p();
0495       double px = tr->get_px();
0496       double py = tr->get_py();
0497       double pt = sqrt(px*px + py*py);
0498       int charge = tr->get_charge();
0499       int pid = it->first;
0500       int tr_id = it->second;
0501       int quality = tr->get_quality();
0502 
0503       double e_cluster = tr->get_cal_cluster_e(SvtxTrack::CAL_LAYER::CEMC);
0504 
0505       double e_cemc_3x3 = tr->get_cal_energy_3x3(SvtxTrack::CAL_LAYER::CEMC);
0506       double e_hcal_in_3x3 = tr->get_cal_energy_3x3(SvtxTrack::CAL_LAYER::HCALIN);
0507       double e_hcal_out_3x3 = tr->get_cal_energy_3x3(SvtxTrack::CAL_LAYER::HCALOUT);
0508 
0509       double EOP = e_cemc_3x3/mom;
0510       double HOM = e_hcal_in_3x3/e_cemc_3x3;
0511 
0512  
0513       unsigned int cemc_clusid = tr->get_cal_cluster_id(SvtxTrack::CAL_LAYER::CEMC);
0514       double cemc_prob = 0.;
0515       double cemc_ecore = 0.;
0516       double cemc_chi2 = 99999.;
0517       if(cemc_clusid<99999) {
0518         RawCluster* cemc_cluster = cemc_cluster_container->getCluster(cemc_clusid);
0519         cemc_prob = cemc_cluster->get_prob();
0520         cemc_ecore = cemc_cluster->get_ecore();
0521         cemc_chi2 = cemc_cluster->get_chi2();
0522       }
0523      
0524       ntp[0] = mom;
0525       ntp[1] = pt;
0526       ntp[2] = e_cemc_3x3;
0527       ntp[3] = e_hcal_in_3x3;
0528       ntp[4] = e_hcal_out_3x3;
0529       ntp[5] = charge;
0530       ntp[6] = pid;
0531       ntp[7] = quality;
0532       ntp[8] = e_cluster;
0533       ntp[9] = EventNumber;
0534       ntp[10] = z;
0535       ntp[11] = vtxid;
0536       ntp[12] = tr_id;
0537       ntp[13] = cemc_prob;
0538       ntp[14] = cemc_ecore;
0539       ntp[15] = EOP;
0540       ntp[16] = HOM;
0541       ntp[17] = cemc_chi2;
0542       //if(output_ntuple) { ntpcutEMOP_HinOEM_Pt_read -> Fill(ntp); }
0543       if(output_ntuple) { ntpcutBDT_read -> Fill(ntp); }
0544       
0545       if(Verbosity() > 1)
0546     std::cout << " pid " << it->first << " track ID " << it->second << " mom " << mom << std::endl;
0547     }
0548 
0549   if(Verbosity() > 1)
0550     std::cout << "Read back the association map hadron entries" << std::endl;
0551   auto hadrons = _track_pid_assoc->getTracks(TrackPidAssoc::hadron);
0552   for(auto it = hadrons.first; it != hadrons.second; ++it)
0553     {
0554       SvtxTrack *tr = _track_map->get(it->second);
0555       double p = tr->get_p();
0556 
0557       if(Verbosity() > 1)
0558     std::cout << " pid " << it->first << " track ID " << it->second << " mom " << p << std::endl;
0559     }
0560 
0561   
0562   return Fun4AllReturnCodes::EVENT_OK;
0563 }
0564 
0565 
0566 int ElectronID::GetNodes(PHCompositeNode* topNode)
0567 {
0568   _track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0569 
0570   // if the TrackPidAssoc node is not present, create it
0571   _track_pid_assoc =  findNode::getClass<TrackPidAssoc>(topNode, "TrackPidAssoc");
0572   if(!_track_pid_assoc)
0573     {
0574       PHNodeIterator iter(topNode);      
0575       PHCompositeNode* dstNode = static_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0576       if (!dstNode)
0577     {
0578       cerr << PHWHERE << "DST Node missing, quit!" << endl;
0579       return Fun4AllReturnCodes::ABORTEVENT;
0580     }
0581 
0582      // Get the SVTX node
0583       PHNodeIterator iter_dst(dstNode);
0584       PHCompositeNode* tb_node =
0585     dynamic_cast<PHCompositeNode*>(iter_dst.findFirst("PHCompositeNode", "SVTX"));
0586       if (!tb_node)
0587     {
0588       cout << PHWHERE << "SVTX node missing, quit!" << endl;
0589       return Fun4AllReturnCodes::ABORTEVENT;
0590     }
0591 
0592       // now add the new node  
0593       _track_pid_assoc = new TrackPidAssoc;
0594       PHIODataNode<PHObject>* assoc_node = new PHIODataNode<PHObject>(_track_pid_assoc, "TrackPidAssoc", "PHObject");
0595       tb_node->addNode(assoc_node);
0596       if (Verbosity() > 0)
0597     cout << PHWHERE << "Svtx/TrackPidAssoc node added" << endl;
0598     }
0599   
0600   return Fun4AllReturnCodes::EVENT_OK;
0601 }
0602 
0603 
0604 PHG4Particle* ElectronID::findMCmatch(SvtxTrack* track, PHG4TruthInfoContainer* truth_container)
0605 {
0606   double px = track->get_px();
0607   double py = track->get_py();
0608   double pz = track->get_pz();
0609   double pt = sqrt(px*px+py*py);
0610   double phi = atan2(py,px);
0611   double eta = asinh(pz/pt);
0612   PHG4TruthInfoContainer::ConstRange range = truth_container->GetPrimaryParticleRange();
0613 //  cout << "number of MC particles = " << truth_container->size() << " " << truth_container->GetNumPrimaryVertexParticles() << endl;
0614 
0615   double thedistance = 9999.;
0616   PHG4Particle* matchedMC = NULL;
0617 
0618     for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0619     {
0620       PHG4Particle* g4particle = iter->second;
0621       int trackid = g4particle->get_track_id();
0622       if(trackid<=truth_container->GetNumPrimaryVertexParticles()-50) continue; // only embedded particles
0623       //int gflavor  = g4particle->get_pid();
0624       double gpx = g4particle->get_px();
0625       double gpy = g4particle->get_py();
0626       double gpz= g4particle->get_pz();
0627       double gpt = sqrt(gpx*gpx+gpy*gpy);
0628       double gphi = atan2(gpy,gpx);
0629       double geta = asinh(gpz/gpt);
0630       //int primid =  g4particle->get_primary_id();
0631       //int parentid = g4particle->get_parent_id();
0632       if(sqrt(pow(gphi-phi,2)+pow(geta-eta,2))<thedistance) {
0633         thedistance = sqrt(pow(gphi-phi,2)+pow(geta-eta,2));
0634         matchedMC = g4particle;
0635       }     
0636     }
0637 
0638   return matchedMC;
0639 }
0640 
0641 
0642 int ElectronID::End(PHCompositeNode * /*topNode*/)
0643 {
0644 if(output_ntuple) {
0645   OutputNtupleFile -> cd();
0646   OutputNtupleFile -> Write();
0647   OutputNtupleFile -> Close();
0648 }
0649 
0650   cout << "************END************" << endl;
0651   return Fun4AllReturnCodes::EVENT_OK;
0652 }
0653
[ 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!