Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 #include "fillSpaceChargeMaps.h"
0002 
0003 #include "Shifter.h"
0004 
0005 #include <g4main/PHG4Hit.h>
0006 #include <g4main/PHG4HitContainer.h>
0007 
0008 #include <fun4all/Fun4AllHistoManager.h>
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010 #include <fun4all/SubsysReco.h>  // for SubsysReco
0011 
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/getClass.h>
0014 #include <phool/sphenix_constants.h>
0015 
0016 #include <TAxis.h>  // for TAxis
0017 #include <TFile.h>
0018 #include <TH1.h>
0019 #include <TH2.h>
0020 #include <TH3.h>
0021 #include <TTree.h>
0022 #include <TVector3.h>
0023 
0024 #include <boost/format.hpp>
0025 
0026 #include <algorithm>  // for max
0027 #include <cstdlib>
0028 #include <fstream>
0029 #include <iostream>
0030 #include <map>
0031 #include <sstream>
0032 #include <string>
0033 #include <utility>  // for pair
0034 #include <vector>
0035 
0036 //____________________________________________________________________________..
0037 fillSpaceChargeMaps::fillSpaceChargeMaps(const std::string &name, const std::string &filename)
0038   : SubsysReco(name)
0039   , _filename(filename)
0040 {
0041   std::cout << "fillSpaceChargeMaps::fillSpaceChargeMaps(const std::string &name) Calling ctor" << std::endl;
0042 }
0043 
0044 //____________________________________________________________________________..
0045 fillSpaceChargeMaps::~fillSpaceChargeMaps()
0046 {
0047   // std::cout << "fillSpaceChargeMaps::~fillSpaceChargeMaps() Calling dtor" << std::endl;
0048   //  delete whatever is created
0049   delete hm;
0050 }
0051 
0052 //____________________________________________________________________________..
0053 int fillSpaceChargeMaps::Init(PHCompositeNode * /*topNode*/)
0054 {
0055   int nz = 72;
0056   double z_rdo = 108 * cm;
0057   outfile = new TFile(_filename.c_str(), "RECREATE");
0058 
0059   hm = new Fun4AllHistoManager("HITHIST");
0060   const int r_bins_N = 66;  // 51;
0061   double r_bins[r_bins_N + 1] = {217.83,
0062                                  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,
0063                                  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,
0064                                  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,
0065                                  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};
0066   double r_bins_edges[r_bins_N + 3] = {217.83 - 2, 217.83,
0067                                        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,
0068                                        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,
0069                                        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,
0070                                        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, 759.11 + 2};
0071 
0072   const int nphi = 205;
0073 
0074   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,
0075                                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,
0076                                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,
0077                                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,
0078                                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,
0079                                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,
0080                                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,
0081                                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,
0082                                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,
0083                                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,
0084                                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,
0085                                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,
0086                                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,
0087                                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,
0088                                6.2764, 2 * M_PI};
0089 
0090   double *z_bins = new double[2 * nz + 1];
0091   for (int z = 0; z <= 2 * nz; z++)
0092   {
0093     z_bins[z] = -z_rdo + z_rdo / nz * z;
0094   }
0095 
0096   _h_R = new TH1F("_h_R", "_h_R;R, [mm]", r_bins_N + 2, r_bins_edges);
0097   _h_hits = new TH1F("_h_hits", "_h_hits;N, [hit]", 1000, 0, 1e6);
0098   _h_DC_E = new TH2F("_h_DC_E", "_h_DC_E;SC;E#times10^{6}", 2000, -100, 2e5 - 100, 1000, -50, 1e3 - 50);
0099   _h_SC_XY = new TH2F("_h_SC_XY", "_h_SC_XY;X, [mm];Y, [mm];ADC;", 4 * 159, -1 * 78 * cm, 78 * cm, 4 * 159, -1 * 78 * cm, 78 * cm);
0100 
0101   for (int iz = 0; iz < 30; iz++)
0102   {
0103     std::string name = (boost::format("_h_SC_ibf_%d") % iz).str();
0104     std::string name_ax = (boost::format("_h_SC_ibf_%d;#phi, [rad];R, [mm];Z, [mm]") % iz).str();
0105     _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);
0106     name = (boost::format("_h_SC_prim_%d") % iz).str();
0107     name_ax = (boost::format("_h_SC_prim_%d;#phi, [rad];R, [mm];Z, [mm]") % iz).str();
0108     _h_SC_prim[iz] = new TH3F(name.c_str(), name_ax.c_str(), nphi, phi_bins, r_bins_N, r_bins, 2 * nz, z_bins);
0109 
0110     hm->registerHisto(_h_SC_prim[iz]);
0111     hm->registerHisto(_h_SC_ibf[iz]);
0112   }
0113   hm->registerHisto(_h_hits);
0114   hm->registerHisto(_h_R);
0115   hm->registerHisto(_h_DC_E);
0116   hm->registerHisto(_h_SC_XY);
0117   //_event_timestamp = 0;
0118   _hit_eion = 0;
0119   _hit_r = 0;
0120   _hit_phi = 0;
0121   _hit_z = 0;
0122   _ibf_vol = 0;
0123   _amp_ele_vol = 0;
0124   if (_fSliming == 1)
0125   {
0126     _rawHits = new TTree("hTree", "tpc hit tree for ionization");
0127     _rawHits->Branch("isOnPlane", &_isOnPlane);
0128     _rawHits->Branch("hit_z", &_hit_z);
0129     _rawHits->Branch("hit_r", &_hit_r);
0130     _rawHits->Branch("hit_phi", &_hit_phi);
0131     _rawHits->Branch("hit_eion", &_hit_eion);
0132     _rawHits->Branch("ibf_vol", &_ibf_vol);
0133     _rawHits->Branch("amp_ele_vol", &_amp_ele_vol);
0134 
0135     _rawHits->Branch("event_timestamp", &_event_timestamp);
0136     _rawHits->Branch("event_bunchXing", &_event_bunchXing);
0137   }
0138   // padplane = new PHG4TpcPadPlaneReadout;
0139   // seggeo = new PHG4TpcCylinderGeomContainer();
0140   return 0;
0141 }
0142 
0143 //____________________________________________________________________________..
0144 int fillSpaceChargeMaps::InitRun(PHCompositeNode * /*topNode*/)
0145 {
0146   // std::cout << "fillSpaceChargeMaps::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0147   std::string line;
0148   // AA collisions timestamps
0149   std::string txt_file = "./data/timestamps_50kHz_1M.txt";
0150   int start_line = 3;
0151   if (_collSyst == 1)
0152   {
0153     // pp collisions timestamps
0154     txt_file = "/phenix/u/hpereira/sphenix/work/g4simulations/timestamps_3MHz.txt";
0155     start_line = 2;
0156   }
0157   std::ifstream InputFile(txt_file);
0158   if (InputFile.is_open())
0159   {
0160     int n_line = 0;
0161     while (getline(InputFile, line))
0162     {
0163       n_line++;
0164       if (n_line > start_line)
0165       {
0166         std::istringstream is(line);
0167         double n[2] = {0, 0};
0168         int i = 0;
0169         while (is >> n[i])
0170         {
0171           i++;
0172         }
0173         _timestamps[n[0]] = n[1];
0174         if (n_line < 10)
0175         {
0176           std::cout << n[1] << std::endl;
0177         }
0178         _keys.push_back(int(n[0]));
0179       }
0180     }
0181     InputFile.close();
0182   }
0183 
0184   else
0185   {
0186     std::cout << "Unable to open file:" << txt_file << std::endl;
0187   }
0188 
0189   if (_fUseIBFMap)
0190   {
0191     TFile *MapsFile = new TFile("./data/IBF_Map.root", "READ");
0192     if (MapsFile->IsOpen())
0193     {
0194       std::cout << "File opened successfully" << std::endl;
0195     }
0196     _h_modules_anode = (TH2F *) MapsFile->Get("h_modules_anode")->Clone("_h_modules_anode");
0197     _h_modules_measuredibf = (TH2F *) MapsFile->Get("h_modules_measuredibf")->Clone("_h_modules_measuredibf");
0198   }
0199 
0200   _mbRate = _freqKhz * kHz;
0201   _xingRate = 9.383 * MHz;
0202   _mean = mbRate / xingRate;
0203   // padplane->CreateReadoutGeometry( PHCompositeNode *, seggeo);
0204   return 0;
0205 }
0206 
0207 //____________________________________________________________________________..
0208 int fillSpaceChargeMaps::process_event(PHCompositeNode *topNode)
0209 {
0210   // double tpc_frame_r3_outer = 759.11*mm;  //758.4;//mm inner edge of larger-r frame of r3
0211   // double tpc_frame_r3_inner = 583.67*mm;  //583.5;//mm outer edge of smaller-r frame of r3
0212 
0213   // double tpc_frame_r2_outer = 574.76*mm;  //574.9;//mm inner edge of larger-r frame of r3
0214   // double tpc_frame_r2_inner = 411.53*mm;  //411.4;//mm outer edge of smaller-r frame of r3
0215 
0216   // double tpc_frame_r1_outer = 402.49*mm;  //402.6;//mm inner edge of larger-r frame of r3
0217   // double tpc_frame_r1_inner = 217.83*mm;  //221.0;//mm outer edge of smaller-r frame of r3
0218 
0219   double z_bias_avg = 0;
0220   int bemxingsInFile = _keys.size();
0221   int key = -1;
0222   if (_fAvg == 1)
0223   {
0224     z_bias_avg = 1.055 * m * (float) rand() / RAND_MAX;
0225   }
0226   if (_evtstart >= bemxingsInFile)
0227   {
0228     _evtstart = _evtstart - bemxingsInFile;
0229   }
0230   key = _keys.at(_evtstart);
0231   _event_timestamp = (float) _timestamps[key] * ns;  // units in seconds
0232   _event_bunchXing = key;
0233 
0234   if (_evtstart % 10 == 0)
0235   {
0236     std::cout << "_evtstart = " << _evtstart << std::endl;
0237   }
0238   _evtstart++;
0239 
0240   Shifter shifter("/sphenix/user/rcorliss/distortion_maps/2021.04/apr07.average.real_B1.4_E-400.0.ross_phi1_sphenix_phislice_lookup_r26xp40xz40.distortion_map.hist.root");
0241 
0242   PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
0243   int n_hits = 0;
0244   if (hits)
0245   {
0246     PHG4HitContainer::ConstRange hit_range = hits->getHits();
0247     for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0248     {
0249       int f_fill_prim[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};
0250       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};
0251 
0252       float hit_x0 = hit_iter->second->get_x(0);
0253       float hit_y0 = hit_iter->second->get_y(0);
0254       float hit_z0 = hit_iter->second->get_z(0);
0255       float hit_x1 = hit_iter->second->get_x(1);
0256       float hit_y1 = hit_iter->second->get_y(1);
0257       float hit_z1 = hit_iter->second->get_z(1);
0258 
0259       float hit_eion = hit_iter->second->get_eion();
0260       float N_electrons = hit_eion * Tpc_ElectronsPerGeV;
0261       float x = (hit_x0 + f * (hit_x1 - hit_x0)) * cm;
0262       float y = (hit_y0 + f * (hit_y1 - hit_y0)) * cm;
0263       float z = (hit_z0 + f * (hit_z1 - hit_z0)) * cm;
0264 
0265       float r = sqrt(x * x + y * y);
0266       float phi = atan2(x, y);
0267       if (phi < 0)
0268       {
0269         phi += 2 * M_PI;
0270       }
0271 
0272       // Shift electrons according to the field maps
0273       TVector3 oldPos(x / cm, y / cm, z / cm);
0274       TVector3 newPos;
0275       if (_shiftElectrons == 1)
0276       {
0277         if (oldPos.z() < 0)
0278         {
0279           oldPos.SetZ(std::abs(oldPos.z()));
0280           newPos = shifter.ShiftForward(oldPos);
0281           newPos.SetZ(newPos.z() * -1);
0282         }
0283         else
0284         {
0285           newPos = shifter.ShiftForward(oldPos);
0286         }
0287       }
0288       else
0289       {
0290         newPos.SetZ(oldPos.Z());
0291         newPos.SetX(oldPos.X());
0292         newPos.SetY(oldPos.Y());
0293       }
0294 
0295       // Reading IBF and Gain weights according to X-Y position
0296       float w_ibf = 1.;
0297       float w_gain = 1.;
0298 
0299       if (_fUseIBFMap)
0300       {
0301         int bin_x = _h_modules_anode->GetXaxis()->FindBin(x / mm);
0302         int bin_y = _h_modules_anode->GetYaxis()->FindBin(y / mm);
0303         w_ibf = _h_modules_measuredibf->GetBinContent(bin_x, bin_y);
0304         w_gain = _h_modules_anode->GetBinContent(bin_x, bin_y);
0305       }
0306       float ionsPerEle = w_gain * _ampGain * w_ibf * _ampIBFfrac;
0307 
0308       // Check that it is on the frame
0309       _isOnPlane = 0;
0310       double dr_bin = -1;
0311       double dphi_bin = -1;
0312       double new_phi = newPos.Phi();
0313       double new_r = newPos.Perp() * cm;
0314       if (new_phi < 0)
0315       {
0316         new_phi += 2 * M_PI;
0317       }
0318       if (!IsOverFrame(new_r, new_phi))
0319       {
0320         _isOnPlane = 1;
0321       }
0322       else
0323       {
0324         std::vector<double> r_phi_bin = putOnPlane(new_r / mm, new_phi);
0325         dr_bin = r_phi_bin[0];
0326         dphi_bin = r_phi_bin[1];
0327       }
0328       _hit_z = z;
0329       _hit_r = r;
0330       _hit_phi = phi;
0331       _hit_eion = hit_eion;
0332       _ibf_vol = N_electrons * ionsPerEle;
0333       _amp_ele_vol = w_gain * _ampGain;
0334       if (new_r < 210 || new_r > 770)
0335       {
0336         continue;
0337       }
0338       if (_fSliming == 1)
0339       {
0340         _rawHits->Fill();
0341       }
0342       double z_prim[30] = {-1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10,
0343                            -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10,
0344                            -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10};
0345       double z_ibf[30] = {-1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10,
0346                           -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10,
0347                           -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10};
0348 
0349       if (_hit_z >= 5 * mm && _hit_z < 1.055 * m)
0350       {
0351         n_hits++;
0352         _h_DC_E->Fill(_ibf_vol, hit_eion * 1e6);
0353         for (int iz = 0; iz < 30; iz++)
0354         {
0355           double bX = _beamxing[iz];
0356           // double bX_end = _beamxing_end[iz];
0357           if (_event_bunchXing <= bX)
0358           {
0359             if (_fAvg == 1)
0360             {
0361               z_prim[iz] = _hit_z - z_bias_avg;
0362               z_ibf[iz] = 1.055 * m - z_bias_avg;
0363             }
0364             else
0365             {
0366               z_prim[iz] = _hit_z - (bX - _event_bunchXing) * sphenix_constants::time_between_crossings * vIon * ns;
0367               z_ibf[iz] = 1.055 * m - (bX - _event_bunchXing) * sphenix_constants::time_between_crossings * vIon * ns;
0368             }
0369             if (z_prim[iz] > 0 && z_prim[iz] < 1.055 * m)
0370             {
0371               f_fill_prim[iz] = 1;
0372             }
0373             if (z_ibf[iz] > 0 && z_ibf[iz] < 1.055 * m)
0374             {
0375               f_fill_ibf[iz] = 1;
0376             }
0377           }
0378         }
0379       }
0380 
0381       if (_hit_z < -5 * mm && _hit_z > -1.055 * m)
0382       {
0383         n_hits++;
0384         _h_DC_E->Fill(_ibf_vol, hit_eion * 1e6);
0385         for (int iz = 0; iz < 30; iz++)
0386         {
0387           double bX = _beamxing[iz];
0388           // double bX_end = _beamxing_end[iz];
0389           if (_event_bunchXing <= bX)
0390           {
0391             if (_fAvg == 1)
0392             {
0393               z_prim[iz] = _hit_z + z_bias_avg;
0394               z_ibf[iz] = -1.055 * m + z_bias_avg;
0395             }
0396             else
0397             {
0398               z_prim[iz] = _hit_z + (bX - _event_bunchXing) * sphenix_constants::time_between_crossings * vIon * ns;
0399               z_ibf[iz] = -1.055 * m + (bX - _event_bunchXing) * sphenix_constants::time_between_crossings * vIon * ns;
0400             }
0401             if (z_prim[iz] < 0 && z_prim[iz] > -1.055 * m)
0402             {
0403               f_fill_prim[iz] = 1;
0404             }
0405             if (z_ibf[iz] < 0 && z_ibf[iz] > -1.055 * m)
0406             {
0407               f_fill_ibf[iz] = 1;
0408             }
0409           }
0410         }
0411       }
0412 
0413       double w_prim = _hit_eion * Tpc_ElectronsPerGeV;
0414       for (int iz = 0; iz < 30; iz++)
0415       {
0416         if (f_fill_prim[iz] == 1)
0417         {
0418           _h_SC_prim[iz]->Fill(_hit_phi, _hit_r, z_prim[iz], w_prim);
0419         }
0420         if (f_fill_ibf[iz] == 1)
0421         {
0422           if (!_isOnPlane)
0423           {
0424             // Redistribute charges
0425             // std::vector<double> newWeights = getNewWeights(_h_SC_ibf[iz], _h_modules_anode, _h_modules_measuredibf, _hit_r, _hit_phi, dr_bin, dphi_bin, _fUseIBFMap);
0426             std::vector<double> newWeights = getNewWeights(_h_SC_ibf[iz], _h_modules_anode, _h_modules_measuredibf, new_r, new_phi, dr_bin, dphi_bin, _fUseIBFMap);
0427 
0428             double w_ibf_tmp = newWeights[0];
0429             double w_gain_tmp = newWeights[1];
0430             _hit_r = newWeights[2];
0431             _hit_phi = newWeights[3];
0432 
0433             _ibf_vol = N_electrons * w_gain_tmp * _ampGain * w_ibf_tmp * _ampIBFfrac;
0434             _h_SC_ibf[iz]->Fill(_hit_phi, _hit_r, z_ibf[iz], _ibf_vol);
0435             if (iz == 0)
0436             {
0437               _h_SC_XY->Fill(_hit_r * cos(_hit_phi), _hit_r * sin(_hit_phi));  //,_ibf_vol);
0438             }
0439           }
0440           else
0441           {
0442             _h_SC_ibf[iz]->Fill(new_phi, new_r, z_ibf[iz], _ibf_vol);
0443             if (iz == 0)
0444             {
0445               _h_SC_XY->Fill(new_r * cos(new_phi), new_r * sin(new_phi));  //,_ibf_vol);
0446             }
0447           }
0448         }
0449       }
0450 
0451       if (f_fill_ibf[0] == 1)
0452       {
0453         _h_R->Fill(_hit_r);
0454       }
0455     }
0456   }
0457   else
0458   {
0459     if (_fSliming == 1)
0460     {
0461       _rawHits->Fill();
0462     }
0463   }
0464   _h_hits->Fill(n_hits);
0465 
0466   return Fun4AllReturnCodes::EVENT_OK;
0467 }
0468 
0469 //____________________________________________________________________________..
0470 int fillSpaceChargeMaps::End(PHCompositeNode * /*topNode*/)
0471 {
0472   std::cout << "fillSpaceChargeMaps::End" << std::endl;
0473   if (_fSliming == 1)
0474   {
0475     outfile->cd();
0476     outfile->Write();
0477     outfile->Close();
0478     delete outfile;
0479     for (int iz = 0; iz < 30; iz++)
0480     {
0481       _h_SC_prim[iz]->Sumw2(false);
0482       _h_SC_ibf[iz]->Sumw2(false);
0483     }
0484 
0485     _h_SC_XY->Sumw2(false);
0486     _h_hits->Sumw2(false);
0487     _h_DC_E->Sumw2(false);
0488     _h_R->Sumw2(false);
0489     hm->dumpHistos(_filename, "UPDATE");
0490   }
0491   else
0492   {
0493     std::cout << "Writing File:" << _filename << std::endl;
0494     hm->dumpHistos(_filename, "RECREATE");
0495   }
0496 
0497   return 0;
0498 }
0499 
0500 void fillSpaceChargeMaps::SetFrequency(int freq)
0501 {
0502   _freqKhz = freq;
0503   std::cout << "Frequency is set to: " << _freqKhz << " kHz" << std::endl;
0504 }
0505 
0506 void fillSpaceChargeMaps::SetBeamXing(const std::vector<int> &beamXs)
0507 {
0508   _beamxing = beamXs;
0509   std::cout << "Initial BeamXing is set to: " << _beamxing[0] << std::endl;
0510 }
0511 // void fillSpaceChargeMaps::SetBeamXingEnd(std::vector<int> beamXs_end){
0512 //   _beamxing_end = beamXs_end;
0513 //   std::cout<<"Last BeamXing to fill TPC is set to: "<<_beamxing_end[0]<<std::endl;
0514 //
0515 // }
0516 void fillSpaceChargeMaps::SetEvtStart(int newEvtStart)
0517 {
0518   _evtstart = newEvtStart;
0519   std::cout << "Starting event is set to: " << newEvtStart << std::endl;
0520 }
0521 
0522 void fillSpaceChargeMaps::SetUseIBFMap(bool useIBFMap)
0523 {
0524   _fUseIBFMap = useIBFMap;
0525   std::cout << "Using IBF and Gain Maps" << std::endl;
0526 }
0527 void fillSpaceChargeMaps::SetGain(float ampGain)
0528 {
0529   _ampGain = ampGain;
0530   std::cout << "Gain is set to: " << _ampGain << std::endl;
0531 }
0532 void fillSpaceChargeMaps::SetIBF(float ampIBFfrac)
0533 {
0534   _ampIBFfrac = ampIBFfrac;
0535   std::cout << "IBF is set to: " << _ampIBFfrac << std::endl;
0536 }
0537 
0538 void fillSpaceChargeMaps::SetCollSyst(int coll_syst)
0539 {
0540   _collSyst = coll_syst;
0541   static const std::string s_syst[2] = {"AA", "pp"};
0542   std::cout << "Collision system is set to: " << s_syst[_collSyst] << std::endl;
0543 }
0544 
0545 void fillSpaceChargeMaps::SetAvg(int fAvg)
0546 {
0547   _fAvg = fAvg;
0548   static const std::string s_avg[2] = {"OFF", "ON"};
0549   std::cout << "Averaging is set to: " << s_avg[_fAvg] << std::endl;
0550 }
0551 void fillSpaceChargeMaps::UseSliming(int fSliming)
0552 {
0553   _fSliming = fSliming;
0554   static const std::string s_sliming[2] = {"OFF", "ON"};
0555   std::cout << "Sliming is set to: " << s_sliming[_fSliming] << std::endl;
0556 }
0557 
0558 void fillSpaceChargeMaps::UseFieldMaps(int shiftElectrons)
0559 {
0560   _shiftElectrons = shiftElectrons;
0561   static const std::string s_shiftElectrons[2] = {"OFF", "ON"};
0562   std::cout << "Use Field-Maps is set to: " << s_shiftElectrons[_shiftElectrons] << std::endl;
0563 }
0564 
0565 std::vector<double> fillSpaceChargeMaps::getNewWeights(TH3 *h_SC_ibf, TH2 *h_modules_anode, TH2 *h_modules_measuredibf, double hit_r, double hit_phi, double dr_bin, double dphi_bin, bool fUseIBFMap)
0566 {
0567   double w_ibf_tmp = 1.0;
0568   double w_gain_tmp = 1.0;
0569   int r_bin = _h_R->GetXaxis()->FindBin(dr_bin);
0570   double r_bin_c = _h_R->GetXaxis()->GetBinCenter(r_bin);
0571   double r_bin_r = _h_R->GetXaxis()->GetBinCenter(r_bin + 1);
0572   double r_bin_l = _h_R->GetXaxis()->GetBinCenter(r_bin - 1);
0573 
0574   int phi_bin = h_SC_ibf->GetXaxis()->FindBin(dphi_bin);
0575 
0576   double phi_bin_c = h_SC_ibf->GetXaxis()->GetBinCenter(phi_bin);
0577   double phi_bin_r = h_SC_ibf->GetXaxis()->GetBinCenter(phi_bin + 1);
0578   double phi_bin_l = h_SC_ibf->GetXaxis()->GetBinCenter(phi_bin - 1);
0579 
0580   if (dr_bin > 0 && dphi_bin < 0)
0581   {
0582     if (hit_r - r_bin_c > 0)
0583     {
0584       hit_r = r_bin_r;  // hit_r + dr_bin;
0585     }
0586     else
0587     {
0588       hit_r = r_bin_l;  // hit_r - dr_bin;
0589     }
0590   }
0591 
0592   if (dr_bin < 0 && dphi_bin > 0)
0593   {
0594     if (hit_phi - phi_bin_c > 0)
0595     {
0596       hit_phi = phi_bin_r;  // hit_phi + dphi_bin;
0597     }
0598     else
0599     {
0600       hit_phi = phi_bin_l;  // hit_phi - dphi_bin;
0601     }
0602   }
0603   if (dr_bin > 0 && dphi_bin > 0)
0604   {
0605     if (hit_phi - phi_bin_c >= 0 && hit_r - r_bin_c >= 0)
0606     {
0607       hit_phi = phi_bin_r;  // hit_phi + dphi_bin;
0608       hit_r = r_bin_r;      // hit_r + dr_bin;
0609     }
0610     if (hit_phi - phi_bin_c <= 0 && hit_r - r_bin_c <= 0)
0611     {
0612       hit_phi = phi_bin_l;  // hit_phi - dphi_bin;
0613       hit_r = r_bin_l;      // hit_r - dr_bin;
0614     }
0615     if (hit_phi - phi_bin_c >= 0 && hit_r - r_bin_c <= 0)
0616     {
0617       hit_phi = phi_bin_r;  // hit_phi + dphi_bin;
0618       hit_r = r_bin_l;      // hit_r - dr_bin;
0619     }
0620     if (hit_phi - phi_bin_c <= 0 && hit_r - r_bin_c >= 0)
0621     {
0622       hit_phi = phi_bin_l;  // hit_phi - dphi_bin;
0623       hit_r = r_bin_r;      // hit_r + dr_bin;
0624     }
0625   }
0626   if (fUseIBFMap)
0627   {
0628     double x_tmp = hit_r * cos(hit_phi);
0629     double y_tmp = hit_r * sin(hit_phi);
0630     int bin_x = h_modules_anode->GetXaxis()->FindBin(x_tmp);
0631     int bin_y = h_modules_anode->GetYaxis()->FindBin(y_tmp);
0632     w_ibf_tmp = h_modules_measuredibf->GetBinContent(bin_x, bin_y);
0633     w_gain_tmp = h_modules_anode->GetBinContent(bin_x, bin_y);
0634   }
0635   // for(int p=0;p<24;p++){
0636   //   if(hit_phi>=phi_dead_bins[p] && hit_phi<=phi_dead_bins[p+1]){
0637   //     std::cout<<
0638   //     " dphi_bin  = "<< dphi_bin <<
0639   //     " phi_bin_c = " <<phi_bin_c <<
0640   //     " phi_bin_r = " <<phi_bin_r <<
0641   //     " phi_bin_l = " <<phi_bin_l
0642   //     << std::endl;
0643   //     std::cout<<"end new Weights hit_phi = "<<hit_phi<<std::endl;
0644   //   }
0645   //   p++;
0646   // }
0647   std::vector<double> newWeights;
0648   newWeights.push_back(w_ibf_tmp);
0649   newWeights.push_back(w_gain_tmp);
0650   newWeights.push_back(hit_r);
0651   newWeights.push_back(hit_phi);
0652   return newWeights;
0653 }
0654 
0655 bool fillSpaceChargeMaps::IsOverFrame(double r, double phi)
0656 {
0657   // these parameters are taken from Feb 12 drawings of frames.
0658   // double tpc_frame_side_gap = 0.8 * mm;    //mm //space between radial line and start of frame
0659   // double tpc_frame_side_width = 2.6 * mm;  //mm //thickness of frame
0660   double tpc_margin = 2.0 * mm;  // mm // extra gap between edge of frame and start of GEM holes
0661 
0662   double tpc_frame_r3_outer = 759.11 * mm;  // 758.4;//mm inner edge of larger-r frame of r3
0663   double tpc_frame_r3_inner = 583.67 * mm;  // 583.5;//mm outer edge of smaller-r frame of r3
0664 
0665   double tpc_frame_r2_outer = 574.76 * mm;  // 574.9;//mm inner edge of larger-r frame of r3
0666   double tpc_frame_r2_inner = 411.53 * mm;  // 411.4;//mm outer edge of smaller-r frame of r3
0667 
0668   double tpc_frame_r1_outer = 402.49 * mm;  // 402.6;//mm inner edge of larger-r frame of r3
0669   double tpc_frame_r1_inner = 217.83 * mm;  // 221.0;//mm outer edge of smaller-r frame of r3
0670 
0671   // double tpc_sec0_phi=0.0;//get_double_param("tpc_sec0_phi");
0672 
0673   // if the coordinate is in the radial spaces of the frames, return true:
0674   if (r <= tpc_frame_r1_inner + tpc_margin)
0675   {
0676     return true;
0677   }
0678   if (r >= tpc_frame_r1_outer - tpc_margin && r <= tpc_frame_r2_inner + tpc_margin)
0679   {
0680     return true;
0681   }
0682   if (r >= tpc_frame_r2_outer - tpc_margin && r <= tpc_frame_r3_inner + tpc_margin)
0683   {
0684     return true;
0685   }
0686   if (r >= tpc_frame_r3_outer - tpc_margin)
0687   {
0688     return true;
0689   }
0690 
0691   for (int p = 0; p < 24; p = p + 2)
0692   {
0693     if (phi >= phi_dead_bins[p] && phi <= phi_dead_bins[p + 1])
0694     {
0695       return true;
0696     }
0697   }
0698   return false;
0699 }
0700 
0701 std::vector<double> fillSpaceChargeMaps::putOnPlane(double r, double phi)
0702 {
0703   // these parameters are taken from Feb 12 drawings of frames.
0704   // double tpc_frame_side_gap = 0.8*mm;    //mm //space between radial line and start of frame
0705   // double tpc_frame_side_width = 2.6*mm;  //mm //thickness of frame
0706   double tpc_margin = 2.0 * mm;  // mm // extra gap between edge of frame and start of GEM holes
0707 
0708   double tpc_frame_r3_outer = 759.11 * mm;  // 758.4;//mm inner edge of larger-r frame of r3
0709   double tpc_frame_r3_inner = 583.67 * mm;  // 583.5;//mm outer edge of smaller-r frame of r3
0710 
0711   double tpc_frame_r2_outer = 574.76 * mm;  // 574.9;//mm inner edge of larger-r frame of r3
0712   double tpc_frame_r2_inner = 411.53 * mm;  // 411.4;//mm outer edge of smaller-r frame of r3
0713 
0714   double tpc_frame_r1_outer = 402.49 * mm;  // 402.6;//mm inner edge of larger-r frame of r3
0715   double tpc_frame_r1_inner = 217.83 * mm;  // 221.0;//mm outer edge of smaller-r frame of r3
0716 
0717   // double tpc_sec0_phi=0.0;//get_double_param("tpc_sec0_phi");
0718   double r_0_bin = -1;
0719   double phi_0_bin = -1;
0720   // if the coordinate is in the radial spaces of the frames, return true:
0721   if (r >= tpc_frame_r1_inner - 2 * tpc_margin && r <= tpc_frame_r1_inner + tpc_margin)
0722   {
0723     // Collect charges from inner margin + 1 mm
0724     r_0_bin = tpc_frame_r1_inner - tpc_margin / 2;  //- r;
0725   }
0726   if (r >= tpc_frame_r1_outer - tpc_margin && r <= tpc_frame_r2_inner + tpc_margin)
0727   {
0728     r_0_bin = tpc_frame_r2_inner - (tpc_frame_r2_inner - tpc_frame_r1_outer) / 2;
0729   }
0730   if (r >= tpc_frame_r2_outer - tpc_margin && r <= tpc_frame_r3_inner + tpc_margin)
0731   {
0732     r_0_bin = tpc_frame_r3_inner - (tpc_frame_r3_inner - tpc_frame_r2_outer) / 2;
0733   }
0734   if (r >= tpc_frame_r3_outer - tpc_margin)
0735   {
0736     // Collect charges from outer margin + 1 mm
0737     r_0_bin = tpc_frame_r3_outer + tpc_margin / 2;
0738   }
0739 
0740   // if the coordinate is within gap+width of a sector boundary, return true:
0741   // note that this is not a line of constant radius, but a linear distance from a radius.
0742 
0743   // find the two spokes we're between:
0744   for (int p = 0; p < 24; p = p + 2)
0745   {
0746     if (phi >= phi_dead_bins[p] && phi <= phi_dead_bins[p + 1])
0747     {
0748       phi_0_bin = phi_dead_bins[p] + (phi_dead_bins[p + 1] - phi_dead_bins[p]) / 2;
0749     }
0750   }
0751 
0752   // // float sectorangle = (M_PI / 6.);
0753   // // float nsectors = (phi+sectorangle/2) / sectorangle;
0754   // // if(nsectors>11) nsectors-=11;
0755   // // int nsec = floor(nsectors);
0756   // // float reduced_phi = phi - nsec * sectorangle;  //between zero and sixty degrees.
0757   // // float dist_to_previous = r * sin(reduced_phi );
0758   // // float dist_to_next = r * sin(sectorangle - reduced_phi );
0759   // // if (dist_to_previous <= tpc_frame_side_gap + tpc_frame_side_width + tpc_margin)
0760   // // {
0761   // //   //phi_0_bin = 0.0136;
0762   // //   phi_0_bin = asin((tpc_frame_side_gap + tpc_frame_side_width + tpc_margin) / r);
0763   // // }
0764   // // if (dist_to_next <= tpc_frame_side_gap + tpc_frame_side_width + tpc_margin)
0765   // // {
0766   // //   //phi_0_bin = 0.0136;
0767   // //   phi_0_bin = asin((tpc_frame_side_gap + tpc_frame_side_width + tpc_margin) / r);
0768   // // }
0769   std::vector<double> r_phi_bin;
0770   r_phi_bin.push_back(r_0_bin);
0771   r_phi_bin.push_back(phi_0_bin);
0772   return r_phi_bin;
0773 }
[ 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!