tpc/fillDigitalCurrentMaps/readDigitalCurrents.cc

0001 #include "readDigitalCurrents.h"
0002
0003 #include <g4detectors/PHG4CylinderCellGeom.h>
0004 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
0005 #include <g4detectors/PHG4TpcCylinderGeom.h>
0006 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0007
0008 #include <trackbase/TpcDefs.h>
0009 #include <trackbase/TrkrDefs.h>
0010 #include <trackbase/TrkrHit.h>
0011 #include <trackbase/TrkrHitSet.h>
0012 #include <trackbase/TrkrHitSetContainer.h>
0013
0014 #include <fun4all/Fun4AllHistoManager.h>
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016 #include <fun4all/SubsysReco.h>  // for SubsysReco
0017
0018 #include <phool/getClass.h>
0019 #include <phool/phool.h>  // for PHWHERE
0020 #include <phool/sphenix_constants.h>
0021
0022 #include <TFile.h>
0023 #include <TH1.h>
0024 #include <TH2.h>
0025 #include <TH3.h>
0026
0027 #include <algorithm>  // for max
0028 #include <cassert>
0029 #include <cmath>      // for M_PI, sin, cos
0030 #include <fstream>
0031 #include <iostream>
0032 #include <set>
0033 #include <sstream>
0034 #include <string>
0035 #include <utility>  // for pair
0036
0037 bool IsOverFrame(double r, double phi);
0038
0039 bool IsOverFrame(double r, double phi)
0040 {
0041   // these parameters are taken from Feb 12 drawings of frames.
0042   double tpc_frame_side_gap = 0.8;    // mm //space between radial line and start of frame
0043   double tpc_frame_side_width = 2.6;  // mm //thickness of frame
0044   double tpc_margin = 0.0;            // mm // extra gap between edge of frame and start of GEM holes
0045
0046   double tpc_frame_r3_outer = 758.4;  // mm inner edge of larger-r frame of r3
0047   double tpc_frame_r3_inner = 583.5;  // mm outer edge of smaller-r frame of r3
0048
0049   double tpc_frame_r2_outer = 574.9;  // mm inner edge of larger-r frame of r2
0050   double tpc_frame_r2_inner = 411.4;  // mm outer edge of smaller-r frame of r2
0051
0052   double tpc_frame_r1_outer = 402.6;  // mm inner edge of larger-r frame of r1
0053   double tpc_frame_r1_inner = 221.0;  // mm outer edge of smaller-r frame of r1
0054
0055   // double tpc_sec0_phi=0.0;//get_double_param("tpc_sec0_phi");
0056
0057   // if the coordinate is in the radial spaces of the frames, return true:
0058   if (r < tpc_frame_r1_inner + tpc_margin)
0059   {
0060     return true;
0061   }
0062   if (r > tpc_frame_r1_outer - tpc_margin && r < tpc_frame_r2_inner + tpc_margin)
0063   {
0064     return true;
0065   }
0066   if (r > tpc_frame_r2_outer - tpc_margin && r < tpc_frame_r3_inner + tpc_margin)
0067   {
0068     return true;
0069   }
0070   if (r > tpc_frame_r3_outer - tpc_margin)
0071   {
0072     return true;
0073   }
0074
0075   // if the coordinate is within gap+width of a sector boundary, return true:
0076   // note that this is not a line of constant radius, but a linear distance from a radius.
0077
0078   // find the two spokes we're between:
0079
0080   double sectorangle = (M_PI / 6);
0081   double nsectors = phi / sectorangle;
0082   int nsec = std::floor(nsectors);
0083   double reduced_phi = phi - nsec * sectorangle;  // between zero and sixty degrees.
0084   double dist_to_previous = r * std::sin(reduced_phi);
0085   double dist_to_next = r * std::sin(sectorangle - reduced_phi);
0086   if (dist_to_previous < tpc_frame_side_gap + tpc_frame_side_width + tpc_margin)
0087   {
0088     return true;
0089   }
0090   if (dist_to_next < tpc_frame_side_gap + tpc_frame_side_width + tpc_margin)
0091   {
0092     return true;
0093   }
0094
0095   return false;
0096 }
0097
0098 //____________________________________________________________________________..
0099 readDigitalCurrents::readDigitalCurrents(const std::string &name, const std::string &filename)
0100   : SubsysReco(name)
0101   , hm(nullptr)
0102   , _filename(filename)
0103   , _ampIBFfrac(0.02)
0104   , _collSyst(0)
0105 // , outfile(nullptr)
0106 {
0107   std::cout << "readDigitalCurrents::readDigitalCurrents(const std::string &name) Calling ctor" << std::endl;
0108 }
0109
0110 //____________________________________________________________________________..
0111 readDigitalCurrents::~readDigitalCurrents()
0112 {
0113   std::cout << "readDigitalCurrents::~readDigitalCurrents() Calling dtor" << std::endl;
0114   delete hm;
0115 }
0116
0117 //____________________________________________________________________________..
0118 int readDigitalCurrents::Init(PHCompositeNode * /*topNode*/)
0119 {
0120   std::cout << "readDigitalCurrents::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0121
0122   int nz = 72;
0123   double z_rdo = 108 * cm;
0124
0125   int nr = 159;
0126   // const int nphi=128*3;
0127   // const int nz=62*2;
0128   // double z_rdo=105.5*cm;
0129   // double rmin=20*cm;
0130   double rmax = 78 * cm;
0131
0132   hm = new Fun4AllHistoManager("HITHIST");
0133   const int r_bins_N = 66;  // 51;
0134   double r_bins[r_bins_N + 1] = {217.83,
0135                                  223.83, 229.83, 235.83, 241.83, 247.83, 253.83, 259.83, 265.83, 271.83, 277.83, 283.83, 289.83, 295.83, 301.83, 306.83,
0136                                  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,
0137                                  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,
0138                                  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};
0139   const int nphi = 205;
0140   double phi_bins[nphi + 1] = {0., 6.3083 - 2 * M_PI, 6.3401 - 2 * M_PI, 6.372 - 2 * M_PI, 6.4039 - 2 * M_PI, 6.4358 - 2 * M_PI, 6.4676 - 2 * M_PI, 6.4995 - 2 * M_PI, 6.5314 - 2 * M_PI,
0141                                0.2618, 0.2937, 0.3256, 0.3574, 0.3893, 0.4212, 0.453, 0.4849, 0.5168, 0.5487, 0.5805, 0.6124, 0.6443, 0.6762, 0.7081,
0142                                0.7399, 0.7718, 0.7854, 0.8173, 0.8491, 0.881, 0.9129, 0.9448, 0.9767, 1.0085, 1.0404, 1.0723, 1.1041, 1.136, 1.1679,
0143                                1.1998, 1.2317, 1.2635, 1.2954, 1.309, 1.3409, 1.3727, 1.4046, 1.4365, 1.4684, 1.5002, 1.5321, 1.564, 1.5959, 1.6277,
0144                                1.6596, 1.6915, 1.7234, 1.7552, 1.7871, 1.819, 1.8326, 1.8645, 1.8963, 1.9282, 1.9601, 1.992, 2.0238, 2.0557, 2.0876,
0145                                2.1195, 2.1513, 2.1832, 2.2151, 2.247, 2.2788, 2.3107, 2.3426, 2.3562, 2.3881, 2.42, 2.4518, 2.4837, 2.5156, 2.5474,
0146                                2.5793, 2.6112, 2.6431, 2.6749, 2.7068, 2.7387, 2.7706, 2.8024, 2.8343, 2.8662, 2.8798, 2.9117, 2.9436, 2.9754, 3.0073,
0147                                3.0392, 3.0711, 3.1029, 3.1348, 3.1667, 3.1986, 3.2304, 3.2623, 3.2942, 3.326, 3.3579, 3.3898, 3.4034, 3.4353, 3.4671,
0148                                3.499, 3.5309, 3.5628, 3.5946, 3.6265, 3.6584, 3.6903, 3.7221, 3.754, 3.7859, 3.8178, 3.8496, 3.8815, 3.9134, 3.927,
0149                                3.9589, 3.9907, 4.0226, 4.0545, 4.0864, 4.1182, 4.1501, 4.182, 4.2139, 4.2457, 4.2776, 4.3095, 4.3414, 4.3732, 4.4051,
0150                                4.437, 4.4506, 4.4825, 4.5143, 4.5462, 4.5781, 4.61, 4.6418, 4.6737, 4.7056, 4.7375, 4.7693, 4.8012, 4.8331, 4.865,
0151                                4.8968, 4.9287, 4.9606, 4.9742, 5.0061, 5.0379, 5.0698, 5.1017, 5.1336, 5.1654, 5.1973, 5.2292, 5.2611, 5.2929, 5.3248,
0152                                5.3567, 5.3886, 5.4204, 5.4523, 5.4842, 5.4978, 5.5297, 5.5615, 5.5934, 5.6253, 5.6572, 5.689, 5.7209, 5.7528, 5.7847,
0153                                5.8165, 5.8484, 5.8803, 5.9122, 5.944, 5.9759, 6.0078, 6.0214, 6.0533, 6.0851, 6.117, 6.1489, 6.1808, 6.2127, 6.2445,
0154                                6.2764, 2 * M_PI};
0155
0156   // R0 sector 0: 2.36679 < phi < 2.86919
0157   // R0 sector 1: 1.8432 < phi < 2.3456
0158   // R0 sector 2: 1.3196 < phi < 1.822
0159   // R0 sector 3: 0.795998 < phi < 1.2984
0160   // R0 sector 4: 0.272399 < phi < 0.774799
0161   // R0 sector 5: -0.2512 < phi < 0.2512
0162   // R0 sector 6: -0.774799 < phi < -0.272399
0163   // R0 sector 7: -1.2984 < phi < -0.795998
0164   // R0 sector 8: -1.822 < phi < -1.3196
0165   // R0 sector 9: -2.3456 < phi < -1.8432
0166   // R0 sector 10: -2.86919 < phi < -2.36679
0167   // R0 sector 11: -3.39279 < phi < -2.89039
0168
0169   double *z_bins = new double[2 * nz + 1];
0170   for (int z = 0; z <= 2 * nz; z++)
0171   {
0172     z_bins[z] = -z_rdo + z_rdo / nz * z;
0173   }
0174
0175   _h_R = new TH1F("_h_R", "_h_R;R, [m]", r_bins_N, r_bins);
0176   _h_hits = new TH1F("_h_hits", "_h_hits;N, [hit]", 1e5, 0 - 0.5, 1e5 - 0.5);
0177   _h_hit_XY = new TH2F("_h_hit_XY", "_h_hit_XY;X, [m];Y, [m]", 4 * nr, -1 * rmax, rmax, 4 * nr, -1 * rmax, rmax);
0178
0179   //_h_SC_ibf  = new TH3F("_h_SC_ibf" ,"_h_SC_ibf;#phi, [rad];R, [mm];Z, [mm]"   ,nphi,phi_bins,r_bins_N ,r_bins,2*nz,z_bins);
0180   _h_DC_SC = new TH3F("_h_DC_SC", "_h_DC_SC;#phi, [rad];R, [mm];Z, [mm]", nphi, phi_bins, r_bins_N, r_bins, 2 * nz, z_bins);
0181   _h_DC_SC_XY = new TH2F("_h_DC_SC_XY", "_h_DC_SC_XY;X, [mm];Y, [mm];ADC;", 4 * nr, -1 * rmax, rmax, 4 * nr, -1 * rmax, rmax);
0182   _h_DC_E = new TH2F("_h_DC_E", "_h_DC_E;ADC;E", 200, -100, 2e3 - 100, 500, -100, 5e3 - 100);
0183
0184   std::string name;
0185   std::string name_ax;
0186   for (int iz = 0; iz < nFrames; iz++)
0187   {
0188     name = "_h_SC_ibf_" + std::to_string(iz);
0189     name_ax = "_h_SC_ibf_" + std::to_string(iz) + ";#phi, [rad];R, [mm];Z, [mm]";
0190     _h_SC_ibf[iz] = new TH3F(name.c_str(), name_ax.c_str(), nphi, phi_bins, r_bins_N, r_bins, 2 * nz, z_bins);
0191
0192     hm->registerHisto(_h_SC_ibf[iz]);
0193   }
0194
0195   hm->registerHisto(_h_R);
0196   hm->registerHisto(_h_hits);
0197   hm->registerHisto(_h_DC_SC);
0198   hm->registerHisto(_h_DC_SC_XY);
0199   hm->registerHisto(_h_hit_XY);
0200   hm->registerHisto(_h_DC_E);
0201
0202   _event_timestamp = 0;
0203
0204   if (_fillCSVFile)
0205   {
0206     myCSVFile.open("./Files/example_1ms_120evts_AA.csv");
0207     myCSVFile << "Event,"
0208               << "T,"
0209               << "Pad,"
0210               << "Radius,"
0211               << "ADC"
0212               << "\n";
0213   }
0214
0215   return Fun4AllReturnCodes::EVENT_OK;
0216 }
0217
0218 //____________________________________________________________________________..
0219 int readDigitalCurrents::InitRun(PHCompositeNode * /*topNode*/)
0220 {
0221   std::cout << "readDigitalCurrents::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0222   std::string line;
0223   // AA collisions timestamps
0224   // std::string txt_file = "/sphenix/user/shulga/Work/IBF/DistortionMap/timestamps_50kHz.txt";
0225   std::string txt_file = "/sphenix/user/shulga/Work/TpcPadPlane_phi_coresoftware/coresoftware/calibrations/tpc/fillSpaceChargeMaps/data/timestamps_50kHz_1M.txt";
0226   int start_line = 3;
0227   if (_collSyst == 1)
0228   {
0229     // pp collisions timestamps
0230     txt_file = "/phenix/u/hpereira/sphenix/work/g4simulations/timestamps_3MHz.txt";
0231     // txt_file = "/sphenix/user/shulga/Work/IBF/DistortionMap/timestamps_50kHz.txt";
0232     start_line = 2;
0233   }
0234   std::ifstream InputFile(txt_file);
0235   if (InputFile.is_open())
0236   {
0237     int n_line = 0;
0238     while (getline(InputFile, line))
0239     {
0240       n_line++;
0241       if (n_line > start_line)
0242       {
0243         std::istringstream is(line);
0244         double n[2] = {0, 0};
0245         int i = 0;
0246         while (is >> n[i])
0247         {
0248           i++;
0249         }
0250         _timestamps[n[0]] = n[1];
0251         if (n_line < 10)
0252         {
0253           std::cout << n[1] << std::endl;
0254         }
0255         _keys.push_back(int(n[0]));
0256       }
0257     }
0258     InputFile.close();
0259   }
0260
0261   else
0262   {
0263     std::cout << "Unable to open file:" << txt_file << std::endl;
0264   }
0265
0266   TFile *MapsFile;
0267   // if(_fUseIBFMap){
0268   MapsFile = new TFile("/sphenix/user/shulga/Work/IBF/DistortionMap/IBF_Map.root", "READ");
0269   if (MapsFile->IsOpen())
0270   {
0271     std::cout << "Gain/IBF Maps File opened successfully" << std::endl;
0272   }
0273   //_h_modules_anode       = (TH2F*)MapsFile ->Get("h_modules_anode")      ->Clone("_h_modules_anode");
0274   _h_modules_measuredibf = (TH2F *) MapsFile->Get("h_modules_measuredibf")->Clone("_h_modules_measuredibf");
0275   //}
0276   return Fun4AllReturnCodes::EVENT_OK;
0277 }
0278
0279 //____________________________________________________________________________..
0280 int readDigitalCurrents::process_event(PHCompositeNode *topNode)
0281 {
0282   // double bX = _beamxing;
0283   // float bX = 1508071;
0284   // double z_bias_avg = 0;
0285   // if (_fAvg==1){
0286   //   z_bias_avg=1.05*(float) rand()/RAND_MAX;
0287   // }
0288   int bemxingsInFile = _keys.size();
0289   if (_evtstart >= bemxingsInFile)
0290   {
0291     _evtstart = _evtstart - bemxingsInFile;
0292   }
0293   int key = _keys.at(_evtstart);
0294   _event_timestamp = (float) _timestamps[key] * ns;  // units in seconds
0295   _event_bunchXing = key;
0296   if (_evtstart % 100 == 0)
0297   {
0298     std::cout << "_evtstart = " << _evtstart << std::endl;
0299   }
0300   _evtstart++;
0301
0302   // std::cout << "readDigitalCurrents::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0303   std::set<std::string>::const_iterator iter;
0304   // //nodename << "G4HIT_TPC";
0305   std::string nodename = "TRKR_HITSET";
0306
0307   // //  SvtxEvaluator
0308   // SvtxEvaluator *hits = findNode::getClass<SvtxEvaluator>(topNode, nodename.str().c_str());
0309   // //int n_hits = 0;
0310   // if (hits){
0311   //   PHG4HitContainer::ConstRange hit_range = hits->getHits();
0312   // }
0313   //===================================
0314   // get node containing the digitized hits
0315   TrkrHitSetContainer *_hitmap = findNode::getClass<TrkrHitSetContainer>(topNode, nodename);
0316   if (!_hitmap)
0317   {
0318     std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
0319     return Fun4AllReturnCodes::ABORTRUN;
0320   }
0321   std::string geo_nodename = "CYLINDERCELLGEOM_SVTX";
0322
0323   PHG4TpcCylinderGeomContainer *_geom_container_ccgc = nullptr;
0324   PHG4CylinderCellGeomContainer *_geom_container_cgc = nullptr;
0325   if (_f_ccgc == 1)
0326   {
0327     _geom_container_ccgc = findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, geo_nodename);
0328     if (!_geom_container_ccgc)
0329     {
0330       std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0331       return Fun4AllReturnCodes::ABORTRUN;
0332     }
0333   }
0334   else
0335   {
0336     _geom_container_cgc = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, geo_nodename);
0337
0338     if (!_geom_container_cgc)
0339     {
0340       std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0341       return Fun4AllReturnCodes::ABORTRUN;
0342     }
0343   }
0344
0345   // loop over all the hits
0346   // hits are stored in hitsets, so have to get the hitset first
0347   int n_hits = 0;
0348   // float _event_bunchXing = 1508071;
0349   TrkrHitSetContainer::ConstRange all_hitsets = _hitmap->getHitSets();
0350   for (TrkrHitSetContainer::ConstIterator iter_hitset = all_hitsets.first; iter_hitset != all_hitsets.second; ++iter_hitset)
0351   {
0352     // checking that the object is inside TPC
0353     if (TrkrDefs::getTrkrId(iter_hitset->first) == TrkrDefs::tpcId)
0354     {
0355       TrkrDefs::hitsetkey hitsetkey = iter_hitset->first;
0356       const unsigned int zside = TpcDefs::getSide(hitsetkey);
0357       TrkrHitSet::ConstRange range = iter_hitset->second->getHits();
0358       unsigned int layer = TrkrDefs::getLayer(iter_hitset->first);
0359       // if(layer>6){
0360       PHG4TpcCylinderGeom *layergeom_ccgc = nullptr;
0361       PHG4CylinderCellGeom *layergeom_cgc = nullptr;
0362       double radius = 0;
0363       if (_f_ccgc == 1)
0364       {
0365         layergeom_ccgc = _geom_container_ccgc->GetLayerCellGeom(layer);
0366         radius = layergeom_ccgc->get_radius() * cm;
0367       }
0368       else
0369       {
0370         layergeom_cgc = _geom_container_cgc->GetLayerCellGeom(layer);
0371         radius = layergeom_cgc->get_radius() * cm;
0372       }
0373       // PHG4TpcCylinderGeom *layergeom = _geom_container->GetLayerCellGeom(layer);
0374       // double radius = layergeom->get_radius()*cm;  // returns center of the layer
0375       int min_phiBin = 1e5;
0376       int max_phiBin = -1;
0377       for (TrkrHitSet::ConstIterator hit_iter = range.first; hit_iter != range.second; ++hit_iter)
0378       {
0379         // int f_fill_ibf=0;
0380         int f_fill_ibf[30] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
0381
0382         unsigned short phibin = TpcDefs::getPad(hit_iter->first);
0383         unsigned short zbin = TpcDefs::getTBin(hit_iter->first);
0384         double _drift_velocity = 8.0e-3;  // ActsGeometry::get_drift_velocity();
0385         double phi_center = 0;
0386         if (_f_ccgc == 1)
0387         {
0388           phi_center = layergeom_ccgc->get_phicenter(phibin, zside);
0389         }
0390         else
0391         {
0392           phi_center = layergeom_cgc->get_phicenter(phibin);
0393         }
0394         if (phi_center < 0)
0395         {
0396           phi_center += 2 * M_PI;
0397         }
0398         if (phi_center < M_PI / 2 + M_PI / 12 && phi_center > M_PI / 2 - M_PI / 12)
0399         {
0400           min_phiBin = std::min<int>(min_phiBin, phibin);
0401           max_phiBin = std::max<int>(max_phiBin, phibin);
0402         }
0403         float x = radius * cos(phi_center);
0404         float y = radius * sin(phi_center);
0405
0406         double z = 0;  // layergeom->get_zcenter(zbin)*cm;
0407         if (_f_ccgc == 1)
0408         {
0409           z = layergeom_ccgc->get_zcenter(zbin) * _drift_velocity * cm;  //*cm/ns;
0410           if (zside == 0)
0411           {
0412             z = -z;
0413           }
0414         }
0415         else
0416         {
0417           z = layergeom_cgc->get_zcenter(zbin) * cm;
0418         }
0419         TrkrHit *hit = hit_iter->second;
0420         int adc = hit->getAdc() - adc_pedestal;
0421         float E = hit->getEnergy();
0422         // double z = 0;
0423         // double z_prim = -1*1e10;
0424         // double z_ibf =  -1*1e10;
0425         double z_ibf[30] = {-1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10,
0426                             -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10,
0427                             -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10};
0428
0429         int RBin = _h_R->GetXaxis()->FindBin(radius);
0430         if ((RBin > 33 && RBin < 50) && z > 0)
0431         {
0432       assert(layergeom_ccgc);
0433           int nRBins = layergeom_ccgc->get_phibins();
0434           if (phibin < nRBins / 12)
0435           {
0436             if (_fillCSVFile)
0437             {
0438               myCSVFile << _evtstart << ","
0439                         << layergeom_ccgc->get_zcenter(zbin) << ","
0440                         << phibin << ","
0441                         << RBin - 34 << ","
0442                         << adc << "\n";
0443             }
0444             _h_hit_XY->Fill(x, y);
0445           }
0446         }
0447         // if(!IsOverFrame(radius/mm,phi_center)){
0448         for (int iz = 0; iz < nFrames; iz++)
0449         {
0450           double bX = _beamxing[iz];
0451           if (_event_bunchXing <= bX)
0452           {
0453             if (z >= 0 && z < 1.055 * m)
0454             {
0455               z_ibf[iz] = 1.055 * m - (bX - _event_bunchXing) * sphenix_constants::time_between_crossings * vIon * ns;
0456               if (z_ibf[iz] > 0 && z_ibf[iz] < 1.055 * m)
0457               {
0458                 f_fill_ibf[iz] = 1;
0459               }
0460             }
0461             if (z < 0 && z > -1.055 * m)
0462             {
0463               z_ibf[iz] = -1.055 * m + (bX - _event_bunchXing) * sphenix_constants::time_between_crossings * vIon * ns;
0464               if (z_ibf[iz] < 0 && z_ibf[iz] > -1.055 * m)
0465               {
0466                 f_fill_ibf[iz] = 1;
0467               }
0468             }
0469           }
0470         }
0471         if (z >= 0 && z < 1.055 * m)
0472         {
0473           if (adc >= 0)
0474           {
0475             n_hits++;
0476           }
0477           if (adc >= 0)
0478           {
0479             _h_DC_E->Fill(adc, E);
0480           }
0481         }
0482         if (z < 0 && z > -1.055 * m)
0483         {
0484           if (adc >= 0)
0485           {
0486             n_hits++;
0487           }
0488           if (adc >= 0)
0489           {
0490             _h_DC_E->Fill(adc, E);
0491           }
0492         }
0493
0494         // Reading IBF and Gain weights according to X-Y position
0495         float w_ibf = 1.;
0496         // float w_gain = 1.;
0497         // if(_fUseIBFMap){
0498         int bin_x = _h_modules_measuredibf->GetXaxis()->FindBin(x / mm);
0499         int bin_y = _h_modules_measuredibf->GetYaxis()->FindBin(y / mm);
0500         w_ibf = _h_modules_measuredibf->GetBinContent(bin_x, bin_y);
0501         // w_gain = _h_modules_anode->GetBinContent(bin_x,bin_y);
0502         w_ibf = 1.;
0503         //}
0504         float w_adc = adc * w_ibf;
0505         _h_DC_SC->Fill(phi_center, radius, z, w_adc);
0506         _h_DC_SC_XY->Fill(x, y, w_adc);
0507         if (f_fill_ibf[0] == 1)
0508         {
0509           _h_R->Fill(radius);
0510         }
0511         for (int iz = 0; iz < nFrames; iz++)
0512         {
0513           if (f_fill_ibf[iz] == 1)
0514           {
0515             _h_SC_ibf[iz]->Fill(phi_center, radius, z_ibf[iz], w_adc);
0516           }
0517         }
0518         //}
0519         // if(n_hits%100==0) std::cout<<radius<<"|"<<phi_center<<"|"<<z<<std::endl;
0520       }
0521       // std::cout<<" min_phiBin"<< min_phiBin <<" max_phiBin"<< max_phiBin<< std::endl;
0522     }
0523   }
0524
0525   // TrkrHitSetContainer, TrkrHitSet, TrkrHit, and TrkrDefs objects used above in offline/packages/trackbase.
0526   _h_hits->Fill(n_hits);
0527
0528   return Fun4AllReturnCodes::EVENT_OK;
0529 }
0530
0531 //____________________________________________________________________________..
0532 int readDigitalCurrents::ResetEvent(PHCompositeNode * /*topNode*/)
0533 {
0534   // std::cout << "readDigitalCurrents::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0535   return Fun4AllReturnCodes::EVENT_OK;
0536 }
0537
0538 //____________________________________________________________________________..
0539 int readDigitalCurrents::EndRun(const int runnumber)
0540 {
0541   std::cout << "readDigitalCurrents::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0542   return Fun4AllReturnCodes::EVENT_OK;
0543 }
0544
0545 //____________________________________________________________________________..
0546 int readDigitalCurrents::End(PHCompositeNode * /*topNode*/)
0547 {
0548   std::cout << "readDigitalCurrents::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0549   _h_R->Sumw2(false);
0550   _h_hits->Sumw2(false);
0551   _h_DC_E->Sumw2(false);
0552   _h_DC_SC->Sumw2(false);
0553   _h_hit_XY->Sumw2(false);
0554   _h_DC_SC_XY->Sumw2(false);
0555   for (auto &iz : _h_SC_ibf)
0556   {
0557     iz->Sumw2(false);
0558   }
0559   hm->dumpHistos(_filename, "RECREATE");
0560   if (_fillCSVFile)
0561   {
0562     myCSVFile.close();
0563   }
0564   return Fun4AllReturnCodes::EVENT_OK;
0565 }
0566
0567 //____________________________________________________________________________..
0568 int readDigitalCurrents::Reset(PHCompositeNode * /*topNode*/)
0569 {
0570   std::cout << "readDigitalCurrents::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0571   return Fun4AllReturnCodes::EVENT_OK;
0572 }
0573
0574 //____________________________________________________________________________..
0575 void readDigitalCurrents::Print(const std::string &what) const
0576 {
0577   std::cout << "readDigitalCurrents::Print(const std::string &what) const Printing info for " << what << std::endl;
0578 }
0579
0580 void readDigitalCurrents::SetEvtStart(int newEvtStart)
0581 {
0582   _evtstart = newEvtStart;
0583   std::cout << "Start event is set to: " << newEvtStart << std::endl;
0584 }
0585 void readDigitalCurrents::FillCSV(int fillCSVFile)
0586 {
0587   _fillCSVFile = fillCSVFile;
0588   std::cout << "Fill CSV file is set to: " << fillCSVFile << std::endl;
0589 }
0590
0591 void readDigitalCurrents::SetBeamXing(const std::vector<int> &beamXs)
0592 {
0593   _beamxing = beamXs;
0594   std::cout << "Initial BeamXing is set to: " << _beamxing[0] << std::endl;
0595 }
0596 void readDigitalCurrents::SetCollSyst(int coll_syst)
0597 {
0598   _collSyst = coll_syst;
0599   std::string s_syst[2] = {"AA", "pp"};
0600   std::cout << "Collision system is set to: " << s_syst[_collSyst] << std::endl;
0601 }
0602
0603 void readDigitalCurrents::SetIBF(double ampIBFfrac)
0604 {
0605   _ampIBFfrac = ampIBFfrac;
0606   std::cout << "IBF is set to: " << _ampIBFfrac << std::endl;
0607 }
0608
0609 void readDigitalCurrents::SetCCGC(double f_ccgc)
0610 {
0611   _f_ccgc = f_ccgc;
0612   std::cout << "IBF is set to: " << _f_ccgc << std::endl;
0613 }