packages/rawtodst/tpc_hits.cc

0001
0002 #include "tpc_hits.h"
0003
0004 #include <trackbase/TpcDefs.h>
0005 #include <trackbase/TrkrDefs.h>
0006 #include <trackbase/TrkrHit.h>
0007 #include <trackbase/TrkrHitSet.h>
0008 #include <trackbase/TrkrHitSetContainer.h>
0009 #include <trackbase/TrkrHitSetContainerv1.h>
0010 #include <trackbase/TrkrHitSetv1.h>
0011 #include <trackbase/TrkrHitv2.h>
0012
0013 #include <fun4all/Fun4AllHistoManager.h>
0014 #include <fun4all/Fun4AllReturnCodes.h>
0015
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/PHIODataNode.h>    // for PHIODataNode
0018 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0019 #include <phool/PHObject.h>        // for PHObject
0020 #include <phool/getClass.h>
0021 #include <phool/phool.h>
0022
0023 #include <Event/Event.h>
0024 #include <Event/EventTypes.h>
0025 #include <Event/packet.h>
0026
0027 #include <TFile.h>
0028 #include <TH2.h>
0029
0030 //____________________________________________________________________________..
0031 tpc_hits::tpc_hits(const std::string &name)
0032   : SubsysReco(name)
0033   , hm(nullptr)
0034   , _filename("./outputfile.root")
0035 {
0036   std::cout << "tpc_hits::tpc_hits(const std::string &name)" << std::endl;
0037   starting_BCO = -1;
0038   rollover_value = 0;
0039   current_BCOBIN = 0;
0040   M.setMapNames("AutoPad-R1-RevA.sch.ChannelMapping.csv", "AutoPad-R2-RevA-Pads.sch.ChannelMapping.csv", "AutoPad-R3-RevA.sch.ChannelMapping.csv");
0041 }
0042
0043 //____________________________________________________________________________..
0044 tpc_hits::~tpc_hits()
0045 {
0046   delete hm;
0047   std::cout << "tpc_hits::~tpc_hits() Calling dtor" << std::endl;
0048 }
0049
0050 //____________________________________________________________________________..
0051 int tpc_hits::Init(PHCompositeNode * /*topNode*/)
0052 {
0053   std::cout << "tpc_hits::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0054   hm = new Fun4AllHistoManager("HITHIST");
0055
0056   _h_hit_XY = new TH2F("_h_hit_XY", "_h_hit_XY;X, [m];Y, [m]", 400, -800, 800, 400, -800, 800);
0057
0058   hm->registerHisto(_h_hit_XY);
0059
0060   return Fun4AllReturnCodes::EVENT_OK;
0061 }
0062
0063 //____________________________________________________________________________..
0064 int tpc_hits::InitRun(PHCompositeNode *topNode)
0065 {
0066   // get dst node
0067   PHNodeIterator iter(topNode);
0068   auto dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0069   if (!dstNode)
0070   {
0071     std::cout << "tpc_hits::InitRun - DST Node missing, doing nothing." << std::endl;
0072     exit(1);
0073   }
0074
0075   // create hitset container if needed
0076   auto hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0077   if (!hitsetcontainer)
0078   {
0079     std::cout << "tpc_hits::InitRun - creating TRKR_HITSET." << std::endl;
0080
0081     // find or create TRKR node
0082     PHNodeIterator dstiter(dstNode);
0083     auto trkrnode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0084     if (!trkrnode)
0085     {
0086       trkrnode = new PHCompositeNode("TRKR");
0087       dstNode->addNode(trkrnode);
0088     }
0089
0090     // create container and add to the tree
0091     hitsetcontainer = new TrkrHitSetContainerv1;
0092     auto newNode = new PHIODataNode<PHObject>(hitsetcontainer, "TRKR_HITSET", "PHObject");
0093     trkrnode->addNode(newNode);
0094   }
0095   topNode->print();
0096
0097   // we reset the BCO for the new run
0098   starting_BCO = -1;
0099   rollover_value = 0;
0100   current_BCOBIN = 0;
0101
0102   // m_hits = new TrkrHitSetContainerv1();
0103
0104   return Fun4AllReturnCodes::EVENT_OK;
0105 }
0106
0107 //____________________________________________________________________________..
0108 int tpc_hits::process_event(PHCompositeNode *topNode)
0109 {
0110   // std::cout << "tpc_hits::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0111   //  load relevant nodes
0112   //  Get the TrkrHitSetContainer node
0113   auto trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0114   assert(trkrhitsetcontainer);
0115
0116   Event *_event = findNode::getClass<Event>(topNode, "PRDF");
0117   assert(_event);
0118
0119   if (_event == nullptr)
0120   {
0121     std::cout << "tpc_hits::Process_Event - Event not found" << std::endl;
0122     return -1;
0123   }
0124   if (_event->getEvtType() >= 8)  /// special events
0125   {
0126     return Fun4AllReturnCodes::DISCARDEVENT;
0127   }
0128
0129   // check all possible TPC packets that we need to analyze
0130   for (int ep = 0; ep < 2; ep++)
0131   {
0132     for (int sector = 0; sector <= 23; sector++)
0133     {
0134       const int packet = 4000 + sector * 10 + ep;
0135
0136       // Reading packet
0137       Packet *p = _event->getPacket(packet);
0138
0139       // Figure out which side
0140       int side = 0;
0141       if (sector > 11)
0142       {
0143         side = 1;
0144       }
0145
0146       if (p)
0147       {
0148         std::cout << "tpc_hits:: Event getPacket: " << packet << "| Sector" << sector << "| EndPoint " << ep << std::endl;
0149       }
0150       else
0151       {
0152         continue;
0153       }
0154
0155       int nr_of_waveforms = p->iValue(0, "NR_WF");
0156
0157       for (auto &l : m_hitset)
0158       {
0159         l = new TrkrHitSetv1();
0160
0161         int wf;
0162         for (wf = 0; wf < nr_of_waveforms; wf++)
0163         {
0164           int current_BCO = p->iValue(wf, "BCO") + rollover_value;
0165
0166           if (starting_BCO < 0)
0167           {
0168             starting_BCO = current_BCO;
0169           }
0170
0171           if (current_BCO < starting_BCO)  // we have a rollover
0172           {
0173             rollover_value += 0x100000;
0174             current_BCO = p->iValue(wf, "BCO") + rollover_value;
0175             starting_BCO = current_BCO;
0176             current_BCOBIN++;
0177           }
0178           int sampa_nr = p->iValue(wf, "SAMPAADDRESS");
0179           int channel = p->iValue(wf, "CHANNEL");
0180
0181           int fee = p->iValue(wf, "FEE");
0182           int samples = p->iValue(wf, "SAMPLES");
0183           // clockwise FEE mapping
0184           // int FEE_map[26]={5, 6, 1, 3, 2, 12, 10, 11, 9, 8, 7, 1, 2, 4, 8, 7, 6, 5, 4, 3, 1, 3, 2, 4, 6, 5};
0185           int FEE_R[26] = {2, 2, 1, 1, 1, 3, 3, 3, 3, 3, 3, 2, 2, 1, 2, 2, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3};
0186           // conter clockwise FEE mapping (From Takao)
0187           int FEE_map[26] = {3, 2, 5, 3, 4, 0, 2, 1, 3, 4, 5, 7, 6, 2, 0, 1, 0, 1, 4, 5, 11, 9, 10, 8, 6, 7};
0188           int pads_per_sector[3] = {96, 128, 192};
0189
0190           // setting the mapp of the FEE
0191           int feeM = FEE_map[fee] - 1;
0192           if (FEE_R[fee] == 2)
0193           {
0194             feeM += 6;
0195           }
0196           if (FEE_R[fee] == 3)
0197           {
0198             feeM += 14;
0199           }
0200           int layer = M.getLayer(feeM, channel);
0201
0202           double R = M.getR(feeM, channel);
0203           double phi = M.getPhi(feeM, channel) + sector * M_PI / 6;
0204
0205           TrkrDefs::hitsetkey tpcHitSetKey = TpcDefs::genHitSetKey(layer, sector, side);
0206           TrkrHitSetContainer::Iterator hitsetit = trkrhitsetcontainer->findOrAddHitSet(tpcHitSetKey);
0207
0208           if (Verbosity())
0209           {
0210             int sampaAddress = p->iValue(wf, "SAMPAADDRESS");
0211             int sampaChannel = p->iValue(wf, "SAMPACHANNEL");
0212             int checksum = p->iValue(wf, "CHECKSUM");
0213             int checksumError = p->iValue(wf, "CHECKSUMERROR");
0214             std::cout << "tpc_hits::Process_Event Samples " << samples
0215                       << " Chn:" << channel
0216                       << " layer: " << layer
0217                       << " sampa: " << sampa_nr
0218                       << " fee: " << fee
0219                       << " Mapped fee: " << feeM
0220                       << " sampaAddress: " << sampaAddress
0221                       << " sampaChannel: " << sampaChannel
0222                       << " checksum: " << checksum
0223                       << " checksumError: " << checksumError
0224                       << " hitsetkey " << tpcHitSetKey
0225                       << " R = " << R
0226                       << " phi = " << phi
0227                       << std::endl;
0228           }
0229           for (int s = 0; s < samples; s++)
0230           {
0231             int pad = M.getPad(feeM, channel);
0232             int t = s + 2 * (current_BCO - starting_BCO);
0233             int adc = p->iValue(wf, s);
0234             // generate hit key
0235             TrkrDefs::hitkey hitkey = TpcDefs::genHitKey((unsigned int) pad + sector * pads_per_sector[FEE_R[sector] - 1], (unsigned int) t);
0236             // find existing hit, or create
0237             auto hit = hitsetit->second->getHit(hitkey);
0238
0239             // create hit, assign adc and insert in hitset
0240             if (!hit)
0241             {
0242               // create a new one
0243               hit = new TrkrHitv2();
0244               hit->setAdc(adc);
0245               hitsetit->second->addHitSpecificKey(hitkey, hit);
0246             }
0247             // else{
0248             //   hit->setAdc(adc);
0249             //   hitsetit->second->addHitSpecificKey(hitkey, hit);
0250             // }
0251
0252             _h_hit_XY->Fill(R * cos(phi), R * sin(phi), adc);
0253           }
0254         }
0255       }
0256
0257     }  // End of run over packets
0258   }    // End of ep loop
0259   // we skip the mapping to real pads at first. We just say
0260   // that we get 16 rows (segment R2) with 128 pads
0261   // so each FEE fills 2 rows. Not right, but one step at a time.
0262
0263   return Fun4AllReturnCodes::EVENT_OK;
0264 }
0265
0266 //____________________________________________________________________________..
0267 // int tpc_hits::ResetEvent(PHCompositeNode * /*topNode*/)
0268 //{
0269 //  std::cout << "tpc_hits::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0270 //  return Fun4AllReturnCodes::EVENT_OK;
0271 //}
0272
0273 //____________________________________________________________________________..
0274 // int tpc_hits::EndRun(const int runnumber)
0275 //{
0276 //  std::cout << "tpc_hits::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0277 //  return Fun4AllReturnCodes::EVENT_OK;
0278 //}
0279
0280 //____________________________________________________________________________..
0281 int tpc_hits::End(PHCompositeNode * /*topNode*/)
0282 {
0283   std::cout << "tpc_hits::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0284   hm->dumpHistos(_filename, "RECREATE");
0285
0286   return Fun4AllReturnCodes::EVENT_OK;
0287 }
0288
0289 //____________________________________________________________________________..
0290 // int tpc_hits::Reset(PHCompositeNode * /*topNode*/)
0291 //{
0292 //  std::cout << "tpc_hits::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0293 //  return Fun4AllReturnCodes::EVENT_OK;
0294 //}
0295
0296 //____________________________________________________________________________..
0297 // void tpc_hits::Print(const std::string &what) const
0298 //{
0299 //  std::cout << "tpc_hits::Print(const std::string &what) const Printing info for " << what << std::endl;
0300 //}