packages/trackreco/PHCosmicSeeder.cc

0001
0002 #include "PHCosmicSeeder.h"
0003
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <fun4all/SubsysReco.h>
0006 #include <phool/PHCompositeNode.h>
0007 #include <phool/PHIODataNode.h>
0008 #include <phool/PHNodeIterator.h>
0009 #include <phool/PHObject.h>
0010 #include <phool/getClass.h>
0011 #include <phool/phool.h>
0012
0013 #include <trackbase/TrkrClusterContainer.h>
0014 #include <trackbase/TrkrDefs.h>
0015 #include <trackbase_historic/TrackSeedContainer.h>
0016 #include <trackbase_historic/TrackSeedContainer_v1.h>
0017 #include <trackbase_historic/TrackSeed_v2.h>
0018 #include <trackbase/TrkrCluster.h>
0019
0020 #include <TFile.h>
0021 #include <TNtuple.h>
0022
0023 #include <algorithm>
0024 #include <cmath>
0025 #include <iostream>
0026 #include <iterator>
0027 #include <memory>
0028 #include <ostream>
0029 #include <set>
0030 #include <stdexcept>
0031 #include <string>
0032 #include <utility>
0033 #include <vector>
0034 namespace
0035 {
0036   template <class T>
0037   inline constexpr T square(const T& x)
0038   {
0039     return x * x;
0040   }
0041   template <class T>
0042   inline constexpr T r(const T& x, const T& y)
0043   {
0044     double sign = 1;
0045     if (y < 0)
0046     {
0047       sign = -1;
0048     }
0049     return std::sqrt(square(x) + square(y)) * sign;
0050   }
0051 }  // namespace
0052 //____________________________________________________________________________..
0053 PHCosmicSeeder::PHCosmicSeeder(const std::string& name)
0054   : SubsysReco(name)
0055 {
0056 }
0057
0058 //____________________________________________________________________________..
0059 int PHCosmicSeeder::Init(PHCompositeNode* /*unused*/)
0060 {
0061   return Fun4AllReturnCodes::EVENT_OK;
0062 }
0063
0064 //____________________________________________________________________________..
0065 int PHCosmicSeeder::InitRun(PHCompositeNode* topNode)
0066 {
0067   int ret = getNodes(topNode);
0068   if (ret != Fun4AllReturnCodes::EVENT_OK)
0069   {
0070     return ret;
0071   }
0072   ret = createNodes(topNode);
0073
0074   if (m_analysis)
0075   {
0076     m_outfile = new TFile("PHCosmicSeeder.root", "recreate");
0077     m_tup = new TNtuple("ntp_seed", "seed",
0078                         "event:seed:nclus:xyint:xyslope:xzint:xzslope:"
0079                         "longestxyint:longestxyslope:longestxzint:longestxzslope");
0080   }
0081
0082   return ret;
0083 }
0084
0085 //____________________________________________________________________________..
0086 int PHCosmicSeeder::process_event(PHCompositeNode* /*unused*/)
0087 {
0088   PHCosmicSeeder::PositionMap clusterPositions;
0089   for (const auto& hitsetkey : m_clusterContainer->getHitSetKeys(m_trackerId))
0090   {
0091     auto range = m_clusterContainer->getClusters(hitsetkey);
0092     for (auto citer = range.first; citer != range.second; ++citer)
0093     {
0094       const auto ckey = citer->first;
0095       const auto cluster = citer->second;
0096       if(cluster->getMaxAdc() < m_adcCut){
0097         continue;
0098       }
0099       const auto global = m_tGeometry->getGlobalPosition(ckey, cluster);
0100       clusterPositions.insert(std::make_pair(ckey, global));
0101     }
0102   }
0103   if (Verbosity() > 2)
0104   {
0105     std::cout << "cluster map size is " << clusterPositions.size() << std::endl;
0106   }
0107   auto seeds = makeSeeds(clusterPositions);
0108
0109   if (Verbosity() > 1)
0110   {
0111     std::cout << "Initial seed candidate size is " << seeds.size() << std::endl;
0112   }
0113   std::sort(seeds.begin(), seeds.end(),
0114             [](const seed& a, const seed& b)
0115             { return a.ckeys.size() > b.ckeys.size(); });
0116
0117   auto prunedSeeds = combineSeeds(seeds, clusterPositions);
0118   if (Verbosity() > 1)
0119   {
0120     std::cout << "Pruned seed size is " << prunedSeeds.size() << std::endl;
0121   }
0122   std::sort(prunedSeeds.begin(), prunedSeeds.end(),
0123             [](const seed& a, const seed& b)
0124             { return a.ckeys.size() > b.ckeys.size(); });
0125
0126   auto finalSeeds = findIntersections(prunedSeeds);
0127   if (Verbosity() > 1)
0128   {
0129     std::cout << "final seeds are " << finalSeeds.size() << std::endl;
0130   }
0131   for (auto& seed1 : finalSeeds)
0132   {
0133     recalculateSeedLineParameters(seed1, clusterPositions, true);
0134
0135     recalculateSeedLineParameters(seed1, clusterPositions, false);
0136   }
0137   auto chainedSeeds = chainSeeds(finalSeeds, clusterPositions);
0138   // auto chainedSeeds = finalSeeds;
0139   if (Verbosity() > 1)
0140   {
0141     std::cout << "Total seeds found is " << chainedSeeds.size() << std::endl;
0142   }
0143   int iseed = 0;
0144
0145   auto longestseedit = chainedSeeds.begin();
0146   seed longestseed;
0147   if (longestseedit != chainedSeeds.end())
0148   {
0149     longestseed = chainedSeeds.front();
0150   }
0151   for (auto& seed_A : chainedSeeds)
0152   {
0153     // if mvtx only, we are interested only in seeds with > 3 hits
0154     if (m_trackerId == TrkrDefs::TrkrId::mvtxId && seed_A.ckeys.size() < 4)
0155     {
0156       continue;
0157     }
0158     if (m_analysis)
0159     {
0160       float seed_data[] = {
0161           (float) m_event,
0162           (float) iseed,
0163           (float) seed_A.ckeys.size(),
0164           seed_A.xyintercept,
0165           seed_A.xyslope,
0166           seed_A.xzintercept,
0167           seed_A.xzslope,
0168           longestseed.xyintercept,
0169           longestseed.xyslope,
0170           longestseed.xzintercept,
0171           longestseed.xzslope};
0172       m_tup->Fill(seed_data);
0173     }
0174     auto svtxseed = std::make_unique<TrackSeed_v2>();
0175     for (auto& key : seed_A.ckeys)
0176     {
0177       svtxseed->insert_cluster_key(key);
0178     }
0179     // if (m_trackerId == TrkrDefs::TrkrId::mvtxId){
0180     //   svtxseed->circleFitByTaubin(clusterPositions, 0, 3);
0181     //   svtxseed->lineFit(clusterPositions, 0, 3);
0182     // }
0183     m_seedContainer->insert(svtxseed.get());
0184     ++iseed;
0185   }
0186   if (Verbosity() > 1)
0187   {
0188     std::cout << "Final n seeds: " << m_seedContainer->size() << std::endl;
0189   }
0190   ++m_event;
0191   return Fun4AllReturnCodes::EVENT_OK;
0192 }
0193 PHCosmicSeeder::SeedVector PHCosmicSeeder::findIntersections(PHCosmicSeeder::SeedVector& initialSeeds)
0194 {
0195   std::set<int> seedsToDelete;
0196   //! combine seeds with common ckeys
0197   for (unsigned int i = 0; i < initialSeeds.size(); ++i)
0198   {
0199     auto& seed1 = initialSeeds[i];
0200     for (unsigned int j = i; j < initialSeeds.size(); ++j)
0201     {
0202       if (i == j)
0203       {
0204         continue;
0205       }
0206       auto& seed2 = initialSeeds[j];
0207       std::vector<TrkrDefs::cluskey> intersection;
0208       unsigned int intersection_size_limit = 3;
0209       if (m_trackerId == TrkrDefs::TrkrId::mvtxId)
0210       {
0211         intersection_size_limit = 2;
0212       }
0213       std::set_intersection(seed1.ckeys.begin(), seed1.ckeys.end(),
0214                             seed2.ckeys.begin(), seed2.ckeys.end(), std::back_inserter(intersection));
0215       if (intersection.size() > intersection_size_limit)
0216       {
0217         //! If they share at least 4 (3 for MVTX only) clusters they are likely the same track,
0218         //! so merge and delete
0219         for (auto key : seed2.ckeys)
0220         {
0221           seed1.ckeys.insert(key);
0222         }
0223         seedsToDelete.insert(j);
0224       }
0225     }
0226   }
0227
0228   SeedVector finalSeeds;
0229   unsigned int clus_per_seed_size = 4;
0230   if (m_trackerId == TrkrDefs::TrkrId::mvtxId)
0231   {
0232     clus_per_seed_size = 3;
0233   }
0234   for (unsigned int i = 0; i < initialSeeds.size(); ++i)
0235   {
0236     if (seedsToDelete.find(i) != seedsToDelete.end())
0237     {
0238       continue;
0239     }
0240     if (initialSeeds[i].ckeys.size() < clus_per_seed_size)
0241     {
0242       continue;
0243     }
0244     finalSeeds.push_back(initialSeeds[i]);
0245   }
0246   return finalSeeds;
0247 }
0248 PHCosmicSeeder::SeedVector PHCosmicSeeder::chainSeeds(PHCosmicSeeder::SeedVector& initialSeeds,
0249                                                       PositionMap& clusterPositions)
0250 {
0251   PHCosmicSeeder::SeedVector returnseeds;
0252   std::set<int> seedsToDelete;
0253   if (Verbosity() > 1)
0254   {
0255     std::cout << "Chaining seeds, n seeds =  " << initialSeeds.size() << std::endl;
0256   }
0257   for (unsigned int i = 0; i < initialSeeds.size(); ++i)
0258   {
0259     auto& seed1 = initialSeeds[i];
0260     for (unsigned int j = i; j < initialSeeds.size(); ++j)
0261     {
0262       if (i == j)
0263       {
0264         continue;
0265       }
0266       auto& seed2 = initialSeeds[j];
0267       recalculateSeedLineParameters(seed1, clusterPositions, true);
0268       recalculateSeedLineParameters(seed2, clusterPositions, true);
0269       recalculateSeedLineParameters(seed1, clusterPositions, false);
0270       recalculateSeedLineParameters(seed2, clusterPositions, false);
0271
0272       float longestxyslope = seed1.xyslope;
0273       float longestxyint = seed1.xyintercept;
0274       float longestxzslope = seed1.xzslope;
0275       float longestxzint = seed1.xzintercept;
0276       if (seed1.ckeys.size() < seed2.ckeys.size())
0277       {
0278         longestxyint = seed2.xyintercept;
0279         longestxyslope = seed2.xyslope;
0280         longestxzint = seed2.xzintercept;
0281         longestxzslope = seed2.xzslope;
0282       }
0283
0284       float pdiff_tol = 1.0;
0285       if (m_trackerId == TrkrDefs::TrkrId::mvtxId)
0286       {
0287         pdiff_tol = 0.25;
0288       }
0289       float const pdiff = std::fabs((seed1.xyslope - seed2.xyslope) / longestxyslope);
0290       float const pdiff2 = std::fabs((seed1.xyintercept - seed2.xyintercept) / longestxyint);
0291       float const pdiff3 = std::fabs((seed1.xzintercept - seed2.xzintercept) / longestxzint);
0292       float const pdiff4 = std::fabs((seed1.xzslope - seed2.xzslope) / longestxzslope);
0293       if (Verbosity() > 1)
0294       {
0295         std::cout << "pdiff1,2,3,4 " << pdiff << ", " << pdiff2 << ", " << pdiff3 << ", " << pdiff4 << std::endl;
0296       }
0297       if (pdiff < pdiff_tol && pdiff2 < pdiff_tol && pdiff3 < pdiff_tol && pdiff4 < pdiff_tol)
0298       {
0299         seedsToDelete.insert(j);
0300         for (auto& key : seed2.ckeys)
0301         {
0302           seed1.ckeys.insert(key);
0303         }
0304       }
0305     }
0306   }
0307   for (unsigned int i = 0; i < initialSeeds.size(); ++i)
0308   {
0309     if (seedsToDelete.find(i) != seedsToDelete.end())
0310     {
0311       continue;
0312     }
0313
0314     // if mvtx only, require each seed have at least one hit in each layer
0315     if (m_trackerId == TrkrDefs::TrkrId::mvtxId)
0316     {
0317       std::vector<bool> contains_layer = {false, false, false};
0318       for (auto& key : initialSeeds[i].ckeys)
0319       {
0320         contains_layer[TrkrDefs::getLayer(key)] = true;
0321       }
0322       if (std::find(contains_layer.begin(), contains_layer.end(), false) != contains_layer.end())
0323       {
0324         continue;
0325       }
0326     }
0327
0328     returnseeds.push_back(initialSeeds[i]);
0329   }
0330   if (returnseeds.size() > 1 && m_trackerId == TrkrDefs::TrkrId::mvtxId)
0331   {
0332     PHCosmicSeeder::SeedVector emptyseeds;
0333     return emptyseeds;
0334   }
0335
0336   return returnseeds;
0337 }
0338 PHCosmicSeeder::SeedVector PHCosmicSeeder::combineSeeds(PHCosmicSeeder::SeedVector& initialSeeds,
0339                                                         PHCosmicSeeder::PositionMap& clusterPositions)
0340 {
0341   SeedVector prunedSeeds;
0342   std::set<int> seedsToDelete;
0343
0344   for (unsigned int i = 0; i < initialSeeds.size(); ++i)
0345   {
0346     for (unsigned int j = i; j < initialSeeds.size(); ++j)
0347     {
0348       if (i == j)
0349       {
0350         continue;
0351       }
0352       auto& seed1 = initialSeeds[i];
0353       auto& seed2 = initialSeeds[j];
0354
0355       // recalculate seed parameters
0356       recalculateSeedLineParameters(seed1, clusterPositions, true);
0357       recalculateSeedLineParameters(seed2, clusterPositions, true);
0358       recalculateSeedLineParameters(seed1, clusterPositions, false);
0359       recalculateSeedLineParameters(seed2, clusterPositions, false);
0360
0361       float slope_tol = 0.1;
0362       float incept_tol = 3.0;
0363       if (m_trackerId == TrkrDefs::TrkrId::mvtxId)
0364       {
0365         slope_tol = 0.1;
0366         incept_tol = 0.5;
0367       }
0368       if (Verbosity() > 4)
0369       {
0370         std::cout << "xy slope diff: " << seed1.xyslope - seed2.xyslope << ", xy int diff   " << seed1.xyintercept - seed2.xyintercept << ", xz slope diff " << seed1.xzslope - seed2.xzslope << ", xz int diff   " << seed1.xzintercept - seed2.xzintercept << std::endl;
0371       }
0372       //! These values are tuned on the cosmic data
0373       if (std::fabs(seed1.xyslope - seed2.xyslope) < slope_tol &&
0374           std::fabs(seed1.xyintercept - seed2.xyintercept) < incept_tol &&
0375           std::fabs(seed1.xzslope - seed2.xzslope) < slope_tol &&
0376           std::fabs(seed1.xzintercept - seed2.xzintercept) < incept_tol)
0377       {
0378         for (auto& key : seed2.ckeys)
0379         {
0380           seed1.ckeys.insert(key);
0381         }
0382         seedsToDelete.insert(j);
0383       }
0384     }
0385   }
0386   if (Verbosity() > 4)
0387   {
0388     std::cout << "seeds to delete size is " << seedsToDelete.size() << std::endl;
0389   }
0390   for (unsigned int i = 0; i < initialSeeds.size(); ++i)
0391   {
0392     if (seedsToDelete.find(i) != seedsToDelete.end())
0393     {
0394       continue;
0395     }
0396     prunedSeeds.push_back(initialSeeds[i]);
0397   }
0398
0399   return prunedSeeds;
0400 }
0401 PHCosmicSeeder::SeedVector
0402 PHCosmicSeeder::makeSeeds(PHCosmicSeeder::PositionMap& clusterPositions)
0403 {
0404   PHCosmicSeeder::SeedVector seeds;
0405   std::set<TrkrDefs::cluskey> keys;
0406   //  int seednum = 0;
0407   double dist_check = 2.0;
0408   if (m_trackerId == TrkrDefs::TrkrId::mvtxId)
0409   {
0410     dist_check = 2.5;
0411   }
0412   for (auto& [key1, pos1] : clusterPositions)
0413   {
0414     for (auto& [key2, pos2] : clusterPositions)
0415     {
0416       if (key1 == key2)
0417       {
0418         continue;
0419       }
0420       // make a cut on clusters to at least be close to each other within a few cm
0421       float const dist = (pos2 - pos1).norm();
0422       if (m_trackerId == TrkrDefs::TrkrId::mvtxId && (TrkrDefs::getLayer(key1) == TrkrDefs::getLayer(key2)))
0423       {
0424         continue;
0425       }
0426       if (dist > dist_check)
0427       {
0428         continue;
0429       }
0430       if (Verbosity() > 5)
0431       {
0432         std::cout << "checking keys " << key1 << ", " << key2 << " with dist apart " << dist << std::endl;
0433         std::cout << "positions are " << pos1.transpose() << "    ,     "
0434                   << pos2.transpose() << std::endl;
0435       }
0436       PHCosmicSeeder::seed doub;
0437
0438       doub.xyslope = (pos2.y() - pos1.y()) / (pos2.x() - pos1.x());
0439       doub.xyintercept = pos1.y() - doub.xyslope * pos1.x();
0440       doub.xzslope = (pos2.z() - pos1.z()) / (pos2.x() - pos1.x());
0441       doub.xzintercept = pos1.z() - pos1.x() * doub.xzslope;
0442       doub.yzslope = (pos2.z() - pos1.z()) / (pos2.y() - pos1.y());
0443       doub.yzintercept = pos1.z() - pos1.y() * doub.yzslope;
0444
0445       keys.insert(key1);
0446       keys.insert(key2);
0447       doub.ckeys = keys;
0448       keys.clear();
0449       seeds.push_back(doub);
0450       //      seednum++;
0451     }
0452   }
0453   if (Verbosity() > 2)
0454   {
0455     std::cout << "doublet sizes " << seeds.size() << std::endl;
0456   }
0457   for (auto& dub : seeds)
0458   {
0459     if (Verbosity() > 2)
0460     {
0461       std::cout << "doublet has " << dub.ckeys.size() << " keys " << std::endl;
0462       for (auto key : dub.ckeys)
0463       {
0464         std::cout << "position is " << clusterPositions.find(key)->second.transpose() << std::endl;
0465       }
0466     }
0467     auto begin = dub.ckeys.begin();
0468     auto pos1 = clusterPositions.find(*(begin))->second;
0469     std::advance(begin, 1);
0470     auto pos2 = clusterPositions.find(*(begin))->second;
0471
0472     for (auto& [key, pos] : clusterPositions)
0473     {
0474       //! skip existing keys
0475       if (std::find(dub.ckeys.begin(), dub.ckeys.end(), key) != dub.ckeys.end())
0476       {
0477         continue;
0478       }
0479       // only look at the cluster that is within 2cm of the doublet clusters
0480       float const dist1 = (pos1 - pos).norm();
0481       float const dist2 = (pos2 - pos).norm();
0482       float dist12_check = 2.;
0483       if (m_trackerId == TrkrDefs::TrkrId::mvtxId)
0484       {
0485         dist12_check = 3.5;
0486       }
0487       if (dist1 < dist2)
0488       {
0489         if (dist1 > dist12_check)
0490         {
0491           continue;
0492         }
0493       }
0494       else
0495       {
0496         if (dist2 > dist12_check)
0497         {
0498           continue;
0499         }
0500       }
0501
0502       float const predy = dub.xyslope * pos.x() + dub.xyintercept;
0503       float const predz = dub.xzslope * pos.x() + dub.xzintercept;
0504       float const predz2 = dub.yzslope * pos.y() + dub.yzintercept;
0505       if (Verbosity() > 2)
0506       {
0507         std::cout << "testing ckey " << key << " with box dca "
0508                   << predy << ", " << pos.transpose() << " and " << predz << ", " << pos.z()
0509                   << std::endl;
0510       }
0511       if (fabs(predy - pos.y()) < m_xyTolerance)
0512       {
0513         if (m_trackerId == TrkrDefs::TrkrId::mvtxId && (fabs(predz - pos.z()) > 0.3 || fabs(predz2 - pos.z()) > 0.3))
0514         {
0515           continue;
0516         }
0517         if (Verbosity() > 2)
0518         {
0519           std::cout << "   adding ckey " << key << " with box dca "
0520                     << predy << ", " << pos.y() << " and " << predz << ", " << pos.z()
0521                     << std::endl;
0522           std::cout << "       with pos " << pos.x() << ", " << pos.y() << ", " << pos.z()
0523                     << std::endl;
0524         }
0525         dub.ckeys.insert(key);
0526       }
0527     }
0528   }
0529
0530   /// erase doublets
0531   std::set<int> seedsToDelete;
0532   for (unsigned int i = 0; i < seeds.size(); ++i)
0533   {
0534     auto seed_A = seeds[i];
0535     if (Verbosity() > 2 && seed_A.ckeys.size() > 2)
0536     {
0537       std::cout << "keys in seed " << std::endl;
0538       for (auto& key : seed_A.ckeys)
0539       {
0540         std::cout << key << ", ";
0541       }
0542       std::cout << "seed xy slope: " << seed_A.xyslope << std::endl;
0543       std::cout << std::endl
0544                 << "seed pos " << std::endl;
0545       for (auto& key : seed_A.ckeys)
0546       {
0547         std::cout << clusterPositions.find(key)->second.transpose() << std::endl;
0548       }
0549       std::cout << "Done printing seed info" << std::endl;
0550     }
0551     if (seed_A.ckeys.size() < 3)
0552     {
0553       seedsToDelete.insert(i);
0554     }
0555   }
0556
0557   PHCosmicSeeder::SeedVector returnSeeds;
0558   for (unsigned int i = 0; i < seeds.size(); i++)
0559   {
0560     if (seedsToDelete.find(i) != seedsToDelete.end())
0561     {
0562       continue;
0563     }
0564     returnSeeds.push_back(seeds[i]);
0565   }
0566   return returnSeeds;
0567 }
0568 void PHCosmicSeeder::recalculateSeedLineParameters(seed& seed_A,
0569                                                    PHCosmicSeeder::PositionMap& clusters, bool isXY)
0570 {
0571   float avgx = 0;
0572   float avgy = 0;
0573   PHCosmicSeeder::PositionMap seedClusters;
0574   for (auto& key : seed_A.ckeys)
0575   {
0576     auto glob = clusters.find(key)->second;
0577     if (isXY)
0578     {
0579       avgx += glob.x();
0580       avgy += glob.y();
0581     }
0582     else
0583     {
0584       avgx += glob.x();
0585       avgy += glob.z();
0586     }
0587     seedClusters.insert(std::make_pair(key, glob));
0588   }
0589
0590   avgx /= seed_A.ckeys.size();
0591   avgy /= seed_A.ckeys.size();
0592   float num = 0;
0593   float denom = 0;
0594   for (auto& [key, glob] : seedClusters)
0595   {
0596     if (isXY)
0597     {
0598       num += (glob.x() - avgx) * (glob.y() - avgy);
0599       denom += square(glob.x() - avgx);
0600     }
0601     else
0602     {
0603       num += (glob.x() - avgx) * (glob.z() - avgy);
0604       denom += square(glob.x() - avgx);
0605     }
0606   }
0607   if (isXY)
0608   {
0609     seed_A.xyslope = num / denom;
0610     seed_A.xyintercept = avgy - seed_A.xyslope * avgx;
0611   }
0612   else
0613   {
0614     seed_A.xzslope = num / denom;
0615     seed_A.xzintercept = avgy - seed_A.xzslope * avgx;
0616   }
0617 }
0618 int PHCosmicSeeder::getNodes(PHCompositeNode* topNode)
0619 {
0620   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0621   if (!m_tGeometry)
0622   {
0623     std::cout << PHWHERE << "No acts reco geometry, bailing."
0624               << std::endl;
0625     return Fun4AllReturnCodes::ABORTEVENT;
0626   }
0627   m_clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0628   if (!m_clusterContainer)
0629   {
0630     std::cout << PHWHERE << "No cluster container on node tree, bailing."
0631               << std::endl;
0632     return Fun4AllReturnCodes::ABORTEVENT;
0633   }
0634   return Fun4AllReturnCodes::EVENT_OK;
0635 }
0636
0637 int PHCosmicSeeder::createNodes(PHCompositeNode* topNode)
0638 {
0639   PHNodeIterator iter(topNode);
0640
0641   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0642
0643   if (!dstNode)
0644   {
0645     std::cerr << "DST node is missing, quitting" << std::endl;
0646     throw std::runtime_error("Failed to find DST node in PHCosmicSeeder::createNodes");
0647   }
0648
0649   PHNodeIterator dstIter(dstNode);
0650   PHCompositeNode* svtxNode = dynamic_cast<PHCompositeNode*>(dstIter.findFirst("PHCompositeNode", "SVTX"));
0651
0652   if (!svtxNode)
0653   {
0654     svtxNode = new PHCompositeNode("SVTX");
0655     dstNode->addNode(svtxNode);
0656   }
0657
0658   m_seedContainer = findNode::getClass<TrackSeedContainer>(topNode, m_trackMapName);
0659   if (!m_seedContainer)
0660   {
0661     m_seedContainer = new TrackSeedContainer_v1;
0662     PHIODataNode<PHObject>* trackNode =
0663         new PHIODataNode<PHObject>(m_seedContainer, m_trackMapName, "PHObject");
0664     svtxNode->addNode(trackNode);
0665   }
0666
0667   return Fun4AllReturnCodes::EVENT_OK;
0668 }
0669 //____________________________________________________________________________..
0670 int PHCosmicSeeder::End(PHCompositeNode* /*unused*/)
0671 {
0672   if (m_outfile)
0673   {
0674     m_outfile->cd();
0675     m_tup->Write();
0676     m_outfile->Close();
0677   }
0678   return Fun4AllReturnCodes::EVENT_OK;
0679 }