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