Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fillSpaceChargeMaps.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file fillSpaceChargeMaps.cc
1 #include "fillSpaceChargeMaps.h"
2 
3 #include "Shifter.h"
4 
5 #include <g4main/PHG4Hit.h>
6 #include <g4main/PHG4HitContainer.h>
7 
8 #include <fun4all/Fun4AllHistoManager.h>
9 #include <fun4all/Fun4AllReturnCodes.h>
10 #include <fun4all/SubsysReco.h> // for SubsysReco
11 
12 #include <phool/getClass.h>
13 #include <phool/PHCompositeNode.h>
14 
15 //#include <g4tpc/PHG4TpcPadPlane.h>
16 //#include <g4tpc/PHG4TpcPadPlaneReadout.h>
17 
18 #include <TAxis.h> // for TAxis
19 #include <TFile.h>
20 #include <TH1.h>
21 #include <TH2.h>
22 #include <TH3.h>
23 #include <TTree.h>
24 #include <TVector3.h>
25 
26 #include <algorithm> // for max
27 //#include <cmath> // for sin, asin, cos, floor, M_PI
28 #include <cstdio> // for sprintf, printf
29 #include <cstdlib>
30 #include <fstream>
31 #include <iostream>
32 #include <map>
33 #include <sstream>
34 #include <string>
35 #include <utility> // for pair
36 #include <vector>
37 
38 //____________________________________________________________________________..
39 fillSpaceChargeMaps::fillSpaceChargeMaps(const std::string &name, const std::string &filename)
40  : SubsysReco(name)
41  , _filename(filename)
42 {
43  std::cout << "fillSpaceChargeMaps::fillSpaceChargeMaps(const std::string &name) Calling ctor" << std::endl;
44 }
45 
46 //____________________________________________________________________________..
47 fillSpaceChargeMaps::~fillSpaceChargeMaps()
48 {
49  //std::cout << "fillSpaceChargeMaps::~fillSpaceChargeMaps() Calling dtor" << std::endl;
50  // delete whatever is created
51  delete hm;
52 }
53 
54 //____________________________________________________________________________..
55 int fillSpaceChargeMaps::Init(PHCompositeNode * /*topNode*/)
56 {
57  int nz = 72;
58  double z_rdo = 108 * cm;
59  outfile = new TFile(_filename.c_str(), "RECREATE");
60 
61  hm = new Fun4AllHistoManager("HITHIST");
62  const int r_bins_N = 66; //51;
63  double r_bins[r_bins_N + 1] = {217.83,
64  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,
65  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,
66  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,
67  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};
68  double r_bins_edges[r_bins_N + 3] = {217.83-2,217.83,
69  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,
70  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,
71  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,
72  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};
73 
74  const int nphi = 205;
75 
76  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,
77  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,
78  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,
79  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,
80  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,
81  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,
82  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,
83  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,
84  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,
85  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,
86  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,
87  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,
88  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,
89  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,
90  6.2764, 2 * M_PI};
91 
92 
93  double z_bins[2 * nz + 1];
94  for (int z = 0; z <= 2 * nz; z++)
95  {
96  z_bins[z] = -z_rdo + z_rdo / nz * z;
97  }
98 
99  _h_R = new TH1F("_h_R", "_h_R;R, [mm]", r_bins_N+2, r_bins_edges);
100  _h_hits = new TH1F("_h_hits", "_h_hits;N, [hit]", 1000, 0, 1e6);
101  _h_DC_E = new TH2F("_h_DC_E", "_h_DC_E;SC;E#times10^{6}", 2000, -100, 2e5 - 100, 1000, -50, 1e3 - 50);
102  _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);
103 
104  char name[100];
105  char name_ax[100];
106  for (int iz = 0; iz < 30; iz++)
107  {
108  sprintf(name, "_h_SC_ibf_%d", iz);
109  sprintf(name_ax, "_h_SC_ibf_%d;#phi, [rad];R, [mm];Z, [mm]", iz);
110  _h_SC_ibf[iz] = new TH3F(name, name_ax, nphi, phi_bins, r_bins_N, r_bins, 2 * nz, z_bins);
111  sprintf(name, "_h_SC_prim_%d", iz);
112  sprintf(name_ax, "_h_SC_prim_%d;#phi, [rad];R, [mm];Z, [mm]", iz);
113  _h_SC_prim[iz] = new TH3F(name, name_ax, nphi, phi_bins, r_bins_N, r_bins, 2 * nz, z_bins);
114 
115  hm->registerHisto(_h_SC_prim[iz]);
116  hm->registerHisto(_h_SC_ibf[iz]);
117  }
118  hm->registerHisto(_h_hits);
119  hm->registerHisto(_h_R);
120  hm->registerHisto(_h_DC_E);
121  hm->registerHisto(_h_SC_XY);
122  //_event_timestamp = 0;
123  _hit_eion = 0;
124  _hit_r = 0;
125  _hit_phi = 0;
126  _hit_z = 0;
127  _ibf_vol = 0;
128  _amp_ele_vol = 0;
129  if (_fSliming == 1)
130  {
131  _rawHits = new TTree("hTree", "tpc hit tree for ionization");
132  _rawHits->Branch("isOnPlane", &_isOnPlane);
133  _rawHits->Branch("hit_z", &_hit_z);
134  _rawHits->Branch("hit_r", &_hit_r);
135  _rawHits->Branch("hit_phi", &_hit_phi);
136  _rawHits->Branch("hit_eion", &_hit_eion);
137  _rawHits->Branch("ibf_vol", &_ibf_vol);
138  _rawHits->Branch("amp_ele_vol", &_amp_ele_vol);
139 
140  _rawHits->Branch("event_timestamp", &_event_timestamp);
141  _rawHits->Branch("event_bunchXing", &_event_bunchXing);
142  }
143  //padplane = new PHG4TpcPadPlaneReadout;
144  //seggeo = new PHG4TpcCylinderGeomContainer();
145  return 0;
146 }
147 
148 //____________________________________________________________________________..
149 int fillSpaceChargeMaps::InitRun(PHCompositeNode * /*topNode*/)
150 {
151  //std::cout << "fillSpaceChargeMaps::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
152  std::string line;
153  //AA collisions timestamps
154  std::string txt_file = "./data/timestamps_50kHz_1M.txt";
155  int start_line = 3;
156  if (_collSyst == 1)
157  {
158  //pp collisions timestamps
159  txt_file = "/phenix/u/hpereira/sphenix/work/g4simulations/timestamps_3MHz.txt";
160  start_line = 2;
161  }
162  std::ifstream InputFile(txt_file);
163  if (InputFile.is_open())
164  {
165  int n_line = 0;
166  while (getline(InputFile, line))
167  {
168  n_line++;
169  if (n_line > start_line)
170  {
171  std::istringstream is(line);
172  double n[2] = {0, 0};
173  int i = 0;
174  while (is >> n[i])
175  {
176  i++;
177  }
178  _timestamps[n[0]] = n[1];
179  if (n_line < 10)
180  {
181  std::cout << n[1] << std::endl;
182  }
183  _keys.push_back(int(n[0]));
184  }
185  }
186  InputFile.close();
187  }
188 
189  else
190  std::cout << "Unable to open file:" << txt_file << std::endl;
191 
192  if (_fUseIBFMap)
193  {
194  TFile *MapsFile = new TFile("./data/IBF_Map.root", "READ");
195  if (MapsFile->IsOpen()) printf("File opened successfully\n");
196  _h_modules_anode = (TH2F *) MapsFile->Get("h_modules_anode")->Clone("_h_modules_anode");
197  _h_modules_measuredibf = (TH2F *) MapsFile->Get("h_modules_measuredibf")->Clone("_h_modules_measuredibf");
198  }
199 
200  _mbRate = _freqKhz * kHz;
201  _xingRate = 9.383 * MHz;
202  _mean = mbRate / xingRate;
203  //padplane->CreateReadoutGeometry( PHCompositeNode *, seggeo);
204  return 0;
205 }
206 
207 //____________________________________________________________________________..
208 int fillSpaceChargeMaps::process_event(PHCompositeNode *topNode)
209 {
210  //double tpc_frame_r3_outer = 759.11*mm; //758.4;//mm inner edge of larger-r frame of r3
211  //double tpc_frame_r3_inner = 583.67*mm; //583.5;//mm outer edge of smaller-r frame of r3
212 
213  //double tpc_frame_r2_outer = 574.76*mm; //574.9;//mm inner edge of larger-r frame of r3
214  //double tpc_frame_r2_inner = 411.53*mm; //411.4;//mm outer edge of smaller-r frame of r3
215 
216  //double tpc_frame_r1_outer = 402.49*mm; //402.6;//mm inner edge of larger-r frame of r3
217  //double tpc_frame_r1_inner = 217.83*mm; //221.0;//mm outer edge of smaller-r frame of r3
218 
219 
220  double z_bias_avg = 0;
221  int bemxingsInFile = _keys.size();
222  int key = -1;
223  if (_fAvg == 1)
224  {
225  z_bias_avg = 1.055 * m * (float) rand() / RAND_MAX;
226  }
227  if (_evtstart >= bemxingsInFile) _evtstart = _evtstart - bemxingsInFile;
228  key = _keys.at(_evtstart);
229  _event_timestamp = (float) _timestamps[key] * ns; //units in seconds
230  _event_bunchXing = key;
231 
232  if (_evtstart % 10 == 0) std::cout << "_evtstart = " << _evtstart << std::endl;
233  _evtstart++;
234 
235  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");
236 
237  PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
238  int n_hits = 0;
239  if (hits)
240  {
241  PHG4HitContainer::ConstRange hit_range = hits->getHits();
242  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
243  {
244  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};
245  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};
246 
247  float hit_x0 = hit_iter->second->get_x(0);
248  float hit_y0 = hit_iter->second->get_y(0);
249  float hit_z0 = hit_iter->second->get_z(0);
250  float hit_x1 = hit_iter->second->get_x(1);
251  float hit_y1 = hit_iter->second->get_y(1);
252  float hit_z1 = hit_iter->second->get_z(1);
253 
254  float hit_eion = hit_iter->second->get_eion();
255  float N_electrons = hit_eion * Tpc_ElectronsPerGeV;
256  float x = (hit_x0 + f * (hit_x1 - hit_x0)) * cm;
257  float y = (hit_y0 + f * (hit_y1 - hit_y0)) * cm;
258  float z = (hit_z0 + f * (hit_z1 - hit_z0)) * cm;
259 
260  float r = sqrt(x * x + y * y);
261  float phi = atan2(x, y);
262  if (phi < 0) phi += 2 * M_PI;
263 
264  // Shift electrons according to the field maps
265  TVector3 oldPos(x / cm, y / cm, z / cm);
266  TVector3 newPos;
267  if (_shiftElectrons == 1)
268  {
269  if (oldPos.z() < 0)
270  {
271  oldPos.SetZ(abs(oldPos.z()));
272  newPos = shifter.ShiftForward(oldPos);
273  newPos.SetZ(newPos.z() * -1);
274  }
275  else
276  {
277  newPos = shifter.ShiftForward(oldPos);
278  }
279  }
280  else
281  {
282  newPos.SetZ(oldPos.Z());
283  newPos.SetX(oldPos.X());
284  newPos.SetY(oldPos.Y());
285  }
286 
287  //Reading IBF and Gain weights according to X-Y position
288  float w_ibf = 1.;
289  float w_gain = 1.;
290 
291  if (_fUseIBFMap)
292  {
293  int bin_x = _h_modules_anode->GetXaxis()->FindBin(x / mm);
294  int bin_y = _h_modules_anode->GetYaxis()->FindBin(y / mm);
295  w_ibf = _h_modules_measuredibf->GetBinContent(bin_x, bin_y);
296  w_gain = _h_modules_anode->GetBinContent(bin_x, bin_y);
297  }
298  float ionsPerEle = w_gain * _ampGain * w_ibf * _ampIBFfrac;
299 
300  //Check that it is on the frame
301  _isOnPlane = 0;
302  double dr_bin = -1;
303  double dphi_bin = -1;
304  double new_phi = newPos.Phi();
305  double new_r = newPos.Perp() * cm;
306  if (new_phi < 0) new_phi += 2 * M_PI;
307  if (!IsOverFrame(new_r, new_phi))
308  {
309  _isOnPlane = 1;
310  }
311  else
312  {
313  std::vector<double> r_phi_bin = putOnPlane(new_r / mm, new_phi);
314  dr_bin = r_phi_bin[0];
315  dphi_bin = r_phi_bin[1];
316  }
317  _hit_z = z;
318  _hit_r = r;
319  _hit_phi = phi;
320  _hit_eion = hit_eion;
321  _ibf_vol = N_electrons * ionsPerEle;
322  _amp_ele_vol = w_gain * _ampGain;
323  if (new_r < 210 || new_r > 770) continue;
324  if (_fSliming == 1) _rawHits->Fill();
325  double z_prim[30] = {-1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10,
326  -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10,
327  -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10};
328  double z_ibf[30] = {-1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10,
329  -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10,
330  -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10, -1 * 1e10};
331 
332  if (_hit_z >= 5 * mm && _hit_z < 1.055 * m)
333  {
334  n_hits++;
335  _h_DC_E->Fill(_ibf_vol, hit_eion * 1e6);
336  for (int iz = 0; iz < 30; iz++)
337  {
338  double bX = _beamxing[iz];
339  //double bX_end = _beamxing_end[iz];
340  if (_event_bunchXing <= bX)
341  {
342  if (_fAvg == 1)
343  {
344  z_prim[iz] = _hit_z - z_bias_avg;
345  z_ibf[iz] = 1.055 * m - z_bias_avg;
346  }
347  else
348  {
349  z_prim[iz] = _hit_z - (bX - _event_bunchXing) * 106 * vIon * ns;
350  z_ibf[iz] = 1.055 * m - (bX - _event_bunchXing) * 106 * vIon * ns;
351  }
352  if (z_prim[iz] > 0 && z_prim[iz] < 1.055 * m)
353  {
354  f_fill_prim[iz] = 1;
355  }
356  if (z_ibf[iz] > 0 && z_ibf[iz] < 1.055 * m)
357  {
358  f_fill_ibf[iz] = 1;
359  }
360  }
361  }
362  }
363 
364  if (_hit_z < -5 * mm && _hit_z > -1.055 * m)
365  {
366  n_hits++;
367  _h_DC_E->Fill(_ibf_vol, hit_eion * 1e6);
368  for (int iz = 0; iz < 30; iz++)
369  {
370  double bX = _beamxing[iz];
371  //double bX_end = _beamxing_end[iz];
372  if (_event_bunchXing <= bX)
373  {
374  if (_fAvg == 1)
375  {
376  z_prim[iz] = _hit_z + z_bias_avg;
377  z_ibf[iz] = -1.055 * m + z_bias_avg;
378  }
379  else
380  {
381  z_prim[iz] = _hit_z + (bX - _event_bunchXing) * 106 * vIon * ns;
382  z_ibf[iz] = -1.055 * m + (bX - _event_bunchXing) * 106 * vIon * ns;
383  }
384  if (z_prim[iz] < 0 && z_prim[iz] > -1.055 * m)
385  {
386  f_fill_prim[iz] = 1;
387  }
388  if (z_ibf[iz] < 0 && z_ibf[iz] > -1.055 * m)
389  {
390  f_fill_ibf[iz] = 1;
391  }
392  }
393  }
394  }
395 
396  double w_prim = _hit_eion * Tpc_ElectronsPerGeV;
397  for (int iz = 0; iz < 30; iz++)
398  {
399  if (f_fill_prim[iz] == 1)
400  {
401  _h_SC_prim[iz]->Fill(_hit_phi, _hit_r, z_prim[iz], w_prim);
402  }
403  if (f_fill_ibf[iz] == 1)
404  {
405  if (!_isOnPlane)
406  {
407  //Redistribute charges
408  //std::vector<double> newWeights = getNewWeights(_h_SC_ibf[iz], _h_modules_anode, _h_modules_measuredibf, _hit_r, _hit_phi, dr_bin, dphi_bin, _fUseIBFMap);
409  std::vector<double> newWeights = getNewWeights(_h_SC_ibf[iz], _h_modules_anode, _h_modules_measuredibf, new_r, new_phi, dr_bin, dphi_bin, _fUseIBFMap);
410 
411  double w_ibf_tmp = newWeights[0];
412  double w_gain_tmp = newWeights[1];
413  _hit_r = newWeights[2];
414  _hit_phi = newWeights[3];
415 
416  _ibf_vol = N_electrons * w_gain_tmp * _ampGain * w_ibf_tmp * _ampIBFfrac;
417  _h_SC_ibf[iz]->Fill(_hit_phi, _hit_r, z_ibf[iz], _ibf_vol);
418  if(iz==0)_h_SC_XY->Fill(_hit_r * cos(_hit_phi),_hit_r * sin(_hit_phi));//,_ibf_vol);
419 
420  }
421  else
422  {
423  _h_SC_ibf[iz]->Fill(new_phi, new_r, z_ibf[iz], _ibf_vol);
424  if(iz==0)_h_SC_XY->Fill(new_r * cos(new_phi),new_r * sin(new_phi));//,_ibf_vol);
425  }
426  }
427  }
428 
429  if (f_fill_ibf[0] == 1)
430  {
431  _h_R->Fill(_hit_r);
432  }
433  }
434  }
435  else
436  {
437  if (_fSliming == 1) _rawHits->Fill();
438  }
439  _h_hits->Fill(n_hits);
440 
441  return Fun4AllReturnCodes::EVENT_OK;
442 }
443 
444 //____________________________________________________________________________..
445 int fillSpaceChargeMaps::End(PHCompositeNode * /*topNode*/)
446 {
447  std::cout << "fillSpaceChargeMaps::End" << std::endl;
448  if (_fSliming == 1)
449  {
450  outfile->cd();
451  outfile->Write();
452  outfile->Close();
453  delete outfile;
454  for (int iz = 0; iz < 30; iz++)
455  {
456  _h_SC_prim[iz]->Sumw2(false);
457  _h_SC_ibf[iz]->Sumw2(false);
458  }
459 
460  _h_SC_XY->Sumw2(false);
461  _h_hits->Sumw2(false);
462  _h_DC_E->Sumw2(false);
463  _h_R->Sumw2(false);
464  hm->dumpHistos(_filename, "UPDATE");
465  }
466  else
467  {
468  std::cout << "Writing File:" << _filename << std::endl;
469  hm->dumpHistos(_filename, "RECREATE");
470  }
471 
472  return 0;
473 }
474 
475 void fillSpaceChargeMaps::SetFrequency(int freq)
476 {
477  _freqKhz = freq;
478  std::cout << "Frequency is set to: " << _freqKhz << " kHz" << std::endl;
479 }
480 
481 void fillSpaceChargeMaps::SetBeamXing(const std::vector<int> &beamXs)
482 {
483  _beamxing = beamXs;
484  std::cout << "Initial BeamXing is set to: " << _beamxing[0] << std::endl;
485 }
486 //void fillSpaceChargeMaps::SetBeamXingEnd(std::vector<int> beamXs_end){
487 // _beamxing_end = beamXs_end;
488 // std::cout<<"Last BeamXing to fill TPC is set to: "<<_beamxing_end[0]<<std::endl;
489 //
490 //}
491 void fillSpaceChargeMaps::SetEvtStart(int newEvtStart)
492 {
493  _evtstart = newEvtStart;
494  std::cout << "Starting event is set to: " << newEvtStart << std::endl;
495 }
496 
497 void fillSpaceChargeMaps::SetUseIBFMap(bool useIBFMap)
498 {
499  _fUseIBFMap = useIBFMap;
500  std::cout << "Using IBF and Gain Maps" << std::endl;
501 }
502 void fillSpaceChargeMaps::SetGain(float ampGain)
503 {
504  _ampGain = ampGain;
505  std::cout << "Gain is set to: " << _ampGain << std::endl;
506 }
507 void fillSpaceChargeMaps::SetIBF(float ampIBFfrac)
508 {
509  _ampIBFfrac = ampIBFfrac;
510  std::cout << "IBF is set to: " << _ampIBFfrac << std::endl;
511 }
512 
513 void fillSpaceChargeMaps::SetCollSyst(int coll_syst)
514 {
515  _collSyst = coll_syst;
516  static const std::string s_syst[2] = {"AA", "pp"};
517  std::cout << "Collision system is set to: " << s_syst[_collSyst] << std::endl;
518 }
519 
520 void fillSpaceChargeMaps::SetAvg(int fAvg)
521 {
522  _fAvg = fAvg;
523  static const std::string s_avg[2] = {"OFF", "ON"};
524  std::cout << "Averaging is set to: " << s_avg[_fAvg] << std::endl;
525 }
526 void fillSpaceChargeMaps::UseSliming(int fSliming)
527 {
528  _fSliming = fSliming;
529  static const std::string s_sliming[2] = {"OFF", "ON"};
530  std::cout << "Sliming is set to: " << s_sliming[_fSliming] << std::endl;
531 }
532 
533 void fillSpaceChargeMaps::UseFieldMaps(int shiftElectrons)
534 {
535  _shiftElectrons = shiftElectrons;
536  static const std::string s_shiftElectrons[2] = {"OFF", "ON"};
537  std::cout << "Use Field-Maps is set to: " << s_shiftElectrons[_shiftElectrons] << std::endl;
538 }
539 
540 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)
541 {
542  double w_ibf_tmp = 1.0;
543  double w_gain_tmp = 1.0;
544  int r_bin = _h_R->GetXaxis()->FindBin(dr_bin);
545  double r_bin_c = _h_R->GetXaxis()->GetBinCenter(r_bin);
546  double r_bin_r = _h_R->GetXaxis()->GetBinCenter(r_bin+1);
547  double r_bin_l = _h_R->GetXaxis()->GetBinCenter(r_bin-1);
548 
549  int phi_bin = h_SC_ibf->GetXaxis()->FindBin(dphi_bin);
550 
551  double phi_bin_c = h_SC_ibf->GetXaxis()->GetBinCenter(phi_bin);
552  double phi_bin_r = h_SC_ibf->GetXaxis()->GetBinCenter(phi_bin+1);
553  double phi_bin_l = h_SC_ibf->GetXaxis()->GetBinCenter(phi_bin-1);
554 
555 
556 
557  if (dr_bin > 0 && dphi_bin < 0)
558  {
559  if (hit_r - r_bin_c > 0)
560  {
561  hit_r = r_bin_r;//hit_r + dr_bin;
562  }
563  else
564  {
565  hit_r = r_bin_l;//hit_r - dr_bin;
566  }
567  }
568 
569  if (dr_bin < 0 && dphi_bin > 0)
570  {
571  if (hit_phi - phi_bin_c > 0)
572  {
573  hit_phi = phi_bin_r; //hit_phi + dphi_bin;
574  }
575  else
576  {
577  hit_phi = phi_bin_l; //hit_phi - dphi_bin;
578  }
579  }
580  if (dr_bin > 0 && dphi_bin > 0)
581  {
582  if (hit_phi - phi_bin_c >= 0 && hit_r - r_bin_c >= 0)
583  {
584  hit_phi = phi_bin_r; //hit_phi + dphi_bin;
585  hit_r = r_bin_r;//hit_r + dr_bin;
586  }
587  if (hit_phi - phi_bin_c <= 0 && hit_r - r_bin_c <= 0)
588  {
589  hit_phi = phi_bin_l; //hit_phi - dphi_bin;
590  hit_r = r_bin_l;//hit_r - dr_bin;
591  }
592  if (hit_phi - phi_bin_c >= 0 && hit_r - r_bin_c <= 0)
593  {
594  hit_phi = phi_bin_r; //hit_phi + dphi_bin;
595  hit_r = r_bin_l;//hit_r - dr_bin;
596  }
597  if (hit_phi - phi_bin_c <= 0 && hit_r - r_bin_c >= 0)
598  {
599  hit_phi = phi_bin_l; //hit_phi - dphi_bin;
600  hit_r = r_bin_r;//hit_r + dr_bin;
601  }
602  }
603  if (fUseIBFMap)
604  {
605  double x_tmp = hit_r * cos(hit_phi);
606  double y_tmp = hit_r * sin(hit_phi);
607  int bin_x = h_modules_anode->GetXaxis()->FindBin(x_tmp);
608  int bin_y = h_modules_anode->GetYaxis()->FindBin(y_tmp);
609  w_ibf_tmp = h_modules_measuredibf->GetBinContent(bin_x, bin_y);
610  w_gain_tmp = h_modules_anode->GetBinContent(bin_x, bin_y);
611  }
612  //for(int p=0;p<24;p++){
613  // if(hit_phi>=phi_dead_bins[p] && hit_phi<=phi_dead_bins[p+1]){
614  // std::cout<<
615  // " dphi_bin = "<< dphi_bin <<
616  // " phi_bin_c = " <<phi_bin_c <<
617  // " phi_bin_r = " <<phi_bin_r <<
618  // " phi_bin_l = " <<phi_bin_l
619  // << std::endl;
620  // std::cout<<"end new Weights hit_phi = "<<hit_phi<<std::endl;
621  // }
622  // p++;
623  //}
624  std::vector<double> newWeights;
625  newWeights.push_back(w_ibf_tmp);
626  newWeights.push_back(w_gain_tmp);
627  newWeights.push_back(hit_r);
628  newWeights.push_back(hit_phi);
629  return newWeights;
630 }
631 
632 bool fillSpaceChargeMaps::IsOverFrame(double r, double phi)
633 {
634  //these parameters are taken from Feb 12 drawings of frames.
635  //double tpc_frame_side_gap = 0.8 * mm; //mm //space between radial line and start of frame
636  //double tpc_frame_side_width = 2.6 * mm; //mm //thickness of frame
637  double tpc_margin = 2.0 * mm; //mm // extra gap between edge of frame and start of GEM holes
638 
639  double tpc_frame_r3_outer = 759.11 * mm; //758.4;//mm inner edge of larger-r frame of r3
640  double tpc_frame_r3_inner = 583.67 * mm; //583.5;//mm outer edge of smaller-r frame of r3
641 
642  double tpc_frame_r2_outer = 574.76 * mm; //574.9;//mm inner edge of larger-r frame of r3
643  double tpc_frame_r2_inner = 411.53 * mm; //411.4;//mm outer edge of smaller-r frame of r3
644 
645  double tpc_frame_r1_outer = 402.49 * mm; //402.6;//mm inner edge of larger-r frame of r3
646  double tpc_frame_r1_inner = 217.83 * mm; //221.0;//mm outer edge of smaller-r frame of r3
647 
648  //double tpc_sec0_phi=0.0;//get_double_param("tpc_sec0_phi");
649 
650  //if the coordinate is in the radial spaces of the frames, return true:
651  if (r <= tpc_frame_r1_inner + tpc_margin)
652  {
653  return true;
654  }
655  if (r >= tpc_frame_r1_outer - tpc_margin && r <= tpc_frame_r2_inner + tpc_margin)
656  {
657  return true;
658  }
659  if (r >= tpc_frame_r2_outer - tpc_margin && r <= tpc_frame_r3_inner + tpc_margin)
660  {
661  return true;
662  }
663  if (r >= tpc_frame_r3_outer - tpc_margin)
664  {
665  return true;
666  }
667 
668  for(int p=0;p<24;p=p+2){
669  if(phi>=phi_dead_bins[p] && phi<=phi_dead_bins[p+1]){
670  return true;
671  }
672  }
673  return false;
674 }
675 
676 std::vector<double> fillSpaceChargeMaps::putOnPlane(double r, double phi)
677 {
678  //these parameters are taken from Feb 12 drawings of frames.
679  //double tpc_frame_side_gap = 0.8*mm; //mm //space between radial line and start of frame
680  //double tpc_frame_side_width = 2.6*mm; //mm //thickness of frame
681  double tpc_margin = 2.0*mm; //mm // extra gap between edge of frame and start of GEM holes
682 
683  double tpc_frame_r3_outer = 759.11*mm; //758.4;//mm inner edge of larger-r frame of r3
684  double tpc_frame_r3_inner = 583.67*mm; //583.5;//mm outer edge of smaller-r frame of r3
685 
686  double tpc_frame_r2_outer = 574.76*mm; //574.9;//mm inner edge of larger-r frame of r3
687  double tpc_frame_r2_inner = 411.53*mm; //411.4;//mm outer edge of smaller-r frame of r3
688 
689  double tpc_frame_r1_outer = 402.49*mm; //402.6;//mm inner edge of larger-r frame of r3
690  double tpc_frame_r1_inner = 217.83*mm; //221.0;//mm outer edge of smaller-r frame of r3
691 
692  //double tpc_sec0_phi=0.0;//get_double_param("tpc_sec0_phi");
693  double r_0_bin = -1;
694  double phi_0_bin = -1;
695  //if the coordinate is in the radial spaces of the frames, return true:
696  if (r >= tpc_frame_r1_inner - 2*tpc_margin && r <= tpc_frame_r1_inner + tpc_margin)
697  {
698  //Collect charges from inner margin + 1 mm
699  r_0_bin = tpc_frame_r1_inner - tpc_margin/2;//- r;
700  }
701  if (r >= tpc_frame_r1_outer - tpc_margin && r <= tpc_frame_r2_inner + tpc_margin)
702  {
703  r_0_bin = tpc_frame_r2_inner-(tpc_frame_r2_inner - tpc_frame_r1_outer)/2;
704  }
705  if (r >= tpc_frame_r2_outer - tpc_margin && r <= tpc_frame_r3_inner + tpc_margin)
706  {
707  r_0_bin = tpc_frame_r3_inner-(tpc_frame_r3_inner - tpc_frame_r2_outer)/2;
708  }
709  if (r >= tpc_frame_r3_outer - tpc_margin)
710  {
711  //Collect charges from outer margin + 1 mm
712  r_0_bin = tpc_frame_r3_outer + tpc_margin/2;
713  }
714 
715  //if the coordinate is within gap+width of a sector boundary, return true:
716  //note that this is not a line of constant radius, but a linear distance from a radius.
717 
718  //find the two spokes we're between:
719  for(int p=0;p<24;p=p+2){
720  if(phi>=phi_dead_bins[p] && phi<=phi_dead_bins[p+1]){
721  phi_0_bin = phi_dead_bins[p] + (phi_dead_bins[p+1]-phi_dead_bins[p])/2;
722  }
723  }
724 
725  // // float sectorangle = (M_PI / 6.);
726  // // float nsectors = (phi+sectorangle/2) / sectorangle;
727  // // if(nsectors>11) nsectors-=11;
728  // // int nsec = floor(nsectors);
729  // // float reduced_phi = phi - nsec * sectorangle; //between zero and sixty degrees.
730  // // float dist_to_previous = r * sin(reduced_phi );
731  // // float dist_to_next = r * sin(sectorangle - reduced_phi );
732  // // if (dist_to_previous <= tpc_frame_side_gap + tpc_frame_side_width + tpc_margin)
733  // // {
734  // // //phi_0_bin = 0.0136;
735  // // phi_0_bin = asin((tpc_frame_side_gap + tpc_frame_side_width + tpc_margin) / r);
736  // // }
737  // // if (dist_to_next <= tpc_frame_side_gap + tpc_frame_side_width + tpc_margin)
738  // // {
739  // // //phi_0_bin = 0.0136;
740  // // phi_0_bin = asin((tpc_frame_side_gap + tpc_frame_side_width + tpc_margin) / r);
741  // // }
742  std::vector<double> r_phi_bin;
743  r_phi_bin.push_back(r_0_bin);
744  r_phi_bin.push_back(phi_0_bin);
745  return r_phi_bin;
746 }