Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:15:37 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!