g4simulation/g4eval/DSTTrackInfoWriter.cc

0001 /*!
0002  * \file DSTTrackInfoWriter.cc
0003  * \author Alex Patton <aopatton@mit.edu>
0004  */
0005
0006 #include "DSTTrackInfoWriter.h"
0007
0008 #include <trackbase_historic/SvtxTrackMap.h>
0009 #include <trackbase_historic/TrackInfoContainer_v1.h>
0010 #include <trackbase_historic/ActsTransformations.h>
0011
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013
0014 #include <phool/PHCompositeNode.h>
0015 #include <phool/PHNodeIterator.h>
0016 #include <phool/getClass.h>
0017
0018 #include <iostream>
0019
0020 //_____________________________________________________________________
0021 DSTTrackInfoWriter::DSTTrackInfoWriter(const std::string& name)
0022   : SubsysReco(name)
0023 {
0024 }
0025
0026 //_____________________________________________________________________
0027 int DSTTrackInfoWriter::InitRun(PHCompositeNode* topNode)
0028 {
0029   if (Verbosity() > 1)
0030   {
0031     std::cout << "Writer Init start" << std::endl;
0032   }
0033   // find DST node
0034   PHNodeIterator iter(topNode);
0035   auto *dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0036   if (!dstNode)
0037   {
0038     std::cout << "DSTTrackInfoWriter::Init - DST Node missing" << std::endl;
0039     return Fun4AllReturnCodes::ABORTEVENT;
0040   }
0041
0042   // get EVAL node
0043   iter = PHNodeIterator(dstNode);
0044   auto *evalNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "EVAL"));
0045   if (!evalNode)
0046   {
0047     // create
0048     std::cout << "DSTTrackInfoWriter::Init - EVAL node missing - creating" << std::endl;
0049     evalNode = new PHCompositeNode("EVAL");
0050     dstNode->addNode(evalNode);
0051   }
0052
0053   // // TClonesArary
0054   // m_container->arrClsDST = new TClonesArray("DSTContainerv3::ClusterStruct");
0055   // m_container->trkrClsDST = new TClonesArray("TrkrClusterv4");
0056
0057   auto *newInfoNode = new PHIODataNode<PHObject>(new TrackInfoContainer_v1, "TrackInfoContainer", "PHObject");
0058   evalNode->addNode(newInfoNode);
0059
0060   // auto newTrackNode = new PHIODataNode<PHObject>(new SvtxTrackArray_v1, "TRACK_ARRAYV1", "PHObject");
0061   // evalNode->addNode(newTrackNode);
0062
0063   // TClonesArary
0064   // fcl = new TFile("dstcl.root", "recreate");
0065   // tcl = new TTree("tcl", "dst tree");
0066   // arrEvt = new TClonesArray("DSTContainerv3::EventStruct");
0067   // arrTrk = new TClonesArray("DSTContainerv3::TrackStruct");
0068   // arrCls = new TClonesArray("DSTContainerv3::ClusterStruct");
0069   // tcl->Branch("evt", &arrEvt);
0070   // tcl->Branch("trk", &arrTrk);
0071   // tcl->Branch("cls", &arrCls);
0072   if (Verbosity() > 1)
0073   {
0074     std::cout << "Writer Init end" << std::endl;
0075   }
0076   return Fun4AllReturnCodes::EVENT_OK;
0077 }
0078
0079 //_____________________________________________________________________
0080 int DSTTrackInfoWriter::process_event(PHCompositeNode* topNode)
0081 {
0082   // make topNode run in Init
0083   // Init(topNode);
0084   //  load nodes
0085   if (Verbosity() > 1)
0086   {
0087     std::cout << __FILE__ << "::" << __func__ << "::" << __LINE__ << std::endl;
0088     std::cout << "DSTTrackInfoWriter::process_event" << std::endl;
0089   }
0090   auto res = load_nodes(topNode);
0091   if (res != Fun4AllReturnCodes::EVENT_OK)
0092   {
0093     return res;
0094   }
0095   if (Verbosity() > 1)
0096   {
0097     std::cout << "Return codes  start" << Fun4AllReturnCodes::EVENT_OK << std::endl;
0098   }
0099   // cleanup output container
0100   if (m_track_info_container) { m_track_info_container->Reset();
0101 }
0102   if (Verbosity() > 1)
0103   {
0104     std::cout << "Evalutate track info" << std::endl;
0105   }
0106   evaluate_track_info();
0107   if (Verbosity() > 1)
0108   {
0109     std::cout << "exiting event"
0110               << "\n";
0111
0112     std::cout << "Return codes end" << Fun4AllReturnCodes::EVENT_OK << std::endl;
0113   }
0114   return Fun4AllReturnCodes::EVENT_OK;
0115 }
0116
0117 //_____________________________________________________________________
0118 int DSTTrackInfoWriter::load_nodes(PHCompositeNode* topNode)
0119 {
0120   // get necessary nodes
0121   m_track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0122
0123   m_track_info_container = findNode::getClass<TrackInfoContainer>(topNode, "TrackInfoContainer");
0124
0125   return Fun4AllReturnCodes::EVENT_OK;
0126 }
0127
0128 //____________________________________________________________________
0129 void DSTTrackInfoWriter::evaluate_track_info()
0130 {
0131   if (!(m_track_info_container))
0132   {
0133     return;
0134   }
0135
0136   m_track_info_container->Reset();
0137   // get track into track info
0138   Int_t iTrk = 0;
0139   // long unsigned int iKey = 0;
0140   if (Verbosity() > 1)
0141   {
0142     std::cout << "Before loop"
0143               << "\n";
0144   }
0145   ActsTransformations transformer;
0146   for (const auto& trackpair : *m_track_map)
0147   {
0148     auto *const track = trackpair.second;
0149     // this track will have a TPC and Silicon seed
0150
0151     uint64_t hitbitmap = 0;
0152
0153     SvtxTrackInfo *trackInfo = new SvtxTrackInfo_v1();
0154     if (Verbosity() > 1)
0155     {
0156       std::cout << "Before seeds"
0157                 << "\n";
0158     }
0159     TrackSeed* TPCSeed = track->get_tpc_seed();
0160     TrackSeed* SiliconSeed = track->get_silicon_seed();
0161     if (TPCSeed)
0162     {
0163       for (auto key_iter = TPCSeed->begin_cluster_keys(); key_iter != TPCSeed->end_cluster_keys(); ++key_iter)
0164       {
0165         const auto& cluster_key = *key_iter;
0166
0167         // store information in track array
0168         // std::cout << "TPC clusterkey: " << cluster_key <<"\n";
0169         // std::cout << "TPC subsurfkey: " << cluster->getSubSurfKey() << std::endl;
0170         uint8_t layer = TrkrDefs::getLayer(cluster_key);
0171         if (Verbosity() > 1)
0172         {
0173           std::cout << "Layer is: " << unsigned(layer) << std::endl;
0174         }
0175         hitbitmap = hitbitmap + ((uint64_t) 1 << layer);
0176
0177         // TrkrDefs::
0178       }
0179     }
0180     if (Verbosity() > 1)
0181     {
0182       std::cout << "Before Silicon seeds"
0183                 << "\n";
0184     }
0185     if (!SiliconSeed)
0186     {
0187       if (Verbosity() > 1)
0188       {
0189         std::cout << "Silicon Seed does not exist" << std::endl;
0190       }
0191     }
0192
0193     if (SiliconSeed)
0194     {
0195       for (auto key_iter = SiliconSeed->begin_cluster_keys(); key_iter != SiliconSeed->end_cluster_keys(); ++key_iter)
0196       {
0197         const auto& cluster_key = *key_iter;
0198
0199         uint8_t layer = TrkrDefs::getLayer(cluster_key);
0200         if (Verbosity() > 1)
0201         {
0202           std::cout << "Layer is: " << unsigned(layer) << std::endl;
0203         }
0204         hitbitmap = hitbitmap + ((uint64_t) 1 << layer);
0205       }
0206     }
0207     if (Verbosity() > 1)
0208     {
0209       std::cout << "After Track seeds"
0210                 << "\n";
0211     }
0212     trackInfo->set_hitbitmap(hitbitmap);
0213     trackInfo->set_chisq(track->get_chisq());
0214     trackInfo->set_ndf(track->get_ndf());
0215     trackInfo->set_crossing(track->get_crossing());
0216
0217     trackInfo->set_x(track->get_x());
0218     trackInfo->set_y(track->get_y());
0219     trackInfo->set_z(track->get_z());
0220     trackInfo->set_phi(track->get_phi());
0221     trackInfo->set_theta(2.* std::atan(std::exp(-1*track->get_eta())));
0222     trackInfo->set_qOp(track->get_charge() / track->get_p());
0223     if (Verbosity() > 1)
0224     {
0225       std::cout << "track crossing: " << track->get_crossing() << std::endl;
0226
0227       std::cout << "track.get_z(): " << track->get_z() << std::endl;
0228       std::cout << "trackInfo.get_z(): " << trackInfo->get_z() << std::endl;
0229       std::cout << "Hitbitmap: " << trackInfo->get_hitbitmap() << std::endl;
0230       std::cout << "crossing: " << trackInfo->get_crossing() << std::endl;
0231       std::cout << "chi^2: " << trackInfo->get_chisq() << std::endl;
0232       std::cout << "ndf: " << unsigned(trackInfo->get_ndf()) << std::endl;
0233     }
0234     auto actsMatrix = transformer.rotateSvtxTrackCovToActs(track);
0235     //! only take the first 5 entries, which correspond to d0, z0, phi, theta, q/p
0236     for (int i = 0; i < 5; i++)
0237     {
0238       for (int j = i; j < 5; j++)
0239       {
0240         trackInfo->set_covariance(i, j, actsMatrix(i,j));
0241       }
0242     }
0243     if (Verbosity() > 1)
0244     {
0245       std::cout << "Right before adding track info" << iTrk << std::endl;
0246     }
0247     m_track_info_container->add_trackinfo(iTrk, trackInfo);
0248     if (Verbosity() > 1)
0249     {
0250       std::cout << "Right after adding track info" << std::endl;
0251     }
0252     delete trackInfo;
0253     ++iTrk;
0254   }
0255
0256   // add trackinfo to trackinfocontainer
0257 }
0258 //_____________________________________________________________________
0259
0260 //_____________________________________________________________________