FilterEvents/src/FilterEvents.cc

0001
0002 #include "FilterEvents.h"
0003
0004 #include <cstdlib>
0005 #include <cstdio>
0006 #include <iostream>
0007 #include <iomanip>
0008 #include <fstream>
0009
0010 #include <TFile.h>
0011
0012 #include <trackbase_historic/SvtxTrack.h>
0013 #include <trackbase_historic/SvtxTrackMap.h>
0014 #include <trackbase_historic/SvtxTrackMap_v2.h>
0015 #include <trackbase_historic/SvtxVertex.h>
0016 #include <trackbase_historic/SvtxVertexMap_v1.h>
0017 #include <trackbase_historic/TrackSeed.h>
0018 #include <trackbase_historic/TrackSeed_v1.h>
0019 #include <trackbase_historic/SvtxTrackSeed_v1.h>
0020 #include <trackbase_historic/TrackSeedContainer.h>
0021 #include <trackbase/TrkrDefs.h>
0022 #include <trackbase/TrkrCluster.h>
0023 #include <trackbase/TrkrClusterv4.h>
0024 #include <trackbase/TrkrClusterContainer.h>
0025 #include <trackbase/TrkrClusterContainerv4.h>
0026 #include <g4vertex/GlobalVertexMap.h>
0027 #include <g4vertex/GlobalVertex.h>
0028 #include <calobase/RawClusterContainer.h>
0029 #include <calobase/RawCluster.h>
0030 #include <calobase/RawClusterv1.h>
0031
0032 #include <KFParticle.h>
0033
0034 #include <phool/getClass.h>
0035 #include <phool/recoConsts.h>
0036 #include <phool/PHCompositeNode.h>
0037 #include <phool/PHIODataNode.h>
0038 #include <phool/PHNodeIterator.h>
0039 #include <phool/PHRandomSeed.h>
0040 #include <fun4all/Fun4AllReturnCodes.h>
0041
0042 typedef PHIODataNode<PHObject> PHObjectNode_t;
0043
0044 using namespace std;
0045
0046 //==============================================================
0047
0048 FilterEvents::FilterEvents(const std::string &name) : SubsysReco(name)
0049 {
0050   outnodename_trackmap = "SvtxTrackMap_ee";
0051   outnodename_cemc_clusters = "CLUSTER_CEMC_ee";
0052   EventNumber=0;
0053   goodEventNumber=0;
0054
0055   pt_cut = 2;
0056   dca_cut = 0;
0057   chi2ndof_cut = 999.9;
0058   CEMC_use = true;
0059 }
0060
0061 //-------------------------------------------------------------------------------
0062
0063 int FilterEvents::Init(PHCompositeNode *topNode)
0064 {
0065
0066 cout << "FilterEvents::Init started..." << endl;
0067
0068   PHNodeIterator iter(topNode);
0069   //PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0070   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0071   if (!dstNode) { cerr << PHWHERE << " ERROR: DST node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0072
0073 //  PHNodeIterator dstiter(dstNode);
0074 //  PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "SVTX"));
0075 //  if (!svtxNode) { cerr << PHWHERE << " ERROR: SVTX node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0076 //
0077 //  PHCompositeNode *cemcNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "CEMC"));
0078 //  if (!cemcNode) { cerr << PHWHERE << " ERROR: CEMC node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0079
0080   SvtxTrackMap_v2* trackmap = new SvtxTrackMap_v2();
0081   PHCompositeNode *trackmapNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", outnodename_trackmap));
0082   if (!trackmapNode)
0083   {
0084     PHObjectNode_t *trackmapNode = new PHIODataNode<PHObject>(trackmap,outnodename_trackmap.c_str(),"PHObject");
0085       dstNode->addNode(trackmapNode);
0086     //svtxNode->addNode(trackmapNode);
0087     cout << PHWHERE << " INFO: added " << outnodename_trackmap << " node." << endl;
0088   }
0089   else { cout << PHWHERE << " INFO: " << outnodename_trackmap << " node already exists." << endl; }
0090
0091   if(CEMC_use){
0092      _cemc_clusters_ee = new RawClusterContainer();
0093      PHIODataNode<PHObject> *clusterNode = new PHIODataNode<PHObject>(_cemc_clusters_ee, outnodename_cemc_clusters, "PHObject");
0094      dstNode->addNode(clusterNode);
0095   //cemcNode->addNode(clusterNode);
0096   }
0097 //  PHCompositeNode *cemcclusNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", outnodename_cemc_clusters));
0098 //  if (!cemcclusNode)
0099 //  {
0100 //    PHObjectNode_t *cemcclusNode = new PHIODataNode<PHObject>(cemc_clusters,outnodename_cemc_clusters.c_str(),"PHObject");
0101 //    dstNode->addNode(cemcclusNode);
0102 //    cout << PHWHERE << " INFO: added " << outnodename_cemc_clusters << " node." << endl;
0103 //  }
0104 //  else { cout << PHWHERE << " INFO: " << outnodename_cemc_clusters << " node already exists." << endl; }
0105
0106 cout << "FilterEvents::Init ended." << endl;
0107   return Fun4AllReturnCodes::EVENT_OK;
0108 }
0109
0110 //--------------------------------------------------------------------------------
0111
0112 int FilterEvents::InitRun(PHCompositeNode *topNode)
0113 {
0114   return Fun4AllReturnCodes::EVENT_OK;
0115 }
0116
0117 //--------------------------------------------------------------------------------
0118
0119 int FilterEvents::GetNodes(PHCompositeNode *topNode)
0120 {
0121
0122   _topNode = topNode;
0123
0124   //topNode->print();
0125
0126   _trackmap_ee = findNode::getClass<SvtxTrackMap>(topNode, outnodename_trackmap);
0127   if(!_trackmap_ee) { cerr << PHWHERE << "ERROR: Output SvtxTrackMap_ee node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0128 //    else {cout << "Found SvtxTrackMap_ee node." << endl; }
0129
0130   if (CEMC_use){
0131     _cemc_clusters_ee = findNode::getClass<RawClusterContainer>(topNode, outnodename_cemc_clusters);
0132     if(!_cemc_clusters_ee) { cerr << PHWHERE << "ERROR: CLUSTER_CEMC_ee node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0133 //    else {cout << "Found CLUSTER_CEMC_ee node." << endl; }
0134   }
0135
0136   _trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0137   if(!_trackmap) { cerr << PHWHERE << "ERROR: SvtxTrackMap node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0138
0139
0140   _vtxmap = findNode::getClass<SvtxVertexMap_v1>(topNode, "SvtxVertexMap");
0141   if(!_vtxmap) { cout << "ERROR: SvtxVertexMap node not found!" << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0142
0143   _trackseedcontainer_svtx = findNode::getClass<TrackSeedContainer>(topNode, "SvtxTrackSeedContainer");
0144   if(!_trackseedcontainer_svtx) { cerr << PHWHERE << "ERROR: SvtxTrackSeedContainer node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0145
0146   _trackseedcontainer_silicon = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
0147   if(!_trackseedcontainer_silicon) { cerr << PHWHERE << "ERROR: SiliconTrackSeedContainer node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0148
0149   _trackseedcontainer_tpc = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
0150   if(!_trackseedcontainer_tpc) { cerr << PHWHERE << "ERROR: TpcTrackSeedContainer node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0151
0152   _trkrclusters  = findNode::getClass<TrkrClusterContainerv4>(topNode, "TRKR_CLUSTER");
0153   if(!_trkrclusters) { cerr << PHWHERE << "ERROR: TRKR_CLUSTER node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0154
0155   if (CEMC_use){
0156     _cemc_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
0157     if(!_cemc_clusters) { cerr << PHWHERE << "ERROR: CLUSTER_CEMC node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0158   }
0159   return Fun4AllReturnCodes::EVENT_OK;
0160 }
0161
0162 //--------------------------------------------------------------------
0163
0164 TVector3 FilterEvents::GetProjectionCEMC(SvtxTrack* track) {
0165
0166   TVector3 projection; // 0,0,0
0167
0168   vector<double> proj;
0169   for (SvtxTrack::StateIter stateiter = track->begin_states(); stateiter != track->end_states(); ++stateiter)
0170   {
0171     SvtxTrackState *trackstate = stateiter->second;
0172     if(trackstate) { proj.push_back(trackstate->get_pathlength()); }
0173   }
0174   double pathlength = proj[proj.size()-3]; // CEMC is next next to last, usually 93.5
0175
0176   SvtxTrackState* trackstate = track->get_state(pathlength); // at CEMC inner face
0177   if(trackstate) {
0178     double track_x = trackstate->get_x();
0179     double track_y = trackstate->get_y();
0180     double track_z = trackstate->get_z();
0181     projection.SetX(track_x);
0182     projection.SetY(track_y);
0183     projection.SetZ(track_z);
0184   }
0185
0186   return projection;
0187 }
0188
0189 //--------------------------------------------------------------------
0190
0191 RawCluster* FilterEvents::MatchClusterCEMC(SvtxTrack* track, RawClusterContainer* cemc_clusters, double &dphi, double &deta, double Zvtx) {
0192
0193   RawCluster* returnCluster = NULL;
0194   double track_eta = 99999.;
0195   double track_phi = 99999.;
0196   dphi = 99999.;
0197   deta = 99999.;
0198
0199   vector<double> proj;
0200   for (SvtxTrack::StateIter stateiter = track->begin_states(); stateiter != track->end_states(); ++stateiter)
0201   {
0202     SvtxTrackState *trackstate = stateiter->second;
0203     if(trackstate) { proj.push_back(trackstate->get_pathlength()); }
0204   }
0205   double pathlength = proj[proj.size()-3]; // CEMC is next next to last
0206
0207   SvtxTrackState* trackstate = track->get_state(pathlength); // at CEMC inner face
0208   if(trackstate) {
0209     double track_x = trackstate->get_x();
0210     double track_y = trackstate->get_y();
0211     double track_z = trackstate->get_z() - Zvtx;
0212     double track_r = sqrt(track_x*track_x+track_y*track_y);
0213     track_eta = asinh( track_z / track_r );
0214     track_phi = atan2( track_y, track_x );
0215   } else { return returnCluster; }
0216
0217   if(track_eta == 99999. || track_phi == 99999.) { return returnCluster; }
0218   double dist = 99999.;
0219
0220     RawClusterContainer::Range begin_end = cemc_clusters->getClusters();
0221     RawClusterContainer::Iterator iter;
0222     for (iter = begin_end.first; iter != begin_end.second; ++iter)
0223     {
0224       RawCluster* cluster = iter->second;
0225       if(!cluster) continue;
0226       else {
0227         double cemc_ecore = cluster->get_ecore();
0228         if(cemc_ecore<1.0) continue;
0229         double cemc_x = cluster->get_x();
0230         double cemc_y = cluster->get_y();
0231         double cemc_z = cluster->get_z() - Zvtx;
0232         double cemc_r = cluster->get_r();
0233         double cemc_eta = asinh( cemc_z / cemc_r );
0234         double cemc_phi = atan2( cemc_y, cemc_x );
0235         double tmpdist = sqrt(pow((cemc_eta-track_eta),2)+pow((cemc_phi-track_phi),2));
0236         if(tmpdist<dist) {
0237           dist = tmpdist; returnCluster = cluster; dphi = fabs(cemc_phi-track_phi); deta = fabs(cemc_eta-track_eta);
0238         }
0239       }
0240     }
0241
0242   return returnCluster;
0243 }
0244
0245 //--------------------------------------------------------------------------------
0246
0247 int FilterEvents::process_event(PHCompositeNode *topNode) {
0248
0249   EventNumber++;
0250   int ngood = 0;
0251   bool goodevent = false;
0252
0253 cout << "getting nodes..." << endl;
0254   GetNodes(topNode);
0255
0256   GlobalVertexMap* _global_vtxmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0257   if(!_global_vtxmap) { cerr << PHWHERE << "ERROR: GlobalVertexMap node not found." << endl; return Fun4AllReturnCodes::ABORTEVENT; }
0258
0259   double Zvtx = 0.;
0260   for (GlobalVertexMap::Iter iter = _global_vtxmap->begin(); iter != _global_vtxmap->end(); ++iter)
0261   {
0262     GlobalVertex *vtx = iter->second;
0263     if(vtx->get_id()==1) { Zvtx = vtx->get_z(); } // BBC vertex
0264   }
0265   cout << "Global BBC vertex Z = " << Zvtx << endl;
0266
0267   std::vector<RawCluster *> goodclusters;
0268   //std::vector<TrkrCluster*> vclussilicon;
0269   //std::vector<TrkrCluster*> vclustpc;
0270   std::map<TrkrDefs::cluskey,TrkrCluster*> m_clusters;
0271   //std::map<TrkrDefs::cluskey> vcluskeytpc;
0272   std::vector<SvtxTrackSeed_v1*> v_svtx_trackseed;
0273   std::vector<TrackSeed*> v_silicon_trackseed;
0274   std::vector<TrackSeed*> v_tpc_trackseed;
0275
0276   cout << "Total number of tracks = " << _trackmap->size() << endl;
0277   if (CEMC_use) cout << "Total number of CEMC clusters = " << _cemc_clusters->size() << endl;
0278   cout << "Total number of svtx track seeds: " << _trackseedcontainer_svtx->size() << endl;
0279   cout << "Total number of silicon track seeds: " << _trackseedcontainer_silicon->size() << endl;
0280   cout << "Total number of tpc track seeds: " << _trackseedcontainer_tpc->size() << endl;
0281   cout << "Total number of TRKR clusters: " << _trkrclusters->size() << endl;
0282
0283   // Start loop over tracks;
0284
0285   for (SvtxTrackMap::Iter iter = _trackmap->begin(); iter != _trackmap->end(); ++iter)
0286   {
0287     SvtxTrack *track = iter->second;
0288
0289     double px = track->get_px();
0290     double py = track->get_py();
0291     double pz = track->get_pz();
0292     double pt = sqrt(px * px + py * py);
0293     double chi2ndf = track->get_chisq()/track->get_ndf();
0294     double dca;
0295     double dcaxy;
0296     double dcaz, dcaxy_sig, dcaz_sig;
0297
0298     get_dca(track, _vtxmap, dca, dcaxy);
0299     get_dca_SvtxEval(track, _vtxmap, dcaxy, dcaz, dcaxy_sig, dcaz_sig);
0300       if(pt < pt_cut) continue;
0301       if(chi2ndf > chi2ndof_cut) continue;
0302       if(abs(dca) < dca_cut) continue;
0303
0304     double mom = sqrt(px * px + py * py + pz * pz);
0305
0306     double cemc_dphi = 99999.;
0307     double cemc_deta = 99999.;
0308     if (CEMC_use)
0309     {
0310       RawCluster* clus;
0311       clus = MatchClusterCEMC(track, _cemc_clusters, cemc_dphi, cemc_deta, Zvtx);
0312       if(!clus && CEMC_use) continue;
0313
0314       double cemc_ecore = 0.;
0315       cemc_ecore = clus->get_ecore();
0316       if(cemc_ecore/mom < 0.7) continue;
0317     }
0318
0319     ngood++;
0320     goodevent = true;
0321     SvtxTrack* tmp =_trackmap_ee->insert(track);
0322     if(!tmp) cout << "ERROR: Failed to insert a track." << endl;
0323
0324     cout << " Track: " << pt << endl;
0325
0326     TrackSeed* trackseed_silicon = track->get_silicon_seed();
0327     double trackseed_silicon_pt = trackseed_silicon->get_pt();
0328     cout << "   Silicon seed: " << trackseed_silicon << " " << trackseed_silicon_pt << " " << trackseed_silicon->size_cluster_keys() << endl;
0329     TrackSeed_v1* tmpseedsilicon = (TrackSeed_v1*)trackseed_silicon->CloneMe();
0330     v_silicon_trackseed.push_back(tmpseedsilicon);
0331
0332     TrackSeed* trackseed_tpc = track->get_tpc_seed();
0333     double trackseed_tpc_pt = trackseed_tpc->get_pt();
0334     cout << "   TPC seed: " << trackseed_tpc << " "  << trackseed_tpc_pt << " " << trackseed_tpc->size_cluster_keys() << endl;
0335     TrackSeed_v1* tmpseedtpc = (TrackSeed_v1*)trackseed_tpc->CloneMe();
0336     v_tpc_trackseed.push_back(tmpseedtpc);
0337
0338 // Find SVTX seed for this track
0339     for(TrackSeedContainer::ConstIter seediter = _trackseedcontainer_svtx->begin(); seediter != _trackseedcontainer_svtx->end(); ++seediter)
0340     {
0341       SvtxTrackSeed_v1* seed = (SvtxTrackSeed_v1*)*seediter;
0342       //TrackSeed* seed = *seediter;
0343       bool foundsiliconseed = false;
0344       bool foundtpcseed = false;
0345       unsigned int siliconid = seed->get_silicon_seed_index();
0346       unsigned int tpcid = seed->get_tpc_seed_index();
0347       //cout << "   SVTX seed: " << seed << " " << siliconid << " " << tpcid << endl;
0348       TrackSeed* tmptpcseed  = _trackseedcontainer_tpc->get(tpcid);
0349       double tmptpcseed_pt = tmptpcseed->get_pt();
0350       TrackSeed* tmpsiliconseed  = _trackseedcontainer_silicon->get(siliconid);
0351       double tmpsiliconseed_pt = tmpsiliconseed->get_pt();
0352       if(tmpsiliconseed_pt == trackseed_silicon_pt) {foundsiliconseed = true; cout << "   Found silicon seed " << tmpsiliconseed << endl;}
0353       if(tmptpcseed_pt == trackseed_tpc_pt) {foundtpcseed = true; cout << "   Found tpc seed " << tmptpcseed << endl;}
0354       if(foundsiliconseed && foundtpcseed) {
0355         cout << "     THIS IS THE ONE: " << seed << " " << siliconid << " " << tpcid << endl;
0356         SvtxTrackSeed_v1* tmpseed = (SvtxTrackSeed_v1*)seed->CloneMe();
0357         v_svtx_trackseed.push_back(tmpseed);
0358       }
0359     }
0360
0361 // Find all TRKR clusters for this track
0362     for(auto clusiter = trackseed_silicon->begin_cluster_keys(); clusiter != trackseed_silicon->end_cluster_keys(); ++clusiter)
0363     {
0364       //auto key = *clusiter;
0365       TrkrDefs::cluskey key = *clusiter;
0366       TrkrCluster* clus = _trkrclusters->findCluster(key);
0367       bool isinserted = false;
0368       for(std::map<TrkrDefs::cluskey,TrkrCluster*>::iterator it = m_clusters.begin(); it != m_clusters.end(); it++) {if(clus==it->second) {isinserted=true; break;}}
0369       if(!isinserted) {TrkrCluster* newclus = (TrkrClusterv4*)clus->CloneMe(); m_clusters.insert(std::make_pair(key,newclus));}
0370     }
0371     for(auto clusiter = trackseed_tpc->begin_cluster_keys(); clusiter != trackseed_tpc->end_cluster_keys(); ++clusiter)
0372     {
0373       //auto key = *clusiter;
0374       TrkrDefs::cluskey key = *clusiter;
0375       TrkrCluster* clus = _trkrclusters->findCluster(key);
0376       bool isinserted = false;
0377       for(std::map<TrkrDefs::cluskey,TrkrCluster*>::iterator it = m_clusters.begin(); it != m_clusters.end(); it++) {if(clus==it->second) {isinserted=true; break;}}
0378       if(!isinserted) {TrkrCluster* newclus = (TrkrClusterv4*)clus->CloneMe(); m_clusters.insert(std::make_pair(key,newclus));}
0379     }
0380
0381     TVector3 proj = GetProjectionCEMC(track);
0382     double track_x = proj(0);
0383     double track_y = proj(1);
0384     double track_z = proj(2) - Zvtx;
0385     double track_r = sqrt(track_x*track_x+track_y*track_y);
0386     double track_eta = asinh( track_z / track_r );
0387     double track_phi = atan2( track_y, track_x );
0388
0389   if (CEMC_use){
0390 // Find all CEMC clusters around this track projection
0391       RawClusterContainer::Range begin_end = _cemc_clusters->getClusters();
0392       RawClusterContainer::Iterator clusiter;
0393       for (clusiter = begin_end.first; clusiter != begin_end.second; ++clusiter)
0394       {
0395         RawCluster* cluster = clusiter->second;
0396         if(!cluster) { cout << "ERROR: bad cluster pointer = " << cluster << endl; continue; }
0397         else {
0398           double cemc_ecore = cluster->get_ecore();
0399           if(cemc_ecore<1.0) continue;
0400           double cemc_x = cluster->get_x();
0401           double cemc_y = cluster->get_y();
0402           double cemc_z = cluster->get_z() - Zvtx;
0403           double cemc_r = cluster->get_r();
0404           double cemc_eta = asinh( cemc_z / cemc_r );
0405           double cemc_phi = atan2( cemc_y, cemc_x );
0406           double dist = sqrt(pow(cemc_phi-track_phi,2)+pow(cemc_eta-track_eta,2));
0407           if(dist<0.1) {
0408             bool isinserted = false;
0409             for(unsigned int i=0; i<goodclusters.size(); i++) {if(cluster==goodclusters[i]) {isinserted=true; break;}}
0410             if(!isinserted) {RawCluster* newcluster = (RawClusterv1*)cluster->CloneMe(); goodclusters.push_back(newcluster);}
0411           }
0412         }
0413       } // end loop over cemc clusters
0414
0415    }
0416   } // end loop over tracks
0417
0418   cout << "   Number of CEMC clusters for output = " << goodclusters.size() << endl;
0419   cout << "   Number of svtx seeds for output = " << v_svtx_trackseed.size() << endl;
0420   cout << "   Number of tpc seeds for output = " << v_tpc_trackseed.size() << endl;
0421   cout << "   Number of silicon seeds for output = " << v_silicon_trackseed.size() << endl;
0422
0423 // Fill in selected CEMC clusters
0424   if (CEMC_use){
0425     _cemc_clusters->Reset();
0426     for (unsigned int cl = 0; cl < goodclusters.size(); cl++) { _cemc_clusters->AddCluster(goodclusters[cl]); }
0427   }
0428
0429   cout << "   Number of TRKR clusters and keys for output: " << m_clusters.size() << " " << m_clusters.size() << endl;
0430
0431   _trkrclusters->Reset();
0432   for(std::map<TrkrDefs::cluskey,TrkrCluster*>::iterator it = m_clusters.begin(); it != m_clusters.end(); it++)
0433      {
0434        //std::cout << it->first << std::endl;
0435        _trkrclusters->addClusterSpecifyKey(it->first,it->second);
0436      }
0437   cout << "   New TRKR_CLUSTER size = " << _trkrclusters->size() << endl;
0438
0439
0440   _trackseedcontainer_svtx->Reset();
0441   for(unsigned int cl = 0; cl < v_svtx_trackseed.size(); cl++)
0442     { _trackseedcontainer_svtx->insert(v_svtx_trackseed[cl]); }
0443
0444 /*
0445   _trackseedcontainer_silicon->Reset();
0446   for(unsigned int cl = 0; cl < v_silicon_trackseed.size(); cl++)
0447     { _trackseedcontainer_silicon->insert(v_silicon_trackseed[cl]); }
0448   _trackseedcontainer_tpc->Reset();
0449   for(unsigned int cl = 0; cl < v_tpc_trackseed.size(); cl++)
0450     { _trackseedcontainer_tpc->insert(v_tpc_trackseed[cl]); }
0451 */
0452
0453   cout << "      New numbers of seeds = " << _trackseedcontainer_svtx->size() << " " << _trackseedcontainer_silicon->size() << " " << _trackseedcontainer_tpc->size() << endl;
0454   if(ngood>=2) { goodEventNumber++; }
0455
0456   cout << ngood << " " << EventNumber << " " << goodEventNumber << endl;
0457   if(goodevent) {return Fun4AllReturnCodes::EVENT_OK;}  else {return Fun4AllReturnCodes::ABORTEVENT;}
0458 }
0459
0460 //---------------------------------------------------------------------------
0461
0462 int FilterEvents::End(PHCompositeNode *topNode)
0463 {
0464   cout << "Number of scanned events = " << EventNumber << endl;
0465   cout << "Number of good events = " << goodEventNumber << endl;
0466   return Fun4AllReturnCodes::EVENT_OK;
0467 }
0468
0469 void FilterEvents::get_dca(SvtxTrack* track, SvtxVertexMap_v1* vertexmap,
0470                                   double& DCA, double& DCAxy)
0471 {
0472     KFParticle::SetField(-1.4e0);
0473
0474     float trackParameters[6] = {track->get_x(),
0475                                 track->get_y(),
0476                                 track->get_z(),
0477                                 track->get_px(),
0478                                 track->get_py(),
0479                                 track->get_pz()};
0480
0481     float trackCovariance[21];
0482     unsigned int iterate = 0;
0483     for (unsigned int i = 0; i < 6; ++i)
0484     {
0485         for (unsigned int j = 0; j <= i; ++j)
0486         {
0487             trackCovariance[iterate] = track->get_error(i, j);
0488             ++iterate;
0489         }
0490     }
0491
0492     KFParticle particle;
0493     particle.Create(trackParameters, trackCovariance, (Int_t) track->get_charge(), -1);
0494     particle.NDF() = track->get_ndf();
0495     particle.Chi2() = track->get_chisq();
0496     particle.SetId(track->get_id());
0497
0498     auto vtxid = track->get_vertex_id();
0499     auto svtxVertex = vertexmap->get(vtxid);
0500     if(!svtxVertex)
0501     {
0502       return;
0503     }
0504
0505
0506     float vertexParameters[6] = {svtxVertex->get_x(),
0507                                  svtxVertex->get_y(),
0508                                  svtxVertex->get_z(), 0, 0, 0};
0509
0510     float vertexCovariance[21];
0511     iterate = 0;
0512     for (unsigned int i = 0; i < 3; ++i)
0513     {
0514         for (unsigned int j = 0; j <= i; ++j)
0515         {
0516             vertexCovariance[iterate] = svtxVertex->get_error(i, j);
0517             ++iterate;
0518         }
0519     }
0520
0521     KFParticle vertex;
0522     vertex.Create(vertexParameters, vertexCovariance, 0, -1);
0523     vertex.NDF() = svtxVertex->get_ndof();
0524     vertex.Chi2() = svtxVertex->get_chisq();
0525
0526     DCA = particle.GetDistanceFromVertex(vertex);
0527     DCAxy = particle.GetDistanceFromVertexXY(vertex);
0528 }
0529
0530 void FilterEvents::get_dca_SvtxEval(SvtxTrack* track, SvtxVertexMap_v1* vertexmap,
0531                 double& dca3dxy, double& dca3dz,
0532                 double& dca3dxysigma, double& dca3dzsigma)
0533 {
0534   Acts::Vector3 pos(track->get_x(),
0535             track->get_y(),
0536             track->get_z());
0537   Acts::Vector3 mom(track->get_px(),
0538             track->get_py(),
0539             track->get_pz());
0540
0541   auto vtxid = track->get_vertex_id();
0542   auto svtxVertex = vertexmap->get(vtxid);
0543   if(!svtxVertex)
0544     { return; }
0545   Acts::Vector3 vertex(svtxVertex->get_x(),
0546                svtxVertex->get_y(),
0547                svtxVertex->get_z());
0548
0549   pos -= vertex;
0550
0551   Acts::ActsSymMatrix<3> posCov;
0552   for(int i = 0; i < 3; ++i)
0553     {
0554       for(int j = 0; j < 3; ++j)
0555     {
0556       posCov(i, j) = track->get_error(i, j);
0557     }
0558     }
0559
0560   Acts::Vector3 r = mom.cross(Acts::Vector3(0.,0.,1.));
0561   float phi = atan2(r(1), r(0));
0562
0563   Acts::RotationMatrix3 rot;
0564   Acts::RotationMatrix3 rot_T;
0565   rot(0,0) = cos(phi);
0566   rot(0,1) = -sin(phi);
0567   rot(0,2) = 0;
0568   rot(1,0) = sin(phi);
0569   rot(1,1) = cos(phi);
0570   rot(1,2) = 0;
0571   rot(2,0) = 0;
0572   rot(2,1) = 0;
0573   rot(2,2) = 1;
0574
0575   rot_T = rot.transpose();
0576
0577   Acts::Vector3 pos_R = rot * pos;
0578   Acts::ActsSymMatrix<3> rotCov = rot * posCov * rot_T;
0579
0580   dca3dxy = pos_R(0);
0581   dca3dz = pos_R(2);
0582   dca3dxysigma = sqrt(rotCov(0,0));
0583   dca3dzsigma = sqrt(rotCov(2,2));
0584
0585 }