track_EEC/src/trackEEC.cc

0001 #include "trackEEC.h"
0002
0003 #include <fastjet/ClusterSequence.hh>
0004
0005 #include "TMath.h"
0006
0007 #include <trackbase_historic/SvtxTrack.h>
0008 #include <trackbase_historic/SvtxTrackMap.h>
0009 #include <trackbase_historic/TrackSeed.h>
0010
0011
0012 //____________________________________________________________________________..
0013 trackEEC::trackEEC(const std::string &name)
0014   : SubsysReco(name)
0015 {
0016 }
0017
0018 //____________________________________________________________________________..
0019 int trackEEC::Init(PHCompositeNode* /*topNode*/)
0020 {
0021
0022     outfile = new TFile(m_outfileName.c_str(), "RECREATE");
0023
0024     const int nBins = 40;
0025     double bins[nBins+1];
0026     double min_dR = 1e-3;
0027     double max_dR = 0.8;
0028     double dR_shift = 0.0;
0029     double dRWidth = (log(max_dR+dR_shift) - log(min_dR + dR_shift))/nBins;
0030     for(int i=0; i<=nBins; i++){
0031         bins[i] = (min_dR + dR_shift)*exp(i*dRWidth) - dR_shift;
0032     }
0033
0034     h_nEvents = new TH1D("h_nEvents","Number of Events",3,-0.5,2.5);
0035     h_crossing = new TH1D("h_crossing",";crossing",651,-150.5,500.5);
0036     h_nintt = new TH1D("h_nintt",";nintt",6,-0.5,5.5);
0037     h_chi2_ndf = new TH1D("h_chi2_ndf",";#chi^{2}/ndf",3000,0,300);
0038     h_nClus = new TH1D("h_nClus",";Number of Clusters on Track",81,-0.5,80.5);
0039     h_trackSpec = new TH1D("h_trackSpec","Track Spectrum;p_{T} [GeV]",200,0,100);
0040     h_track_etaPhi = new TH2D("h_track_etaPhi","Tracks;#phi;#eta",100,-TMath::Pi(),TMath::Pi(),50,-1.1,1.1);
0041     h_jetCons = new TH1D("h_jetCons",";Number of Jet Constituents",31,-0.5,30.5);
0042     h_jetSpec = new TH1D("h_jetSpec","Jet Spectrum;p_{T} [GeV]",100,0,100);
0043     h_jetConsSpec = new TH2D("h_jetConsSpec",";p_{T} [GeV];Number of Jet Constituents",100,0,100,31,-0.5,30.5);
0044     h_jet_etaPhi = new TH2D("h_jet_etaPhi","Jets;#phi;#eta",100,-TMath::Pi(),TMath::Pi(),50,-1.1,1.1);
0045     h_jetMass = new TH1D("h_jetMass","Jet Mass;M [GeV]",1000,0,100);
0046     h_inJetAll = new TH1D("hEEC_inJetAll","In-jet EEC all jet p_{T};#Delta R",nBins,&bins[0]);
0047     h_inJetAll_Q2 = new TH1D("hEEC_inJetAll_Q2","In-jet EEC all jet p_{T} with Q^{2};#Delta R",nBins,&bins[0]);
0048     h_wEEC = new TH1D("h_wEEC","wEEC, no track jet requirement;#Delta#phi",22,0.0,TMath::Pi());
0049     h_wEEC_jet = new TH1D("h_wEEC_jet","wEEC at least one jet;#Delta#phi",22,0.0,TMath::Pi());
0050     h_wEEC_jet_Q2 = new TH1D("h_wEEC_jet_Q2","wEEC at least one jet with Q^{2};#Delta#phi",22,0.0,TMath::Pi());
0051     h_wEEC_dijet = new TH1D("h_wEEC_dijet","wEEC dijet All Jet p_{T};#Delta#phi",22,0.0,TMath::Pi());
0052     h_wEEC_dijet_Q2 = new TH1D("h_wEEC_dijet_Q2","wEEC dijet with Q^{2} All Jet p_{T};#Delta#phi",22,0.0,TMath::Pi());
0053
0054
0055     for(int i=0; i<5; i++)
0056     {
0057         h_inJet[i] = new TH1D(std::format("hEEC_inJet_{:.0f}_{:.0f}",jet_pTs[i],jet_pTs[i+1]).c_str(),std::format("In-jet EEC {:.0f} < p_{{T}}^{{jet}} < {:.0f};#Delta R",jet_pTs[i],jet_pTs[i+1]).c_str(),nBins,&bins[0]);
0058         h_inJet_Q2[i] = new TH1D(std::format("hEEC_inJetAll_Q2_{:.0f}_{:.0f}",jet_pTs[i],jet_pTs[i+1]).c_str(),std::format("In-jet EEC with Q^2 {:.0f} < p_{{T}}^{{jet}} < {:.0f};#Delta R",jet_pTs[i],jet_pTs[i+1]).c_str(),nBins,&bins[0]);
0059     }
0060
0061     return Fun4AllReturnCodes::EVENT_OK;
0062 }
0063
0064 int trackEEC::process_event(PHCompositeNode *topNode)
0065 {
0066
0067
0068     auto *trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0069     if(!trackmap)
0070     {
0071         return Fun4AllReturnCodes::ABORTEVENT;
0072     }
0073
0074     h_nEvents->Fill(0);
0075
0076     std::vector<fastjet::PseudoJet> parts;
0077
0078     for (const auto& [key, track] : *trackmap)
0079     {
0080         if (!track)
0081         {
0082             continue;
0083         }
0084
0085         //track->identify(std::cout);
0086
0087         h_crossing->Fill(track->get_crossing());
0088         if(track->get_crossing() != 0)
0089         {
0090             continue;
0091         }
0092
0093         float px = track->get_px();
0094         float py = track->get_py();
0095         float pz = track->get_pz();
0096
0097         float pt = sqrt(px*px + py*py);
0098         if(pt < 0.5) continue;
0099         //if(pt > 5) continue;
0100         float eta = std::atanh(pz / std::sqrt(pt*pt + pz*pz));
0101         float phi = std::atan2(py, px);
0102
0103         int n_intt = 0;
0104
0105         std::vector<TrkrDefs::cluskey> clus_keys;
0106         //std::copy(track->begin_cluster_keys(), track->end_cluster_keys(), std::back_inserter(clus_keys));
0107
0108         int nClus = 0;
0109         TrackSeed *sil_seed = track->get_silicon_seed();
0110         if(sil_seed)
0111         {
0112             for(auto iter = sil_seed->begin_cluster_keys(); iter != sil_seed->end_cluster_keys(); ++iter)
0113             {
0114                 TrkrDefs::cluskey ckey = *iter;
0115                 TrkrDefs::TrkrId id = static_cast<TrkrDefs::TrkrId>(TrkrDefs::getTrkrId(ckey));
0116                 //std::cout << "working on sil seed " << nClus << " out of " << sil_seed->size_cluster_keys() << " cluskey: " << ckey << " with TrkrId: " << id << " with name: " << TrkrDefs::TrkrNames.at(id)<< std::endl;
0117                 if(id == TrkrDefs::inttId) n_intt++;
0118                 nClus++;
0119             }
0120         }
0121         int nTPC = 0;
0122         TrackSeed *tpc_seed = track->get_tpc_seed();
0123         if(tpc_seed)
0124         {
0125             for(auto iter = tpc_seed->begin_cluster_keys(); iter != tpc_seed->end_cluster_keys(); ++iter)
0126             {
0127                 TrkrDefs::cluskey ckey = *iter;
0128                 TrkrDefs::TrkrId id = static_cast<TrkrDefs::TrkrId>(TrkrDefs::getTrkrId(ckey));
0129                 //std::cout << "working on tpc seed " << nTPC << " out of " << tpc_seed->size_cluster_keys() << " cluskey: " << ckey << " with TrkrId: " << id << " with name: " << TrkrDefs::TrkrNames.at(id)<< std::endl;
0130                 if(id == TrkrDefs::inttId) n_intt++;
0131                 nClus++;
0132                 nTPC++;
0133             }
0134         }
0135
0136         h_nintt->Fill(n_intt);
0137         h_nClus->Fill(nClus);
0138
0139         if(n_intt <= 0) continue;
0140
0141         h_chi2_ndf->Fill(1.0 * track->get_chisq() / track->get_ndf());
0142
0143         if(1.0 * track->get_chisq() / track->get_ndf() > 20) continue;
0144
0145         h_trackSpec->Fill(pt);
0146         h_track_etaPhi->Fill(phi,eta);
0147
0148         float p = sqrt(px*px + py*py + pz*pz);
0149         float E = sqrt(p*p + m_pi*m_pi);
0150
0151         fastjet::PseudoJet tmp_pj(px, py, pz, E);
0152         parts.push_back(tmp_pj);
0153     }
0154
0155     fastjet::ClusterSequence cs(parts, jetDef);
0156     std::vector<fastjet::PseudoJet> jets = sorted_by_pt( jetSel( cs.inclusive_jets() ) );
0157
0158     bool goodDijet = false;
0159     if((int)jets.size() >= 2)
0160     {
0161         double dPhiJets = std::abs(jets[0].phi_std() - jets[1].phi_std());
0162         if(dPhiJets > TMath::Pi()) dPhiJets = 2*TMath::Pi() - dPhiJets;
0163
0164         if(dPhiJets > 3.0*TMath::Pi()/4.0) goodDijet = true;
0165     }
0166
0167     if((int)jets.size() >= 1) h_nEvents->Fill(1);
0168     if(goodDijet) h_nEvents->Fill(2);
0169
0170     for(int i=0; i<(int)parts.size(); i++)
0171     {
0172        for(int j=i+1; j<(int)parts.size(); j++)
0173        {
0174             double dPhi = std::abs(parts[i].phi_std() - parts[j].phi_std());
0175             if(dPhi > TMath::Pi()) dPhi = 2*TMath::Pi() - dPhi;
0176
0177             h_wEEC->Fill(dPhi, parts[i].perp() * parts[j].perp());
0178             if(jets.size() >= 1)
0179             {
0180                 h_wEEC_jet->Fill(dPhi, parts[i].perp() * parts[j].perp());
0181                 h_wEEC_jet_Q2->Fill(dPhi, parts[i].perp() * parts[j].perp()/pow(jets[0].perp(),2));
0182             }
0183             if(goodDijet)
0184             {
0185                 h_wEEC_dijet->Fill(dPhi, parts[i].perp() * parts[j].perp());
0186                 h_wEEC_dijet_Q2->Fill(dPhi, parts[i].perp() * parts[j].perp() / pow(0.5*(jets[0].perp() + jets[1].perp()),2));
0187             }
0188        }
0189     }
0190
0191     std::vector<fastjet::PseudoJet> EECCons;
0192     for(auto jet : jets)
0193     {
0194
0195         //if(jet.m() < 0.5) continue;
0196
0197         h_jetCons->Fill((int) jet.constituents().size());
0198
0199         h_jetSpec->Fill(jet.perp());
0200         h_jetConsSpec->Fill(jet.perp(), (int) jet.constituents().size());
0201         h_jet_etaPhi->Fill(jet.phi_std(), jet.pseudorapidity());
0202         h_jetMass->Fill(jet.m());
0203
0204         int jet_pT_index = -1;
0205         for(int i=0; i<5; i++)
0206         {
0207             if(jet.perp() >= jet_pTs[i] && jet.perp() < jet_pTs[i+1])
0208             {
0209                 jet_pT_index = i;
0210                 break;
0211             }
0212         }
0213
0214         EECCons.clear();
0215         for(auto part : parts)
0216         {
0217             double dPhi = std::abs(jet.phi_std() - part.phi_std());
0218             if(dPhi > TMath::Pi()) dPhi = 2*TMath::Pi() - dPhi;
0219
0220             double dEta = jet.pseudorapidity() - part.pseudorapidity();
0221
0222             double dR = sqrt(dPhi*dPhi + dEta*dEta);
0223
0224             if(dR < 0.4) EECCons.push_back(part);
0225         }
0226
0227         for(int i=0; i<(int)EECCons.size(); i++)
0228         {
0229             for(int j=i+1; j<(int)EECCons.size(); j++)
0230             {
0231                 double dPhi = std::abs(EECCons[i].phi_std() - EECCons[j].phi_std());
0232                 if(dPhi > TMath::Pi()) dPhi = 2*TMath::Pi() - dPhi;
0233                 double dEta = EECCons[i].pseudorapidity() - EECCons[j].pseudorapidity();
0234                 double dR = sqrt(dPhi*dPhi + dEta*dEta);
0235
0236                 h_inJetAll->Fill(dR, EECCons[i].perp() * EECCons[j].perp());
0237                 h_inJetAll_Q2->Fill(dR, EECCons[i].perp() * EECCons[j].perp() / pow(jet.perp(),2));
0238
0239                 if(jet_pT_index >= 0)
0240                 {
0241                     h_inJet[jet_pT_index]->Fill(dR, EECCons[i].perp() * EECCons[j].perp());
0242                     h_inJet_Q2[jet_pT_index]->Fill(dR, EECCons[i].perp() * EECCons[j].perp() / pow(jet.perp(),2));
0243                 }
0244             }
0245         }
0246     }
0247
0248     return Fun4AllReturnCodes::EVENT_OK;
0249 }
0250
0251 int trackEEC::End(PHCompositeNode* /*topNode */)
0252 {
0253
0254     outfile->cd();
0255
0256     h_nEvents->Write();
0257     h_crossing->Write();
0258     h_nintt->Write();
0259     h_chi2_ndf->Write();
0260     h_nClus->Write();
0261     h_trackSpec->Write();
0262     h_track_etaPhi->Write();
0263     h_jetCons->Write();
0264     h_jetSpec->Write();
0265     h_jetConsSpec->Write();
0266     h_jet_etaPhi->Write();
0267     h_jetMass->Write();
0268     h_inJetAll->Write();
0269     h_inJetAll_Q2->Write();
0270     for(int i=0; i<5; i++)
0271     {
0272         h_inJet[i]->Write();
0273         h_inJet_Q2[i]->Write();
0274     }
0275     h_wEEC->Write();
0276     h_wEEC_jet->Write();
0277     h_wEEC_jet_Q2->Write();
0278     h_wEEC_dijet->Write();
0279     h_wEEC_dijet_Q2->Write();
0280
0281     outfile->Close();
0282
0283     return Fun4AllReturnCodes::EVENT_OK;
0284 }