g4simulation/g4tpc/PHG4TpcPadBaselineShift.cc

0001 #include "PHG4TpcPadBaselineShift.h"
0002
0003 #include <trackbase/TpcDefs.h>
0004
0005 #pragma GCC diagnostic push
0006 #pragma GCC diagnostic ignored "-Wshadow"
0007 #include <trackbase/ActsSurfaceMaps.h>       // for ActsSurfaceMaps
0008 #include <trackbase/ActsTrackingGeometry.h>  // for ActsTrackingG...
0009 #pragma GCC diagnostic pop
0010
0011 #include <trackbase/TrkrClusterContainer.h>  // for TrkrClusterCo...
0012 #include <trackbase/TrkrClusterHitAssoc.h>   // for TrkrClusterHi...
0013 #include <trackbase/TrkrHit.h>               // for TrkrHit
0014
0015 #include <trackbase/TrkrDefs.h>  // for hitkey, getLayer
0016 #include <trackbase/TrkrHitSet.h>
0017 #include <trackbase/TrkrHitSetContainer.h>
0018
0019 #include <fun4all/Fun4AllReturnCodes.h>
0020 #include <fun4all/SubsysReco.h>  // for SubsysReco
0021
0022 #include <g4detectors/PHG4TpcCylinderGeom.h>
0023 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0024
0025 #include <phool/PHCompositeNode.h>
0026 #include <phool/PHNode.h>  // for PHNode
0027 #include <phool/PHNodeIterator.h>
0028 #include <phool/getClass.h>
0029 #include <phool/phool.h>  // for PHWHERE
0030
0031 #include <TF1.h>
0032 #include <TFile.h>
0033 #include <TTree.h>
0034
0035 #include <cmath>  // for sqrt, cos, sin
0036 #include <iostream>
0037 #include <map>  // for _Rb_tree_cons...
0038 #include <string>
0039 #include <utility>  // for pair
0040 #include <vector>
0041
0042 namespace
0043 {
0044   template <class T>
0045   inline constexpr T square(const T &x)
0046   {
0047     return x * x;
0048   }
0049 }  // namespace
0050
0051 int findRBin(float R)
0052 {
0053   // Finding pad number from the center (bin) for hits
0054   int binR = -1;
0055   // Realistic binning
0056   // double r_bins[r_bins_N+1] = {217.83,
0057   //                             311.05,317.92,323.31,329.27,334.63,340.59,345.95,351.91,357.27,363.23,368.59,374.55,379.91,385.87,391.23,397.19,402.49,
0058   //                             411.53,421.70,431.90,442.11,452.32,462.52,472.73,482.94,493.14,503.35,513.56,523.76,533.97,544.18,554.39,564.59,574.76,
0059   //                             583.67,594.59,605.57,616.54,627.51,638.48,649.45,660.42,671.39,682.36,693.33,704.30,715.27,726.24,737.21,748.18,759.11};
0060   const int r_bins_N = 53;
0061   double r_bins[r_bins_N + 1];
0062   r_bins[0] = 30.3125;
0063   double bin_width = 0.625;
0064   for (int i = 1; i < r_bins_N; i++)
0065   {
0066     if (i == 16)
0067     {
0068       bin_width = 0.9375;
0069     }
0070     if (i > 16)
0071     {
0072       bin_width = 1.25;
0073     }
0074     if (i == 31)
0075     {
0076       bin_width = 1.1562;
0077     }
0078     if (i > 31)
0079     {
0080       bin_width = 1.0624;
0081     }
0082
0083     r_bins[i] = r_bins[i - 1] + bin_width;
0084   }
0085
0086   double R_min = 30;
0087   while (R > R_min)
0088   {
0089     binR += 1;
0090     R_min = r_bins[binR];
0091   }
0092   return binR;
0093 }
0094
0095 //____________________________________________________________________________..
0096 PHG4TpcPadBaselineShift::PHG4TpcPadBaselineShift(const std::string &name)
0097   : SubsysReco(name)
0098 {
0099   std::cout << "PHG4TpcPadBaselineShift::PHG4TpcPadBaselineShift(const std::string &name) Calling ctor" << std::endl;
0100 }
0101
0102 bool PHG4TpcPadBaselineShift::is_in_sector_boundary(int phibin, int sector, PHG4TpcCylinderGeom *layergeom)
0103 {
0104   bool reject_it = false;
0105
0106   // sector boundaries occur every 1/12 of the full phi bin range
0107   int PhiBins = layergeom->get_phibins();
0108   int PhiBinsSector = PhiBins / 12;
0109
0110   double radius = layergeom->get_radius();
0111   double PhiBinSize = 2.0 * radius * M_PI / (double) PhiBins;
0112
0113   // sector starts where?
0114   int sector_lo = sector * PhiBinsSector;
0115   int sector_hi = sector_lo + PhiBinsSector - 1;
0116
0117   int sector_fiducial_bins = (int) (SectorFiducialCut / PhiBinSize);
0118
0119   if (phibin < sector_lo + sector_fiducial_bins || phibin > sector_hi - sector_fiducial_bins)
0120   {
0121     reject_it = true;
0122   }
0123
0124   return reject_it;
0125 }
0126 //____________________________________________________________________________..
0127 PHG4TpcPadBaselineShift::~PHG4TpcPadBaselineShift()
0128 {
0129   std::cout << "PHG4TpcPadBaselineShift::~PHG4TpcPadBaselineShift() Calling dtor" << std::endl;
0130 }
0131
0132 //____________________________________________________________________________..
0133 int PHG4TpcPadBaselineShift::Init(PHCompositeNode * /*topNode*/)
0134 {
0135   // outfile = new TFile(_filename.c_str(), "RECREATE");
0136   _hit_z = 0;
0137   _hit_r = 0;
0138   _hit_phi = 0;
0139   _hit_e = 0;
0140   _hit_adc = 0;
0141   _hit_adc_bls = 0;
0142   _hit_layer = -1;
0143   _hit_sector = -1;
0144   if (_writeTree == 1)
0145   {
0146     outfile = new TFile(_filename.c_str(), "RECREATE");
0147     _rawHits = new TTree("hTree", "tpc hit tree for base-line shift tests");
0148     _rawHits->Branch("z", &_hit_z);
0149     _rawHits->Branch("r", &_hit_r);
0150     _rawHits->Branch("phi", &_hit_phi);
0151     _rawHits->Branch("e", &_hit_e);
0152     _rawHits->Branch("adc", &_hit_adc);
0153     _rawHits->Branch("adc_BLS", &_hit_adc_bls);
0154     _rawHits->Branch("hit_layer", &_hit_layer);
0155     _rawHits->Branch("_hit_sector", &_hit_sector);
0156   }
0157   return 0;
0158   // std::cout << "PHG4TpcPadBaselineShift::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0159   // return Fun4AllReturnCodes::EVENT_OK;
0160 }
0161
0162 //____________________________________________________________________________..
0163 int PHG4TpcPadBaselineShift::InitRun(PHCompositeNode *topNode)
0164 {
0165   PHNodeIterator iter(topNode);
0166
0167   // Looking for the DST node
0168   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0169   if (!dstNode)
0170   {
0171     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0172     return Fun4AllReturnCodes::ABORTRUN;
0173   }
0174
0175   auto geom =
0176       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0177   if (!geom)
0178   {
0179     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0180     return Fun4AllReturnCodes::ABORTRUN;
0181   }
0182 AdcClockPeriod = geom->GetFirstLayerCellGeom()->get_zstep();
0183
0184   std::cout << "PHG4TpcPadBaselineShift::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0185   return Fun4AllReturnCodes::EVENT_OK;
0186 }
0187
0188 //____________________________________________________________________________..
0189 int PHG4TpcPadBaselineShift::process_event(PHCompositeNode *topNode)
0190 {
0191   //  int print_layer = 18;
0192
0193   if (Verbosity() > 1000)
0194   {
0195     std::cout << "PHG4TpcPadBaselineShift::Process_Event" << std::endl;
0196   }
0197
0198   PHNodeIterator iter(topNode);
0199   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0200   if (!dstNode)
0201   {
0202     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0203     return Fun4AllReturnCodes::ABORTRUN;
0204   }
0205
0206   // get node containing the digitized hits
0207   m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0208   if (!m_hits)
0209   {
0210     std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
0211     return Fun4AllReturnCodes::ABORTRUN;
0212   }
0213
0214   // get node for clusters
0215   m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0216   if (!m_clusterlist)
0217   {
0218     std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTER." << std::endl;
0219     return Fun4AllReturnCodes::ABORTRUN;
0220   }
0221
0222   // get node for cluster hit associations
0223   m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0224   if (!m_clusterhitassoc)
0225   {
0226     std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERHITASSOC" << std::endl;
0227     return Fun4AllReturnCodes::ABORTRUN;
0228   }
0229
0230   PHG4TpcCylinderGeomContainer *geom_container =
0231       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0232   if (!geom_container)
0233   {
0234     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0235     return Fun4AllReturnCodes::ABORTRUN;
0236   }
0237
0238   m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
0239                                                          "ActsTrackingGeometry");
0240   if (!m_tGeometry)
0241   {
0242     std::cout << PHWHERE
0243               << "ActsTrackingGeometry not found on node tree. Exiting"
0244               << std::endl;
0245     return Fun4AllReturnCodes::ABORTRUN;
0246   }
0247
0248   m_surfMaps = findNode::getClass<ActsSurfaceMaps>(topNode,
0249                                                    "ActsSurfaceMaps");
0250   if (!m_surfMaps)
0251   {
0252     std::cout << PHWHERE
0253               << "ActsSurfaceMaps not found on node tree. Exiting"
0254               << std::endl;
0255     return Fun4AllReturnCodes::ABORTRUN;
0256   }
0257
0258   // The hits are stored in hitsets, where each hitset contains all hits in a given TPC readout (layer, sector, side), so clusters are confined to a hitset
0259   // The TPC clustering is more complicated than for the silicon, because we have to deal with overlapping clusters
0260
0261   // loop over the TPC HitSet objects
0262   TrkrHitSetContainer::ConstRange hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId);
0263   // const int num_hitsets = std::distance(hitsetrange.first,hitsetrange.second);
0264
0265   // Loop over R positions & sectors
0266   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0267        hitsetitr != hitsetrange.second;
0268        ++hitsetitr)
0269   {
0270     TrkrHitSet *hitset = hitsetitr->second;
0271     unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
0272     int side = TpcDefs::getSide(hitsetitr->first);
0273     unsigned int sector = TpcDefs::getSectorId(hitsetitr->first);
0274     PHG4TpcCylinderGeom *layergeom = geom_container->GetLayerCellGeom(layer);
0275
0276     _hit_sector = sector;
0277     _hit_layer = layer;
0278
0279     double radius = layergeom->get_radius();  // in cm
0280
0281     _hit_r = radius;
0282
0283     unsigned short NPhiBins = (unsigned short) layergeom->get_phibins();
0284     unsigned short NPhiBinsSector = NPhiBins / 12;
0285     unsigned short NZBins = (unsigned short) layergeom->get_zbins();
0286     unsigned short NZBinsSide = NZBins / 2;
0287     unsigned short NZBinsMin = 0;
0288     unsigned short PhiOffset = NPhiBinsSector * sector;
0289
0290     if (side == 0)
0291     {
0292       NZBinsMin = 0;
0293       NZBinsMax = NZBins / 2 - 1;
0294     }
0295     else
0296     {
0297       NZBinsMin = NZBins / 2;
0298       NZBinsMax = NZBins;
0299     }
0300     unsigned short ZOffset = NZBinsMin;
0301     // Gives per pad ADC for particular R and sector in event
0302     int perPadADC = 0;
0303     unsigned short phibins = NPhiBinsSector;
0304     unsigned short phioffset = PhiOffset;
0305     unsigned short zbins = NZBinsSide;
0306     unsigned short zoffset = ZOffset;
0307     float sumADC = perPadADC;  // 14 lines later just set to zero
0308
0309     // phibins - number of pads in the sector
0310     // The Sampa clock time is 18.8 MHz, so the sampling time is 53.2 ns = 1 Z bin.
0311     // The Z bin range for one side of the TPC is 0-248.
0312     // Check: 249 bins x 53.2 ns = 13.2 microseconds.
0313     // The maximum drift time in the TPC is 13.2 microseconds.
0314     // So, 53.2 ns / Z bin.
0315
0316     TF1 *f1 = new TF1("f1", "[0]*exp(-(x-[1])/[2])", 0, 1000);
0317     f1->SetParameter(0, 0.005);
0318     f1->SetParameter(1, 0);
0319     f1->SetParameter(2, 60);  // in terms of 50nsec time bins
0320
0321     //    sumADC=0.; // this is set to zero 14 lines up (perPadADC is zero)
0322     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0323
0324     std::vector<unsigned short> adcval(zbins, 0);
0325     //    std::multimap<unsigned short, ihit> all_hit_map;
0326     //    std::vector<ihit> hit_vect;
0327     // Loop over phi & z
0328     for (TrkrHitSet::ConstIterator hitr = hitrangei.first; hitr != hitrangei.second; ++hitr)
0329     {
0330       unsigned short phibin = TpcDefs::getPad(hitr->first) - phioffset;
0331       unsigned short zbin = TpcDefs::getTBin(hitr->first) - zoffset;
0332       float_t fadc = (hitr->second->getAdc());  // proper int rounding +0.5
0333       unsigned short adc = 0;
0334       if (fadc > 0)
0335       {
0336         adc = (unsigned short) fadc;
0337       }
0338       if (phibin >= phibins)
0339       {
0340         continue;
0341       }
0342       if (zbin >= zbins)
0343       {
0344         continue;  // zbin is unsigned int, <0 cannot happen
0345       }
0346       adcval[zbin] = (unsigned short) adc;
0347       sumADC += adc;
0348     }
0349     // Define ion-induced charge
0350     sumADC /= phibins;
0351     float ind_charge = -0.5 * sumADC * _CScale;  // CScale is the coefficient related to the capacitance of the bottom layer of the bottom GEM
0352     double pi = M_PI;
0353
0354     for (TrkrHitSet::ConstIterator hitr = hitrangei.first; hitr != hitrangei.second; ++hitr)
0355     {
0356       unsigned short phibin = TpcDefs::getPad(hitr->first);
0357       unsigned short tbin = TpcDefs::getTBin(hitr->first);
0358       // Get the hitkey
0359       TrkrDefs::hitkey hitkey = TpcDefs::genHitKey(phibin, tbin);
0360       TrkrHit *hit = nullptr;
0361       hit = hitsetitr->second->getHit(hitkey);
0362
0363       tbin = TpcDefs::getTBin(hitr->first);
0364       phibin = TpcDefs::getPad(hitr->first);
0365       double phi_center = layergeom->get_phicenter(phibin, side);
0366       if (phi_center < 0)
0367       {
0368         phi_center += 2 * pi;
0369       }
0370       _hit_phi = phi_center;
0371       _hit_z = AdcClockPeriod * MaxTBins * _drift_velocity / 2.0 - layergeom->get_zcenter(tbin) * _drift_velocity;
0372       if (side == 0)
0373       {
0374         _hit_z *= -1.0;
0375       }
0376
0377       if (hit)
0378       {
0379         _hit_e = hit->getEnergy();
0380       }
0381       _hit_adc = 0;
0382       _hit_adc_bls = 0;
0383
0384       float_t fadc = (hitr->second->getAdc());  // proper int rounding +0.5
0385       _hit_adc = fadc;
0386       _hit_adc_bls = _hit_adc + int(ind_charge);
0387
0388       if (hit && _hit_adc > 0)
0389       {
0390         // Trkr hit has only one value m_adc which is energy and ADC at the same time
0391         if (_hit_adc_bls > 0)
0392         {
0393           hit->setAdc(_hit_adc_bls);
0394         }
0395         else
0396         {
0397           hit->setAdc(0);
0398         }
0399       }
0400       if (_writeTree == 1)
0401       {
0402         _rawHits->Fill();
0403       }
0404     }
0405
0406     // hitsetitr++;
0407   }
0408
0409   // pthread_attr_destroy(&attr);
0410
0411   // wait for completion of all threads
0412   // for( const auto& thread_pair:threads )
0413   //{
0414   //  int rc2 = pthread_join(thread_pair.thread, nullptr);
0415   //  if (rc2)
0416   //  { std::cout << "Error:unable to join," << rc2 << std::endl; }
0417   //}
0418
0419   if (Verbosity() > 0)
0420   {
0421     std::cout << "TPC Clusterizer found " << m_clusterlist->size() << " Clusters " << std::endl;
0422   }
0423   std::cout << "PHG4TpcPadBaselineShift::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0424   return Fun4AllReturnCodes::EVENT_OK;
0425 }
0426
0427 //____________________________________________________________________________..
0428 // int PHG4TpcPadBaselineShift::ResetEvent(PHCompositeNode *topNode)
0429 //{
0430 //  std::cout << "PHG4TpcPadBaselineShift::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0431 //  return Fun4AllReturnCodes::EVENT_OK;
0432 //}
0433
0434 //____________________________________________________________________________..
0435 // int PHG4TpcPadBaselineShift::EndRun(const int runnumber)
0436 //{
0437 //  std::cout << "PHG4TpcPadBaselineShift::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0438 //  return Fun4AllReturnCodes::EVENT_OK;
0439 //}
0440
0441 //____________________________________________________________________________..
0442 int PHG4TpcPadBaselineShift::End(PHCompositeNode * /*topNode*/)
0443 {
0444   if (_writeTree == 1)
0445   {
0446     outfile->cd();
0447     outfile->Write();
0448     outfile->Close();
0449     std::cout << "PHG4TpcPadBaselineShift::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0450   }
0451   // return 0;
0452   return Fun4AllReturnCodes::EVENT_OK;
0453 }
0454
0455 //____________________________________________________________________________..
0456 // int PHG4TpcPadBaselineShift::Reset(PHCompositeNode *topNode)
0457 //{
0458 // std::cout << "PHG4TpcPadBaselineShift::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0459 //  return Fun4AllReturnCodes::EVENT_OK;
0460 //}
0461
0462 //____________________________________________________________________________..
0463 // void PHG4TpcPadBaselineShift::Print(const std::string &what) const
0464 //{
0465 //  std::cout << "PHG4TpcPadBaselineShift::Print(const std::string &what) const Printing info for " << what << std::endl;
0466 //}
0467 void PHG4TpcPadBaselineShift::setScale(float CScale)
0468 {
0469   _CScale = CScale;
0470   std::cout << "PHG4TpcPadBaselineShift::setFileName: Scale factor is set to:" << CScale << std::endl;
0471 }
0472 void PHG4TpcPadBaselineShift::setFileName(const std::string &filename)
0473 {
0474   _filename = filename;
0475   std::cout << "PHG4TpcPadBaselineShift::setFileName: Output file name for PHG4TpcPadBaselineShift is set to:" << filename << std::endl;
0476 }
0477 void PHG4TpcPadBaselineShift::writeTree(int f_writeTree)
0478 {
0479   _writeTree = f_writeTree;
0480   std::cout << "PHG4TpcPadBaselineShift::writeTree: True" << std::endl;
0481 }