StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StRefMultCorr.cxx
1 // C++ headers
2 #include <algorithm>
3 #include <fstream>
4 #include <iostream>
5 #include <iomanip>
6 #include <string>
7 #include <sstream>
8 #include <limits>
9 
10 // StRefMultCorr headers
11 #include "StRefMultCorr.h"
12 #include "Param.h"
13 
14 // ROOT headers
15 #include "TError.h"
16 #include "TRandom.h"
17 #include "TMath.h"
18 
19 ClassImp(StRefMultCorr)
20 
21 namespace {
22  typedef std::pair<Double_t, Int_t> keys;
23 }
24 
25 //_________________
26 // Default constructor
27 StRefMultCorr::StRefMultCorr(const TString name, const TString subname,
28  const TString libname) :
29  mName(name), mSubName(subname), mLibName(libname),
30  mRefX(1), mVerbose(kFALSE) {
31  mRefMult = 0 ;
32  mVz = -9999. ;
33  mRefMult_corr = -1.0;
34  mIsRu = false;
35  mIsZr = false;
36 
37  std::cout << mSubName.Data() <<" "<< mLibName.Data() << std::endl;
38 
39  // Clear all data members
40  clear() ;
41 
42  readHeaderFile() ;
43  readBadRunsFromHeaderFile() ;
44 }
45 
46 //_________________
47 // Default destructor
49  /* empty */
50 }
51 
52 //_________________
53 Int_t StRefMultCorr::getBeginRun(const Double_t energy, const Int_t year) {
54  keys key(std::make_pair(energy, year));
55 
56  // Make sure key exists
57  std::multimap<keys, Int_t>::iterator iterCheck = mBeginRun.find(key);
58  if ( iterCheck == mBeginRun.end() ) {
59  Error("StRefMultCorr::getBeginRun", "can't find energy = %1.1f, year = %d", energy, year);
60  return -1;
61  }
62 
63  std::pair<std::multimap<keys, Int_t>::iterator, std::multimap<keys, Int_t>::iterator> iterRange = mBeginRun.equal_range(key);
64 
65  return (*(iterRange.first)).second ;
66 }
67 
68 //________________
69 Int_t StRefMultCorr::getEndRun(const Double_t energy, const Int_t year) {
70  keys key(std::make_pair(energy, year));
71 
72  // Make sure key exists
73  std::multimap<keys, Int_t>::iterator iterCheck = mEndRun.find(key);
74  if (iterCheck == mEndRun.end()) {
75  Error("StRefMultCorr::getEndRun", "can't find energy = %1.1f, year = %d", energy, year);
76  return -1;
77  }
78 
79  std::pair<std::multimap<keys, Int_t>::iterator, std::multimap<keys, Int_t>::iterator> iterRange = mEndRun.equal_range(key);
80  std::multimap<keys, Int_t>::iterator iter = iterRange.second;
81  iter--;
82 
83  return (*iter).second;
84 }
85 
86 //_________________
87 void StRefMultCorr::clear() {
88  // Clear all arrays, and set parameter index = -1
89 
90  mYear.clear() ;
91  mStart_runId.clear() ;
92  mStop_runId.clear() ;
93  mStart_zvertex.clear() ;
94  mStop_zvertex.clear() ;
95  mNormalize_stop.clear() ;
96 
97  for(Int_t i=0; i<mNCentrality; i++) {
98  mCentrality_bins[i].clear() ;
99  }
100 
101  mParameterIndex = -1 ;
102 
103  for(Int_t i=0;i<mNPar_z_vertex;i++) {
104  mPar_z_vertex[i].clear() ;
105  }
106 
107  for(Int_t i=0;i<mNPar_weight;i++) {
108  mPar_weight[i].clear();
109  }
110 
111  for(Int_t i=0;i<mNPar_luminosity;i++) {
112  mPar_luminosity[i].clear();
113  }
114 
115  mBeginRun.clear() ;
116  mEndRun.clear() ;
117  mBadRun.clear() ;
118 
119  mnVzBinForWeight = 0 ;
120  mVzEdgeForWeight.clear();
121  mgRefMultTriggerCorrDiffVzScaleRatio.clear() ;
122 }
123 
124 //_________________
125 Bool_t StRefMultCorr::isBadRun(const Int_t RunId) {
126  // Return true if a given run id is bad run
127  std::vector<Int_t>::iterator iter = std::find(mBadRun.begin(), mBadRun.end(), RunId);
128 #if 0
129  if ( iter != mBadRun.end() ) {
130  // QA
131  std::cout << "StRefMultCorr::isBadRun Find bad run = " << (*iter) << std::endl;
132  }
133 #endif
134 
135  return ( iter != mBadRun.end() ) ;
136 }
137 
138 //_________________
139 void StRefMultCorr::initEvent(const UShort_t RefMult, const Double_t z,
140  const Double_t zdcCoincidenceRate) {
141  if (mVerbose) {
142  std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
143  std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
144  std::cout << "+++++++++++++++++++++ New Event +++++++++++++++++++\n";
145  std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
146  std::cout << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
147  }
148  // Set refmult, vz and corrected refmult if current (refmult,vz) are different from inputs
149  // User must call this function event-by-event before
150  // calling any other public functions
151  if ( mRefMult != RefMult || mVz != z || mZdcCoincidenceRate != zdcCoincidenceRate ) {
152  mRefMult = RefMult ;
153  mVz = z ;
154  mZdcCoincidenceRate = zdcCoincidenceRate ;
155  mRefMult_corr = getRefMultCorr(mRefMult, mVz, mZdcCoincidenceRate) ;
156  }
157 }
158 
159 //_________________
160 Bool_t StRefMultCorr::passnTofMatchRefmultCut(Double_t refmult, Double_t ntofmatch,
161  Double_t vz /* default = 0*/) const {
162 
163 
164  if (mVerbose) {
165  std::cout << "Pile up check...";
166  }
167  Double_t a0{}, a1{}, a2{}, a3{}, a4{};
168  Double_t b0{}, b1{}, b2{}, b3{}, b4{};
169  Double_t c0{}, c1{}, c2{}, c3{}, c4{};
170  Double_t refmultcutmax{};
171  Double_t refmultcutmin{};
172 
173  Bool_t notPileUp = kFALSE;
174 
175  // TODO:
176  // Reference multiplicity dependent pile up rejection parameter selection
177  // Be aware that the parameters are hardcoded (should be replaced with the
178  // arrays stored in Param.h)
179 
180  if (mRefX == 1) { // refMult
181  if( mParameterIndex>=30 && mParameterIndex<=35 ) { // Au+Au 27 GeV 2018
182  const Double_t min = 4.0;
183  const Double_t max = 5.0;
184 
185  if(ntofmatch<=2) return false;
186 
187  a0 = -0.704787625248525;
188  a1 = 0.99026234637141;
189  a2 = -0.000680713101607504;
190  a3 = 2.77035215460421e-06;
191  a4 = -4.04096185674966e-09;
192  b0 = 2.52126730672253;
193  b1 = 0.128066911940844;
194  b2 = -0.000538959206681944;
195  b3 = 1.21531743671716e-06;
196  b4 = -1.01886685404478e-09;
197  c0 = 4.79427731664144;
198  c1 = 0.187601372159186;
199  c2 = -0.000849856673886957;
200  c3 = 1.9359155975421e-06;
201  c4 = -1.61214724626684e-09;
202 
203  double refmultCenter = a0+a1*(ntofmatch)+a2*pow(ntofmatch,2)+a3*pow(ntofmatch,3)+a4*pow(ntofmatch,4);
204  double refmultLower = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
205  double refmultUpper = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
206 
207  refmultcutmin = refmultCenter - min*refmultLower;
208  refmultcutmax = refmultCenter + max*refmultUpper;
209  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
210  }
211  else if( mParameterIndex>=36 && mParameterIndex<=37 ) { // Isobars 200 GeV 2018
212 
213  // Zr+Zr
214  if(mParameterIndex==36) {
215  b0=13.5244327901538;
216  b1=1.4429201808933;
217  b2=-0.002873496957537;
218  b3=7.29172798142226e-06;
219  b4=-7.45759942317285e-09;
220  c0=-11.2781454979572;
221  c1=0.420728494449501;
222  c2=0.00184005031913895;
223  c3=-4.61008765754698e-06;
224  c4=4.28291905929182e-09;
225  }
226  // Ru+Ru
227  else if(mParameterIndex==37) {
228  b0=13.5426221755897;
229  b1=1.44261201539344;
230  b2=-0.00288428931227279;
231  b3=7.35384541646783e-06;
232  b4=-7.53407759526067e-09;
233  c0=-11.2591376113937;
234  c1=0.419541462167548;
235  c2=0.00185578651291454;
236  c3=-4.68933832723005e-06;
237  c4=4.4151761900593e-09;
238  }
239  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
240  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
241  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
242  }
243  else if ( mParameterIndex==38 ) { // Au+Au 19.6 GeV 2019
244 
245  if ( -145. <= vz && vz < -87. ) {
246  b0=33.7732676854599;
247  b1=1.75937881368933;
248  b2=-0.00285868075552296;
249  b3=8.50344260510873e-06;
250  b4=-1.19215380537174e-08;
251  c0=-22.7060232139884;
252  c1=0.863809402986806;
253  c2=-0.000368119767293671;
254  c3=5.97122714011036e-06;
255  c4=-1.27438638584224e-08;
256  }
257  else if ( -87. <= vz && vz < -29. ) {
258  b0=20.5875336291458;
259  b1=1.67371122493626;
260  b2=-0.00307534477962496;
261  b3=7.93755518827246e-06;
262  b4=-8.46257293600085e-09;
263  c0=-15.5923736275086;
264  c1=0.604206551537668;
265  c2=0.00131805594121643;
266  c3=-2.04753779335401e-06;
267  c4=5.73181898751325e-10;
268  }
269  else if ( -29. <= vz && vz < 29. ) {
270  b0=15.1015102672534;
271  b1=1.53929151189229;
272  b2=-0.00269203062814483;
273  b3=6.488759952638e-06;
274  b4=-6.06073586314757e-09;
275  c0=-13.1157864223955;
276  c1=0.504707692972168;
277  c2=0.00187997948645203;
278  c3=-4.7317012773039e-06;
279  c4=4.8234194091071e-09;
280  }
281  else if ( 29. <= vz && vz < 87. ) {
282  b0=20.7718852504153;
283  b1=1.67316129891511;
284  b2=-0.00315093393202473;
285  b3=8.35823870487966e-06;
286  b4=-9.14822467807924e-09;
287  c0=-15.9411138444366;
288  c1=0.61506063963685;
289  c2=0.0011824174541949;
290  c3=-1.48902496972716e-06;
291  c4=-2.29371463231934e-10;
292  }
293  else if ( 87. <= vz && vz <= 145. ) {
294  b0=33.4926150575549;
295  b1=1.79372677959986;
296  b2=-0.00319461487211403;
297  b3=9.56612691680003e-06;
298  b4=-1.31049413530369e-08;;
299  c0=-22.4679773305418;
300  c1=0.83906220918966;
301  c2=-0.000106213494253586;
302  c3=4.93946486222714e-06;
303  c4=-1.1450298089717e-08;
304  }
305 
306  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
307  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
308  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
309  }
310  else if ( mParameterIndex==39 ) { // Au+Au 14.6 GeV 2019
311  b0 = 36.4811873257854;
312  b1 = 1.96363692967013;
313  b2 = -0.00491528146300182;
314  b3 = 1.45179464078414e-05;;
315  b4 = -1.82634741809226e-08;
316  c0 = -16.176117733536;;
317  c1 = 0.780745107634961;
318  c2 = -2.03347057620351e-05;
319  c3 = 3.80646723724747e-06;
320  c4 = -9.43403282145648e-09;
321  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
322  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
323  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
324  }
325  else if ( mParameterIndex==40 ) { // Au+Au 200 GeV 2019
326 
327  b0=18.0459;
328  b1=1.32913;
329  b2=-0.000929385;
330  b3=1.53176e-06;
331  b4=-9.911e-10;
332  c0=-18.7481;
333  c1=0.785467;
334  c2=2.12757e-05;
335  c3=3.4805e-07;
336  c4=-3.80776e-10;
337 
338  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
339  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
340  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
341  }
342  else if ( mParameterIndex==41 ) { // Au+Au 7.7 GeV 2020
343 
344  if ( -145. <= vz && vz < -87. ) {
345  b0=39.578630496797;
346  b1=1.46561577132993;
347  b2=0.006515367058115;
348  b3=-4.06391982010589e-05;
349  b4=5.51203917383809e-08;
350  c0=-14.8817460248614;
351  c1=0.764539480062978;
352  c2=0.00368901349656326;
353  c3=-1.27602217700865e-05;
354  c4=8.02618485000158e-10;
355  }
356  else if ( -87. <= vz && vz < -29. ) {
357  b0=26.1841414192908;
358  b1=1.73354655107464;
359  b2=-0.00280668326418846;
360  b3=1.22370803379957e-05;
361  b4=-3.15068617200212e-08;
362  c0=-13.1831127837376;
363  c1=0.760227210117286;
364  c2=0.00195873375843822;
365  c3=-2.69378951644624e-06;
366  c4=-1.05344843941749e-08;
367  }
368  else if ( -29. <= vz && vz < 29. ) {
369  b0=23.3635904884101;
370  b1=1.58179764458174;
371  b2=-0.00100184372825271;
372  b3=7.76378744751984e-07;
373  b4=-6.46469867000365e-09;
374  c0=-11.4340781454132;
375  c1=0.72398407747444;
376  c2=0.00121092416745035;
377  c3=1.17875404059176e-07;
378  c4=-9.81658682040738e-09;
379  }
380  else if ( 29. <= vz && vz < 87. ) {
381  b0=29.4343991835005;
382  b1=1.48353715105631;
383  b2=0.00106271734149745;
384  b3=-9.07835076338586e-06;
385  b4=6.7722581625238e-09;
386  c0=-9.97159163811459;
387  c1=0.591000613390771;
388  c2=0.00449768928484714;
389  c3=-1.71667412152202e-05;
390  c4=1.6467383813372e-08;
391  }
392  else if ( 87. <= vz && vz <= 145. ) {
393  b0=37.0772875081557;
394  b1=1.53484162926915;
395  b2=0.00471873506675937;
396  b3=-2.94958548877277e-05;
397  b4=3.60887574265838e-08;
398  c0=-13.3927733032856;
399  c1=0.704319390196747;
400  c2=0.00485360248820988;
401  c3=-2.10416804123978e-05;
402  c4=1.92342533435503e-08;
403  }
404 
405  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
406  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
407  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
408  }
409  else if ( mParameterIndex==42 ) { // Au+Au 9.2 GeV 2020 TriggerID = 780020
410 
411  if ( -145. <= vz && vz < -87. ) {
412  b0=25.6055790979197;
413  b1=2.02528136596901;
414  b2=-0.0058370984051939;
415  b3=2.59602314466234e-05;
416  b4=-5.3014743584261e-08;
417  c0=-17.7059596791057;
418  c1=0.614538168662738;
419  c2=0.00534180935164814;
420  c3=-1.79582873880806e-05;
421  c4=1.01623054170579e-08;
422  }
423  else if ( -87. <= vz && vz < -29. ) {
424  b0=23.0160060308621;
425  b1=1.61885832757588;
426  b2=-0.00275873189631398;
427  b3=1.31262550392554e-05;
428  b4=-2.94368020941846e-08;
429  c0=-17.3591842617911;
430  c1=0.796170989774258;
431  c2=0.000670722514533827;
432  c3=3.26258075150876e-06;
433  c4=-1.60611460182112e-08;
434  }
435  else if ( -29. <= vz && vz < 29. ) {
436  b0=16.4277056306649;
437  b1=1.71652229539398;
438  b2=-0.00406847684302521;
439  b3=1.65203560938885e-05;
440  b4=-2.96250329214512e-08;
441  c0=-15.7887025834219;
442  c1=0.789786364309292;
443  c2=-0.000637115144252616;
444  c3=1.00019972792727e-05;
445  c4=-2.45208851616324e-08;
446  }
447  else if ( 29. <= vz && vz < 87. ) {
448  b0=21.2024767158778;
449  b1=1.70521848381614;
450  b2=-0.00352260930859763;
451  b3=1.60905730948817e-05;
452  b4=-3.37443468806432e-08;
453  c0=-17.1166088395929;
454  c1=0.814739436616432;
455  c2=0.000227197779215977;
456  c3=6.55397838050604e-06;
457  c4=-2.28812912596058e-08;
458  }
459  else if ( 87. <= vz && vz <= 145. ) {
460  b0=26.0970905882739;
461  b1=1.88889714311734;
462  b2=-0.00195374948885512;
463  b3=-6.14244087431038e-06;
464  b4=1.99930095058841e-08;
465  c0=-15.6624325989392;
466  c1=0.52385751891358;
467  c2=0.00794996911844969;
468  c3=-4.09239155250494e-05;
469  c4=6.40163739983216e-08;
470  }
471 
472  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
473  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
474  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
475  }
476  else if ( mParameterIndex==43 ) { // Au+Au 17.3 GeV 2021
477 
478  if ( -145. <= vz && vz < -87. ) {
479  b0=25.8023785946209;
480  b1=1.80974818833103;
481  b2=-0.00230107205687879;
482  b3=1.04069753338853e-05;
483  b4=-2.43265995270951e-08;
484  c0=-25.7628397848;
485  c1=1.15844463977968;
486  c2=-0.00285234327923795;
487  c3=1.68279361312683e-05;
488  c4=-2.89872992178789e-08;
489  }
490  else if ( -87. <= vz && vz < -29. ) {
491  b0=26.2142811336132;
492  b1=1.40180659301151;
493  b2=-0.000197781802002694;
494  b3=1.02666189094347e-06;
495  b4=-5.52762010064236e-09;
496  c0=-21.4352021999217;
497  c1=1.01067273031472;
498  c2=-0.00160328567162831;
499  c3=8.94486444751978e-06;
500  c4=-1.46093145276812e-08;
501  }
502  else if ( -29. <= vz && vz < 29. ) {
503  b0=20.1361585417616;
504  b1=1.54339163322734;
505  b2=-0.00277257992675217;
506  b3=1.01670412308599e-05;
507  b4=-1.4564482074994e-08;
508  c0=-18.0093218064881;
509  c1=0.858263071231256;
510  c2=-0.000411359635522234;
511  c3=4.21562873026016e-06;
512  c4=-8.07993954642765e-09;
513  }
514  else if ( 29. <= vz && vz < 87. ) {
515  b0=25.8570023358432;
516  b1=1.37245590215625;
517  b2=-5.45184310087876e-05;
518  b3=6.25643605701836e-07;
519  b4=-4.90542835006027e-09;
520  c0=-20.7158089395719;
521  c1=1.00148007639466;
522  c2=-0.00138806953636318;
523  c3=7.92595642206008e-06;
524  c4=-1.32107375325913e-08;
525  }
526  else if ( 87. <= vz && vz <= 145. ) {
527  b0=28.2036847494035;
528  b1=1.640750436652;
529  b2=-0.000569887807630565;
530  b3=3.95821109316978e-06;
531  b4=-1.60367555403757e-08;
532  c0=-26.3129222166004;
533  c1=1.21481523017369;
534  c2=-0.00341644731702994;
535  c3=1.84782571448044e-05;
536  c4=-3.03333077890128e-08;
537  }
538 
539  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
540  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
541  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
542  }
543  else if ( mParameterIndex==44 ) { // Au+Au 11.5 GeV 2020
544 
545  if ( -145. <= vz && vz < -87. ) {
546  b0=18.0402708948567;
547  b1=2.09478604674414;
548  b2=-0.00685576746251115;
549  b3=3.88333589216404e-05;
550  b4=-8.12179090437804e-08;
551  c0=-12.7515169659501;
552  c1=0.705235205311516;
553  c2=0.00321598985910965;
554  c3=-1.56896265545575e-05;
555  c4=2.97072869656044e-08;
556  }
557  else if ( -87. <= vz && vz < -29. ) {
558  b0=14.2601983060724;
559  b1=1.71255613728895;
560  b2=-0.00383919825526746;
561  b3=1.7756145374654e-05;
562  b4=-3.19509246865534e-08;
563  c0=-10.9408282877465;
564  c1=0.617024824873745;
565  c2=0.00264576299008488;
566  c3=-1.158420066816e-05;
567  c4=2.01763088491799e-08;
568  }
569  else if ( -29. <= vz && vz < 29. ) {
570  b0=11.1331231719184;
571  b1=1.69710478538775;
572  b2=-0.00464826171041643;
573  b3=2.02639545153783e-05;
574  b4=-3.4169236655577e-08;
575  c0=-8.82209022882564;
576  c1=0.524312884632579;
577  c2=0.00321682247003759;
578  c3=-1.35894996081641e-05;
579  c4=2.26005417512409e-08;
580  }
581  else if ( 29. <= vz && vz < 87. ) {
582  b0=14.615141872526;
583  b1=1.69217111894767;
584  b2=-0.00377600546419821;
585  b3=1.83551619792816e-05;
586  b4=-3.48332786210067e-08;
587  c0=-11.0113966446419;
588  c1=0.616128886729022;
589  c2=0.00278642638292705;
590  c3=-1.3124493295967e-05;
591  c4=2.44388293439677e-08;
592  }
593  else if ( 87. <= vz && vz <= 145. ) {
594  b0=17.988224598148;
595  b1=2.07853473508418;
596  b2=-0.00668791264313384;
597  b3=3.61562317906595e-05;
598  b4=-7.30405696800251e-08;
599  c0=-12.6730707166176;
600  c1=0.709713827776669;
601  c2=0.00318794623382361;
602  c3=-1.47530903374243e-05;
603  c4=2.55638251982488e-08;
604  }
605 
606  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
607  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
608  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
609  }
610  else {
611  notPileUp = kTRUE;
612  }
613 
614  } // if (mRefX == 1) { // refMult
615  else if ( mRefX == 5 ) { // fxtMult
616  if (mParameterIndex == 1) { // Run 19 Au+Au 4.59 GeV (sqrt(s_NN)=3.2 GeV)
617  b0=19.48;
618  b1=5.428;
619  b2=-0.007;
620  b3=-2.428e-4;
621  b4=1.197e-7;
622  c0=-13.59;
623  c1=1.515;
624  c2=0.02816;
625  c3=-1.195e-4;
626  c4=-9.639e-7;
627 
628  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
629  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
630  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
631  }
632  else if (mParameterIndex == 2) { // Run 20 Au+Au 5.75 GeV (sqrt(s_NN)=3.5 GeV)
633  b0=23.28;
634  b1=5.247;
635  b2=0.04037;
636  b3=-1.206e-3;
637  b4=5.792e-6;
638  c0=-14.82;
639  c1=1.583;
640  c2=0.02684;
641  c3=4.605e-5;
642  c4=-2.410e-6;
643 
644  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
645  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
646  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
647  }
648  else if (mParameterIndex == 3) { // Run 19 Au+Au 7.3 GeV (sqrt(s_NN)=3.9 GeV)
649  b0=31.5858;
650  b1=4.29274;
651  b2=0.0485779;
652  b3=-1.039e-3;
653  b4=4.408e-6;
654  c0=-22.7905;
655  c1=2.6578;
656  c2=-0.03112;
657  c3=1.031e-3;
658  c4=-7.434e-6;
659 
660  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
661  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
662  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
663  }
664  else if (mParameterIndex == 4) { // Run 20 Au+Au 7.3 GeV (sqrt(s_NN)=3.9 GeV)
665  b0=29.74;
666  b1=4.421;
667  b2=0.09139;
668  b3=-1.977e-3;
669  b4=9.435e-6;
670  c0=-20.53;
671  c1=2.557;
672  c2=-0.02094;
673  c3=8.943e-4;
674  c4=-6.879e-6;
675 
676  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
677  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
678  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
679  }
680  else if (mParameterIndex == 5) { // Run 20 Au+Au 31.2 GeV (sqrt(s_NN)=7.7 GeV)
681  b0=24.7299053323955;
682  b1=7.79546460550082;
683  b2=-0.000336278299464254;
684  b3=-0.000549204114892259;
685  b4=1.89274000668251e-06;
686  c0=-24.3335976220474;
687  c1=1.93207052914575;
688  c2=0.0042539477677528;
689  c3=0.000725893349545147;
690  c4=-6.29726910263091e-06;
691 
692  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
693  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
694  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
695  }
696  else if (mParameterIndex == 6) { // Run 20 Au+Au 5.75 GeV (sqrt(s_NN)=3.5 GeV)
697  b0=23.28;
698  b1=5.247;
699  b2=0.04037;
700  b3=-1.206e-3;
701  b4=5.792e-6;
702  c0=-14.82;
703  c1=1.583;
704  c2=0.02684;
705  c3=4.605e-5;
706  c4=-2.410e-6;
707 
708  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
709  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
710  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
711  }
712  else if (mParameterIndex == 7) { // Run 20 Au+Au 13.5 GeV (sqrt(s_NN)=5.2 GeV)
713  b0=18.6707;
714  b1=6.92307;
715  b2=-0.0293523;
716  b3=0.000412261;
717  b4=-4.74922e-06;
718  c0=-14.4436;
719  c1=-0.047413;
720  c2=0.100793;
721  c3=-0.00121203;
722  c4=5.59521e-06;
723 
724  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
725  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
726  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
727  }
728  else if (mParameterIndex == 8) { // Run 20 Au+Au 19.5 GeV (sqrt(s_NN)=6.2 GeV)
729  b0=25.0191;
730  b1=5.51924;
731  b2=0.0694824;
732  b3=-0.00121388;
733  b4=3.44057e-06;
734  c0=-16.9132;
735  c1=2.35278;
736  c2=-0.0341491;
737  c3=0.00131257;
738  c4=-9.0295e-06;
739 
740  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
741  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
742  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
743  }
744  else {
745  notPileUp = kTRUE;
746  }
747  }
748 
749  /*
750  else if ( mRefX == 6 ) { // refMult6
751  if ( mParameterIndex==0 ) { // d+Au 200 GeV 2021
752  b0 = 2.4412914662443033;
753  b1 = 5.523540420923605;
754  b2 = -0.16458436958697667;
755  b3 = 0.002805908341435613;
756  b4 = -1.6300934820294975e-05;
757  c0 = -0.86595124167792;
758  c1 = 0.44263208748354943;
759  c2 = 0.06024976895762696;
760  c3 = -0.0013523620327006189;
761  c4 = 1.0553696607739253e-05;
762  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
763  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
764  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
765  }
766  else {
767  notPileUp = kTRUE;
768  }
769  }
770  */
771 
772  if (mVerbose) {
773  std::cout << "\t notPileUp: ";
774  if (notPileUp) {
775  std::cout << "TRUE";
776  }
777  else {
778  std::cout << "FALSE";
779  }
780  std::cout << "\t[DONE]\n";
781  }
782 
783  return notPileUp;
784 }
785 
786 //_________________
787 Bool_t StRefMultCorr::isIndexOk() const {
788  // mParameterIndex not initialized (-1)
789  if ( mParameterIndex == -1 ) {
790  Error("StRefMultCorr::isIndexOk", "mParameterIndex = -1. Call init(const Int_t RunId) function to initialize centrality bins, corrections");
791  Error("StRefMultCorr::isIndexOk", "mParameterIndex = -1. or use valid run numbers defined in Centrality_def_%s.txt", mName.Data());
792  Error("StRefMultCorr::isIndexOk", "mParameterIndex = -1. exit");
793  std::cout << std::endl;
794  // Stop the process if invalid run number found
795  // exit(0);
796  return kFALSE;
797  }
798 
799  // Out of bounds
800  if ( mParameterIndex >= (Int_t)mStart_runId.size() ) {
801  Error("StRefMultCorr::isIndexOk",
802  Form("mParameterIndex = %d > max number of parameter set = %d. Make sure you put correct index for this energy",
803  mParameterIndex, mStart_runId.size()));
804  return kFALSE ;
805  }
806 
807  return kTRUE ;
808 }
809 
810 //_________________
811 Bool_t StRefMultCorr::isZvertexOk() const {
812  // Primary z-vertex check
813  return (mVz > mStart_zvertex[mParameterIndex] &&
814  mVz < mStop_zvertex[mParameterIndex]);
815 }
816 
817 //_________________
818 Bool_t StRefMultCorr::isRefMultOk() const {
819  // Invalid index
820  if ( !isIndexOk() ) return kFALSE ;
821 
822  // select 0-80%
823  return (mRefMult_corr > mCentrality_bins[0][mParameterIndex] &&
824  mRefMult_corr < mCentrality_bins[mNCentrality][mParameterIndex]);
825 }
826 
827 //_________________
828 Bool_t StRefMultCorr::isCentralityOk(const Int_t icent) const {
829  // Invalid centrality id
830  if ( icent < -1 || icent >= mNCentrality+1 ) return kFALSE ;
831 
832  // Invalid index
833  if ( !isIndexOk() ) return kFALSE ;
834 
835  // Special case
836  // 1. 80-100% for icent=-1
837  if ( icent == -1 ) return ( mRefMult_corr <= mCentrality_bins[0][mParameterIndex] );
838 
839  // 2. icent = mNCentrality
840  if ( icent == mNCentrality ) return ( mRefMult_corr <= mCentrality_bins[mNCentrality][mParameterIndex] );
841 
842  const Bool_t isOK = ( mRefMult_corr > mCentrality_bins[icent][mParameterIndex] &&
843  mRefMult_corr <= mCentrality_bins[icent+1][mParameterIndex] );
844  if(mVerbose && isOK) {
845  std::cout << "StRefMultCorr::isCentralityOk refmultcorr = " << mRefMult_corr
846  << " min. bin = " << mCentrality_bins[icent][mParameterIndex]
847  << " max. bin = " << mCentrality_bins[icent+1][mParameterIndex]
848  << std::endl;
849  }
850  return isOK ;
851 }
852 
853 //_________________
854 void StRefMultCorr::init(const Int_t RunId) {
855  // Reset mParameterIndex
856  mParameterIndex = -1 ;
857 
858  // call setParameterIndex
859  setParameterIndex(RunId) ;
860 }
861 
862 //_________________
863 Int_t StRefMultCorr::setParameterIndex(const Int_t RunId) {
864  // Determine the corresponding parameter set for the input RunId
865  for(UInt_t npar = 0; npar < mStart_runId.size(); npar++) {
866  if(RunId >= mStart_runId[npar] && RunId <= mStop_runId[npar]) {
867  mParameterIndex = npar ;
868  // std::cout << "StRefMultCorr::setParameterIndex Parameter set = " << mParameterIndex << " for RUN " << RunId << std::endl;
869  break ;
870  }
871  }
872 
873  // Multiple parameters/definitions for Run14/16 data
874  // Set mParameterIndex by hand
875  // For Run14 P16id production
876  // For Run16 P16ij production
877  if ( mName.CompareTo("grefmult", TString::kIgnoreCase) == 0 ) {
878  if ( mSubName.CompareTo("Run14_AuAu200_VpdMB5", TString::kIgnoreCase) == 0 ) {
879  if(RunId/1000000==15 && mLibName.CompareTo("P16id", TString::kIgnoreCase) == 0 ) {
880  mParameterIndex = 0;
881  if (mVzEdgeForWeight.empty()) {
882  setVzForWeight(nWeightVzBin_Run14_P16id,
883  WeightVzEdgeMin_Run14_P16id,
884  WeightVzEdgeMax_Run14_P16id);
885  }
886  if (mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
887  readScaleForWeight(nWeightgRefmultBin_Run14_P16id,
888  weight_VpdMB5ToVpdMB30_Run14_P16id);
889  }
890  } // if(RunId/1000000==15 && mLibName.CompareTo("P16id", TString::kIgnoreCase) == 0 )
891  } // if ( mSubName.CompareTo("Run14_AuAu200_VpdMB5", TString::kIgnoreCase) == 0 )
892  else if ( mSubName.CompareTo("Run16_AuAu200_VpdMB5", TString::kIgnoreCase) == 0 ) {
893  if(mLibName.CompareTo("P16ij", TString::kIgnoreCase) == 0) {
894  // prod.1
895  if (RunId / 1000 >= 17039 && RunId / 1000 <= 17130) {
896  mParameterIndex = 4;
897  if (mVzEdgeForWeight.empty()) {
898  setVzForWeight(nWeightVzBin_Run16_P16ij_prod1,
899  WeightVzEdgeMin_Run16_P16ij_prod1,
900  WeightVzEdgeMax_Run16_P16ij_prod1);
901  }
902  if (mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
903  readScaleForWeight(nWeightgRefmultBin_Run16_P16ij_prod1,
904  weight_VpdMB5ToVpdMBnoVtx_Run16_P16ij_prod1);
905  }
906  } // if (RunId / 1000 >= 17039 && RunId / 1000 <= 17130)
907  // prod.2
908  else if (RunId / 1000 >= 17169 && RunId / 1000 <= 17179) {
909  mParameterIndex = 5;
910  if (mVzEdgeForWeight.empty()) {
911  setVzForWeight(nWeightVzBin_Run16_P16ij_prod2,
912  WeightVzEdgeMin_Run16_P16ij_prod2,
913  WeightVzEdgeMax_Run16_P16ij_prod2);
914  }
915  if (mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
916  readScaleForWeight(nWeightgRefmultBin_Run16_P16ij_prod2,
917  weight_VpdMB5ToVpdMBnoVtx_Run16_P16ij_prod2);
918  }
919  } // else if (RunId / 1000 >= 17169 && RunId / 1000 <= 17179)
920  } // if(mLibName.CompareTo("P16ij", TString::kIgnoreCase) == 0)
921  } // else if ( mSubName.CompareTo("Run16_AuAu200_VpdMB5", TString::kIgnoreCase) == 0 )
922  else if ( mSubName.CompareTo("Run14_AuAu200_VpdMB30", TString::kIgnoreCase) == 0 ) {
923  if (mLibName.CompareTo("P16id", TString::kIgnoreCase) == 0) {
924  mParameterIndex = 1;
925  }
926  }
927  else if ( mSubName.CompareTo("Run14_AuAu200_VpdMBnoVtx_LowMid", TString::kIgnoreCase) == 0 ) {
928  if(mLibName.CompareTo("P16id", TString::kIgnoreCase) == 0) {
929  mParameterIndex = 2;
930  }
931  }
932  else if ( mSubName.CompareTo("Run14_AuAu200_VpdMBnoVtx_High", TString::kIgnoreCase) == 0 ) {
933  if(mLibName.CompareTo("P15ic", TString::kIgnoreCase) == 0) {
934  mParameterIndex = 3;
935  }
936  }
937  else if ( mSubName.CompareTo("Run16_AuAu200_VpdMBnoVtx", TString::kIgnoreCase) == 0 ) {
938  if(mLibName.CompareTo("P16ij", TString::kIgnoreCase) == 0) {
939  mParameterIndex = 6;
940  }
941  }
942  else{
943  mParameterIndex = -1;
944  }
945  }
946 
947  // Run numbers for isobar data are not successive
948  if ( mName.CompareTo("refmult", TString::kIgnoreCase) == 0 &&
949  mSubName.CompareTo("Isobar", TString::kIgnoreCase) == 0 ) {
950 
951  // std::cout <<"This is isobaric runs"<< std::endl;
952  mIsZr = mIsRu = kFALSE;
953  for (Int_t i = 0; i < nRunId_Zr; i++) {
954  if (RunId == IsobarRunId_Zr[i]) {
955  mIsZr = kTRUE;
956  mParameterIndex = 36;
957  // std::cout <<"--> Zr+Zr"<< std::endl;
958  break;
959  }
960  }
961  for (Int_t i = 0; i < nRunId_Ru; i++) {
962  if (RunId == IsobarRunId_Ru[i]) {
963  mIsRu = kTRUE;
964  mParameterIndex = 37;
965  // std::cout <<"--> Ru+Ru"<< std::endl;
966  break;
967  }
968  }
969  if(!mIsZr && !mIsRu) {
970  Error("StRefMultCorr::setParameterIndex", "RUN %d is not isobaric data", RunId);
971  }
972  }
973 
974  // std::cout <<"mParameterIndex = "<< mParameterIndex << std::endl;
975 
976  if(mParameterIndex == -1) {
977  Error("StRefMultCorr::setParameterIndex", "Parameter set does not exist for RUN %d", RunId);
978  }
979  //else std::cout << "Parameter set = " << npar_set << std::endl;
980 
981  if (mVerbose) {
982  std::cout << "Parameter index set to: " << mParameterIndex << std::endl;
983  }
984 
985  return mParameterIndex ;
986 }
987 
988 //_________________
990  // Call initEvent() first
991  return mRefMult_corr ;
992 }
993 
994 //________________
995 Double_t StRefMultCorr::luminosityCorrection(Double_t zdcCoincidenceRate) const {
996 
997  Double_t lumiCorr = 1.;
998 
999  if ( mVerbose ) {
1000  std::cout << "Estimation of luminosity correction factor..." << std::endl;
1001  std::cout << "\t ZDC coincidence rate: " << zdcCoincidenceRate << " BBC coincidence rate: " << 0 << std::endl;
1002  }
1003  // 200 GeV only. correction = 1 for all the other energies for BES-I
1004  // the above statement may not true for BES-II, since the luminosity is much higher than BES-I, add by Zaochen
1005  // better to check the <Refmult> vs ZDCX to see whether they are flat or not, add by Zaochen
1006  const Double_t par0_lum = mPar_luminosity[0][mParameterIndex] ;
1007  const Double_t par1_lum = mPar_luminosity[1][mParameterIndex] ;
1008 
1009  if( mParameterIndex==36 || mParameterIndex==37 || mParameterIndex==40 ||
1010  ( mRefX==5 && mParameterIndex==7 ) ||
1011  ( mRefX==5 && mParameterIndex==8 ) ) {
1012  // if(mYear[mParameterIndex] == 2018 && mIsZr) zdcmean = 96.9914;
1013  // if(mYear[mParameterIndex] == 2018 && mIsRu) zdcmean = 97.9927;
1014  Double_t b_prime = 1.;
1015  if(mParameterIndex==36) b_prime = 96.9914; // Zr
1016  if(mParameterIndex==37) b_prime = 97.9927; // Ru
1017  if(mParameterIndex==40) b_prime = 213.383; // AuAu 200GeV Run19
1018  if(mParameterIndex==7 ) b_prime = 106.245; // AuAu 5.2GeV FXT Run20
1019  if(mParameterIndex==8 ) b_prime = 114.041; // AuAu 6.2GeV FXT Run20
1020  lumiCorr = (par0_lum<std::numeric_limits<double>::epsilon() ) ? 1.0 : b_prime/(par0_lum+zdcCoincidenceRate*par1_lum);
1021  }
1022  else {
1023  lumiCorr = (par0_lum < std::numeric_limits<double>::epsilon() ) ? 1.0 : 1.0/(1.0 + par1_lum/par0_lum*zdcCoincidenceRate/1000.);
1024  }
1025 
1026  // from Run14, P16id, for VpdMB5/VPDMB30/VPDMB-noVtx, use refMult at ZdcX=30, other is at ZdcX=0;
1027  // -->changed by xlchen@lbl.gov, Run16 ~ 50kHz
1028  if(
1029  ( mSubName.CompareTo("Run14_AuAu200_VpdMB5", TString::kIgnoreCase) == 0 && mLibName.CompareTo("P16id", TString::kIgnoreCase) == 0 )
1030  || ( mSubName.CompareTo("Run14_AuAu200_VpdMB30", TString::kIgnoreCase) == 0 && mLibName.CompareTo("P16id", TString::kIgnoreCase) == 0 )
1031  || ( mSubName.CompareTo("Run14_AuAu200_VpdMBnoVtx_LowMid", TString::kIgnoreCase) == 0 && mLibName.CompareTo("P16id", TString::kIgnoreCase) == 0 )
1032  || ( mSubName.CompareTo("Run14_AuAu200_VpdMBnoVtx_High", TString::kIgnoreCase) == 0 && mLibName.CompareTo("P15ic", TString::kIgnoreCase) == 0 )
1033  || ( mSubName.CompareTo("Run16_AuAu200_VpdMB5", TString::kIgnoreCase) == 0 && mLibName.CompareTo("P16ij", TString::kIgnoreCase) == 0 )
1034  || ( mSubName.CompareTo("Run16_AuAu200_VpdMBnoVtx", TString::kIgnoreCase) == 0 && mLibName.CompareTo("P16ij", TString::kIgnoreCase) == 0 )
1035  ) {
1036  float zdcmean = 0;
1037  if(mYear[mParameterIndex] == 2014) zdcmean = 30.;
1038  if(mYear[mParameterIndex] == 2016) zdcmean = 50.;
1039  lumiCorr = (par0_lum==0.0) ? lumiCorr : lumiCorr*(par0_lum+par1_lum*zdcmean)/par0_lum;
1040  }
1041 
1042  if (mVerbose) {
1043  std::cout << "\tLuminosity correction factor: " << lumiCorr << std::endl;
1044  std::cout << "\t[DONE]\n";
1045  }
1046  return lumiCorr;
1047 }
1048 
1049 //________________
1050 Double_t StRefMultCorr::vzCorrection(Double_t z) const {
1051 
1052  if ( mVerbose ) {
1053  std::cout << "Estimation of acceptance correction factor..." << std::endl;
1054  std::cout << "\t Vertex z position: " << z << " mParameterIndex: " << mParameterIndex << std::endl;
1055 
1056  }
1057  Double_t vzCorr = 1.;
1058  if ( mParameterIndex < 38 ) {
1059  // Old correction based on the 6th-order polynomial fit of the high-end point
1060  // fit of refMult for the given Vz range
1061 
1062  // par0 to par5 define the parameters of a polynomial to parametrize z_vertex dependence of RefMult
1063  const Double_t par0 = mPar_z_vertex[0][mParameterIndex];
1064  const Double_t par1 = mPar_z_vertex[1][mParameterIndex];
1065  const Double_t par2 = mPar_z_vertex[2][mParameterIndex];
1066  const Double_t par3 = mPar_z_vertex[3][mParameterIndex];
1067  const Double_t par4 = mPar_z_vertex[4][mParameterIndex];
1068  const Double_t par5 = mPar_z_vertex[5][mParameterIndex];
1069  const Double_t par6 = mPar_z_vertex[6][mParameterIndex];
1070  // This parameter is usually 0, it takes care for an additional efficiency,
1071  // usually difference between phase A and phase B parameter 0
1072  const Double_t par7 = mPar_z_vertex[7][mParameterIndex];
1073 
1074  const Double_t RefMult_ref = par0; // Reference mean RefMult at z=0
1075  const Double_t RefMult_z = par0 + par1*z + par2*z*z + par3*z*z*z + par4*z*z*z*z + par5*z*z*z*z*z + par6*z*z*z*z*z*z; // Parametrization of mean RefMult vs. z_vertex position
1076  if(RefMult_z > 0.0) {
1077  vzCorr = (RefMult_ref + par7)/RefMult_z;
1078  }
1079  }
1080  else if ( mParameterIndex == 38 ) { // Au+Au 19.6 GeV Run 19
1081  // New Vz correction. All vz bins bins are normalize to the one at the center
1082  vzCorr = auau19_run19_vzCorr[ getVzWindowForVzDepCentDef() ];
1083  } // else if ( mParameterIndex == 38 )
1084  else if ( mParameterIndex == 39 ) {
1085  // New Vz correction. All vz bins bins are normalize to that at the center
1086  vzCorr = auau14_run19_vzCorr[ getVzWindowForVzDepCentDef() ];
1087  }
1088  else if ( mParameterIndex == 40 ) { // Au+Au 200 GeV Run 19
1089  // New Vz correction. All vz bins bins are normalize to that at the center
1090  vzCorr = auau200_run19_vzCorr[ getVzWindowForVzDepCentDef() ];
1091  }
1092  else if ( mParameterIndex == 41 ) { // Au+Au 7.7 GeV Run 21
1093  // New Vz correction. All vz bins bins are normalize to that at the center
1094  vzCorr = auau7_run21_vzCorr[ getVzWindowForVzDepCentDef() ];
1095  }
1096  else if ( mParameterIndex == 42 ) { // Au+Au 9.2 GeV Run 20 TriggerID = 780020
1097  // New Vz correction. All vz bins bins are normalize to that at the center
1098  vzCorr = auau9_trig2_run20_vzCorr[ getVzWindowForVzDepCentDef() ];
1099  }
1100  else if ( mParameterIndex == 43 ) { // Au+Au 17.3 GeV Run 21
1101  // New Vz correction. All vz bins bins are normalize to that at the center
1102  vzCorr = auau17_run21_vzCorr[ getVzWindowForVzDepCentDef() ];
1103  }
1104  else if ( mParameterIndex == 44 ) { // Au+Au 11.5 GeV Run 20
1105  // New Vz correction. All vz bins bins are normalize to that at the center
1106  vzCorr = auau11_run20_vzCorr[ getVzWindowForVzDepCentDef() ];
1107  }
1108 
1109  if (mVerbose) {
1110  std::cout << "\t Acceptance correction factor: " << vzCorr << std::endl;
1111  std::cout << "\t[DONE]\n";
1112  }
1113  return vzCorr;
1114 }
1115 
1116 //________________
1117 Double_t StRefMultCorr::sampleRefMult(Int_t refMult) const {
1118 
1119  if (mVerbose) {
1120  std::cout << "Sampling refMult value: " << refMult << std::endl;
1121  }
1122 
1123  Double_t refMult_d = -9999.;
1124  if( mParameterIndex>=30 && mParameterIndex<=44 ) {
1125  refMult_d = (Double_t)refMult - 0.5 + gRandom->Rndm();
1126  }
1127  else {
1128  refMult_d = (Double_t)refMult + gRandom->Rndm();
1129  }
1130 
1131  if (mVerbose) {
1132  std::cout << "\tSampled refMult value: " << refMult_d << std::endl;
1133  std::cout << "\t[DONE]\n";
1134  }
1135  return refMult_d;
1136 }
1137 
1138 //________________
1139 Double_t StRefMultCorr::getRefMultCorr(const UShort_t refMult,
1140  const Double_t z,
1141  const Double_t zdcCoincidenceRate,
1142  const UInt_t flag) const {
1143 
1144  if (mVerbose) {
1145  std::cout << "Start refMultCorr calculations" << std::endl
1146  << "Initial values refMult / vz / zdcX / flag: "
1147  << refMult << " / " << z << " / " << zdcCoincidenceRate << " / "
1148  << flag << std::endl;
1149  }
1150  // Apply correction if parameter index & z-vertex are ok
1151  if (!isIndexOk() || !isZvertexOk()) return refMult ;
1152 
1153  // Correction function for RefMult, takes into account z_vertex dependence
1154 
1155  Double_t lumiCorr = luminosityCorrection(zdcCoincidenceRate);
1156  Double_t vzCorr = vzCorrection(z);
1157  Double_t refMult_d = sampleRefMult( refMult );
1158  Double_t refMultCorr = -9999. ;
1159  switch (flag) {
1160  case 0: refMultCorr = refMult_d * lumiCorr; break;
1161  case 1: refMultCorr = refMult_d * vzCorr; break;
1162  case 2: refMultCorr = refMult_d * vzCorr * lumiCorr; break;
1163  default: {
1164  Error("StRefMultCorr::getRefMultCorr", "invalid flag, flag=%d, should be 0, 1 or 2", flag);
1165  refMultCorr = -9999.;
1166  }
1167  } // switch ( flag )
1168 
1169  if (mVerbose) {
1170  std::cout << "Final refMultCorr value: " << refMultCorr << std::endl;
1171  }
1172 
1173  return refMultCorr;
1174 }
1175 
1176 //_________________
1177 void StRefMultCorr::readScaleForWeight(const Char_t* input) {
1178  std::ifstream fin(input) ;
1179  if(!fin) {
1180  Error("StRefMultCorr::readScaleForWeight", "can't open %s", input);
1181  return;
1182  }
1183 
1184  // Users must set the vz bin size by setVzForWeight() (see below)
1185  if(mnVzBinForWeight==0) {
1186  Error("StRefMultCorr::readScaleForWeight",
1187  "Please call setVzForWeight() to set vz bin size");
1188  return;
1189  }
1190 
1191  // Do not allow multiple calls
1192  if(!mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
1193  Error("StRefMultCorr::readScaleForWeight",
1194  "scale factor has already set in the array");
1195  return;
1196  }
1197 
1198  std::cout << "StRefMultCorr::readScaleForWeight Read scale factor ..." << std::flush;
1199 
1200  while(fin) {
1201  Double_t scale[mnVzBinForWeight] ;
1202  for(Int_t i=0; i<mnVzBinForWeight; i++) {
1203  fin >> scale[i] ;
1204  }
1205 
1206  if(fin.eof()) break ;
1207 
1208  for(Int_t i=0; i<mnVzBinForWeight; i++) {
1209  mgRefMultTriggerCorrDiffVzScaleRatio.push_back(scale[i]);
1210  }
1211  }
1212  std::cout << " [OK]" << std::endl;
1213 }
1214 
1215 // NEW version to read Vz dependent weights from header
1216 // Implemented inside StRefMultCorr::setParameterIndex(RunId)
1217 //_________________
1218 void StRefMultCorr::readScaleForWeight(const Int_t nRefmultbin, const Double_t *weight) {
1219 
1220  // Users must set the vz bin size by setVzForWeight() (see below)
1221  if( mnVzBinForWeight==0 ) {
1222  Error("StRefMultCorr::readScaleForWeight",
1223  "Please call setVzForWeight() to set vz bin size");
1224  return;
1225  }
1226 
1227  // Do not allow multiple calls
1228  if(!mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
1229  Error("StRefMultCorr::readScaleForWeight",
1230  "scale factor has already set in the array");
1231  return;
1232  }
1233 
1234  std::cout << "StRefMultCorr::readScaleForWeight Read scale factor ..." << std::flush;
1235 
1236  for(Int_t i=0; i<nRefmultbin*mnVzBinForWeight; i++) {
1237  mgRefMultTriggerCorrDiffVzScaleRatio.push_back(weight[i]);
1238  }
1239 
1240  std::cout << " [OK]" << std::endl;
1241 }
1242 
1243 
1244 // In NEW version, setVzForWeight() is implemented inside StRefMultCorr::setParameterIndex(RunId)
1245 // It does not need to be called by users.
1246 //_________________
1247 void StRefMultCorr::setVzForWeight(const Int_t nbin, const Double_t min, const Double_t max) {
1248  // Do not allow multiple calls
1249  if ( !mVzEdgeForWeight.empty() ) {
1250  Error("StRefMultCorr::setVzForWeight",
1251  "z-vertex range for weight has already been defined");
1252  return;
1253  }
1254 
1255  mnVzBinForWeight = nbin ;
1256  // calculate increment size
1257  const Double_t step = (max-min)/(Double_t)nbin;
1258 
1259  for(Int_t i=0; i<mnVzBinForWeight+1; i++) {
1260  mVzEdgeForWeight.push_back( min + step*i );
1261  }
1262 
1263  // Debug
1264  for(Int_t i=0; i<mnVzBinForWeight; i++) {
1265  std::cout << i << " " << step << " " << mVzEdgeForWeight[i]
1266  << ", " << mVzEdgeForWeight[i+1] << std::endl;
1267  }
1268 }
1269 
1270 //_________________
1271 Double_t StRefMultCorr::getScaleForWeight() const {
1272  // Special scale factor for global refmult in Run14 (Run16)
1273  // to account for the relative difference of VPDMB5 w.r.t VPDMB30 (VPDMBnoVtx)
1274 
1275  // return 1 if mgRefMultTriggerCorrDiffVzScaleRatio array is empty
1276  if(mgRefMultTriggerCorrDiffVzScaleRatio.empty()) return 1.0 ;
1277 
1278  // const Int_t nVzBins =6;
1279  // Double_t VzEdge[nVzBins+1]={-6., -4., -2., 0., 2., 4., 6.};
1280 
1281  Double_t VPD5weight=1.0;
1282  for(Int_t j=0;j<mnVzBinForWeight;j++) {
1283  if(mVz>mVzEdgeForWeight[j] && mVz<=mVzEdgeForWeight[j+1]) {
1284  /*
1285  //refMultbin=mgRefMultTriggerCorrDiffVzScaleRatio_2[j]->FindBin(mRefMult_corr+1e-6);
1286  //VPD5weight=mgRefMultTriggerCorrDiffVzScaleRatio_2[j]->GetBinContent(refMultbin);
1287  const Int_t refMultbin=static_cast<Int_t>(mRefMult_corr);
1288  //VPD5weight=mgRefMultTriggerCorrDiffVzScaleRatio[j][refMultbin];
1289  VPD5weight=mgRefMultTriggerCorrDiffVzScaleRatio[refMultbin*mnVzBinForWeight + j];
1290  const Double_t tmpContent=VPD5weight;
1291  if(tmpContent==0 || (mRefMult_corr>500 && tmpContent<=0.65)) VPD5weight=1.15;//Just because the value of the weight is around 1.15
1292  if(mRefMult_corr>500 && tmpContent>=1.35) VPD5weight=1.15;//Remove those Too large weight factor,gRefmult>500
1293  // this weight and reweight should be careful, after reweight(most peripheral),Then weight(whole range)
1294  */
1295 
1296  const Int_t refMultbin=static_cast<Int_t>(mRefMult_corr);
1297  VPD5weight=mgRefMultTriggerCorrDiffVzScaleRatio[refMultbin*mnVzBinForWeight + j];
1298  const Double_t tmpContent=VPD5weight;
1299  // 1) Ratios fluctuate too much at very high gRefmult due to low statistics
1300  // 2) Avoid some events with too high weight
1301  if(mRefMult_corr>550 && (tmpContent>3.0||tmpContent<0.3)) VPD5weight=1.0;
1302  // this weight and reweight should be careful, after reweight(most peripheral),Then weight(whole range)
1303  }
1304  }
1305 
1306  return 1.0/VPD5weight;
1307 }
1308 
1309 //_________________
1310 // For Run18 27 GeV
1312 
1313  if (mVerbose) {
1314  std::cout << "Estimating shape correction factor" << std::endl;
1315  }
1316 
1317  Double_t weight = 1.;
1318  Int_t iVzBinIndex = getVzWindowForVzDepCentDef();
1319 
1320  if ( mParameterIndex>=30 && mParameterIndex<=35 ) { // Run18 27 GeV MB
1321 
1322  if(mRefMult_corr>=500.) return 1.0; // almost no Refmult>500 for this collision energy
1323 
1324  Int_t iRunPdIndex = mParameterIndex-30;
1325  Int_t iRefmultBin = (Int_t)(mRefMult_corr/2.); //find the refmult bin matching to the Parameter bin, if binWidth=2
1326  //Int_t iRefmultBin = (Int_t)(mRefMult_corr); //find the refmult bin matching to the Parameter bin, if binWidth=1
1327 
1328  if(iRunPdIndex<0 || iRunPdIndex>5) return 1.0;
1329  if(iVzBinIndex<0 || iVzBinIndex>13) return 1.0;
1330 
1331  //----------------------------------------------------------
1332  //load the ShapeReweight factors in given RunPd and VzBin
1333  //----------------------------------------------------------
1334  std::vector<std::string> sParam_ShapeWeight = StringSplit(getParamX_ShapeWeight(iRunPdIndex,iVzBinIndex),',');
1335  if( iRefmultBin>=(Int_t)sParam_ShapeWeight.size() ) {
1336  std::cout<<"ERROR: sParam_ShapeWeight is out of ranges!!!!!"<<std::endl;
1337  return 1.0;
1338  }
1339 
1340  //std::cout<<"sParam_ShapeWeight.size(): "<<sParam_ShapeWeight.size()<<std::endl;
1341  //for(UInt_t is=0; is<sParam_ShapeWeight.size(); is++) std::cout<<"sParam_ShapeWeight[is]: "<<sParam_ShapeWeight[is]<<std::endl;
1342 
1343  Double_t ShapeReweight = 1.0;
1344 
1345  Double_t tem_ShapeReweight = std::stod( sParam_ShapeWeight[iRefmultBin] );
1346  //prevent the crazy numbers for the large fluctuations
1347  if( tem_ShapeReweight<0.1 ) {
1348  ShapeReweight = 0.1;
1349  }
1350  else if(tem_ShapeReweight>10) {
1351  ShapeReweight = 10.;
1352  }
1353  else if(tem_ShapeReweight>=0.1 && tem_ShapeReweight<=10.) {
1354  ShapeReweight = tem_ShapeReweight;
1355  }
1356 
1357  weight = 1./ShapeReweight;
1358 
1359  } // if ( mParameterIndex>=30 && mParameterIndex<=35 ) { // Run18 27 GeV MB
1360  else if ( mParameterIndex>=36 && mParameterIndex<=37 ) { // Isobar collision 200 GeV 2018
1361 
1362  if (mVz >= -9 && mVz <= 9) {
1363  return 1.;
1364  }
1365 
1366  Int_t mShapeIndex = 0;
1367  if (mIsZr) {
1368  mShapeIndex = 1;
1369  }
1370 
1371  //retrive shape weight
1372  if (iVzBinIndex >= 22) {
1373  weight = ShapeWeightArray[mShapeIndex][iVzBinIndex - 9][TMath::Nint(mRefMult_corr)];
1374  }
1375  else {
1376  weight = ShapeWeightArray[mShapeIndex][iVzBinIndex][TMath::Nint(mRefMult_corr)];
1377  }
1378  //handle bad weight
1379  if (weight == 0 || TMath::IsNaN(weight)) {
1380  weight = 1.;
1381  }
1382  } // else if ( mParameterIndex>=36 && mParameterIndex<=37 ) { // Isobar collision 200 GeV 2018
1383  else if (mParameterIndex == 38) { // Au+Au 19.6 GeV 2019
1384 
1385  if (iVzBinIndex < 0 || iVzBinIndex > auau19_run19_nVzBins) return 1.0;
1386 
1387  weight = auau19_run19_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1388  // Handle bad weight
1389  if (weight == 0 || TMath::IsNaN(weight)) {
1390  weight = 1.;
1391  }
1392  } // else if ( mParameterIndex == 38 ) { // Au+Au 19.6 GeV 2019
1393  else if (mParameterIndex == 39) { // Au+Au 14.6 GeV 2019
1394 
1395  if (iVzBinIndex < 0 || iVzBinIndex > auau14_run19_nVzBins) return 1.0;
1396 
1397  weight = auau14_run19_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1398  // Handle bad weight
1399  if (weight == 0 || TMath::IsNaN(weight)) {
1400  weight = 1.;
1401  }
1402  } // else if ( mParameterIndex == 39 ) { // Au+Au 14.6 GeV 2019
1403  else if (mParameterIndex == 40) { // Au+Au 200 GeV 2019
1404 
1405  if (iVzBinIndex < 0 || iVzBinIndex > auau200_run19_nVzBins) return 1.0;
1406 
1407  weight = auau200_run19_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1408  // Handle bad weight
1409  if (weight == 0 || TMath::IsNaN(weight)) {
1410  weight = 1.;
1411  }
1412  }
1413  else if (mParameterIndex == 41) { // Au+Au 7.7 GeV 2020
1414 
1415  if (iVzBinIndex < 0 || iVzBinIndex > auau7_run21_nVzBins) return 1.0;
1416 
1417  weight = auau7_run21_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1418  // Handle bad weight
1419  if (weight == 0 || TMath::IsNaN(weight)) {
1420  weight = 1.;
1421  }
1422  }
1423  else if (mParameterIndex == 42) { // Au+Au 9.2 GeV 2020 TrigerID = 780020
1424 
1425  if (iVzBinIndex < 0 || iVzBinIndex > auau9_run20_nVzBins) return 1.0;
1426 
1427  weight = auau9_trig2_run20_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1428  // Handle bad weight
1429  if (weight == 0 || TMath::IsNaN(weight)) {
1430  weight = 1.;
1431  }
1432  }
1433  else if (mParameterIndex == 43) { // Au+Au 17.3 GeV 2021
1434 
1435  if (iVzBinIndex < 0 || iVzBinIndex > auau17_run21_nVzBins) return 1.0;
1436 
1437  weight = auau17_run21_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1438  // Handle bad weight
1439  if (weight == 0 || TMath::IsNaN(weight)) {
1440  weight = 1.;
1441  }
1442  }
1443  else if (mParameterIndex == 44) { // Au+Au 11.5 GeV 2020
1444 
1445  if (iVzBinIndex < 0 || iVzBinIndex > auau11_run20_nVzBins) return 1.0;
1446 
1447  weight = auau11_run20_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1448  // Handle bad weight
1449  if (weight == 0 || TMath::IsNaN(weight)) {
1450  weight = 1.;
1451  }
1452  }
1453 
1454  /*
1455  else if (mRefX == 6 && mParameterIndex == 0) { // d+Au 200 GeV 2021
1456 
1457  if (iVzBinIndex < 0 || iVzBinIndex > dau200_run21_nVzBins) return 1.0;
1458 
1459  weight = dau200_run21_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1460  // Handle bad weight
1461  if (weight == 0 || TMath::IsNaN(weight)) {
1462  weight = 1.;
1463  }
1464  }
1465  */
1466 
1467  else {
1468  weight = 1.0;
1469  }
1470 
1471  if (mVerbose) {
1472  std::cout << "\tShape weight: " << weight << std::endl;
1473  std::cout << "\t[DONE]\n";
1474  }
1475 
1476  return weight;
1477 }
1478 
1479 //________________
1481 
1482  if (mVerbose) {
1483  std::cout << "Estimating the trigger weight" << std::endl;
1484  }
1485 
1486  Double_t weight = 1.;
1487 
1488  const Double_t par0 = mPar_weight[0][mParameterIndex];
1489  const Double_t par1 = mPar_weight[1][mParameterIndex];
1490  const Double_t par2 = mPar_weight[2][mParameterIndex];
1491  const Double_t par3 = mPar_weight[3][mParameterIndex];
1492  const Double_t par4 = mPar_weight[4][mParameterIndex];
1493  const Double_t A = mPar_weight[5][mParameterIndex]; // Set to 0 if no z-dependent trigger weight is set
1494  const Double_t par6 = mPar_weight[6][mParameterIndex];
1495  const Double_t par7 = mPar_weight[7][mParameterIndex];
1496 
1497  if (mVerbose) {
1498  std::cout << "Trigger efficiency parameters: "
1499  << " [0]: " << par0
1500  << " [1]: " << par1
1501  << " [2]: " << par2
1502  << " [3]: " << par3
1503  << " [4]: " << par4
1504  << " A: " << A
1505  << " [6]: " << par6
1506  << " [7]: " << par7 << std::endl;
1507  }
1508 
1509  // Additional z-vetex dependent correction
1510  //const Double_t A = ((1.27/1.21))/(30.0*30.0); // Don't ask...
1511  //const Double_t A = (0.05/0.21)/(30.0*30.0); // Don't ask...
1512 
1513 
1514  if ( isRefMultOk() // 0-80%
1515  && mRefMult_corr < mNormalize_stop[mParameterIndex] // reweighting only apply up to normalization point
1516  && mRefMult_corr != -(par3/par2) ) { // avoid denominator = 0
1517 
1518  // Parametrization of MC/data RefMult ratio
1519  if (mRefX == 5 && mParameterIndex == 0) { // Run 18 Au+Au 3.85 GeV (sqrt(s_NN)=3 GeV)
1520  // Trigger efficiency correction does not exist. Temporarily set weight to 1
1521  weight = 1.;
1522  } // else if (mRefX == 5 && mParameterIndex == 0)
1523  else {
1524  weight = ( par0 +
1525  par1 / (par2 * mRefMult_corr + par3) +
1526  par4 * (par2 * mRefMult_corr + par3) +
1527  par6 / ( (par2 * mRefMult_corr + par3) * (par2 * mRefMult_corr + par3) ) +
1528  par7 * ( (par2 * mRefMult_corr + par3) * (par2 * mRefMult_corr + par3) ) );
1529  }
1530  /*
1531  std::cout << "par0: " << par0 << " par1: " << par1 << " par2: " << par2
1532  << " par3: " << par3 << " par4: " << par4 << " A: " << A
1533  << " par6: " << par6 << " par7: " << par7 << "\n"
1534  << "refMultCorr: " << mRefMult_corr << " weight: " << weight << std::endl;
1535  */
1536 
1537  weight = weight + (weight - 1.0) * (A * mVz * mVz); // z-dependent weight correction
1538  }
1539  if (mVerbose) {
1540  std::cout << "\tTrigger weight: " << weight << std::endl;
1541  std::cout << "\t[DONE]\n";
1542  }
1543  return weight;
1544 }
1545 
1546 //_________________
1547 Double_t StRefMultCorr::getWeight() const {
1548 
1549  Double_t Weight = 1.0;
1550 
1551  // Invalid index
1552  if( !isIndexOk() ) return Weight ;
1553  // Invalid z-vertex
1554  if( !isZvertexOk() ) return Weight ;
1555 
1556  // Trigger efficiency
1557  Weight = triggerWeight();
1558 
1559  //------------for Run14 and Run16----------------
1560  // Special scale factor for global refmult depending on Vz window
1561  // for others, scale factor = 1
1562  Weight *= getScaleForWeight();
1563 
1564  // Shape correction
1565  Weight *= getShapeWeight_SubVz2Center();
1566 
1567  return Weight;
1568 }
1569 
1570 //_________________
1572 
1573  if (mVerbose) {
1574  std::cout << "Finding centrality bin (out of 16)\n";
1575  }
1576  Int_t CentBin16 = -1;
1577 
1578  // Invalid index
1579  if( !isIndexOk() ) return CentBin16 ;
1580 
1581  while( CentBin16 < mNCentrality && !isCentralityOk(CentBin16) ) {
1582  CentBin16++;
1583  }
1584 
1585  // Vz dependent centrality definition for Vz<-27 and Vz>25
1586  // Run17 54.4 GeV, trigid = 580001
1587  if( mParameterIndex==28&&(mVz<-27||mVz>25) ) {
1588  CentBin16 = getCentralityBin16VzDep();
1589  // if(CentBin16==-1) std::cout <<"Vz="<< mVz <<" RefMult="<< mRefMult <<" RefMultCorr="<< mRefMult_corr << std::endl;
1590  // if(CentBin16==9999) std::cout << "Invalide number"<< std::endl;
1591  // std::cout <<"Vz dependent centrality definition for Vz<-27 and Vz>25 ... Vz="<< mVz <<" RefMult="<< mRefMult <<" RefMultCorr="<< mRefMult_corr <<" iCent="<< CentBin16 << std::endl;
1592  }
1593 
1594  if (mVerbose) {
1595  std::cout << "\tCentrality bin (out of 16): " << CentBin16 << "\t[DONE]" << std::endl;
1596  }
1597 
1598  // return -1 if CentBin16 = 16 (very large refmult, refmult>5000)
1599  return ( CentBin16==16 ) ? -1 : CentBin16;
1600 }
1601 
1602 //_________________
1604 
1605  if (mVerbose) {
1606  std::cout << "Finding centrality bin (out of 9)\n";
1607  }
1608 
1609  Int_t CentBin9 = -1;
1610 
1611  // Invalid index
1612  if ( !isIndexOk() ) return CentBin9 ;
1613 
1614  const Int_t CentBin16 = getCentralityBin16(); // Centrality bin 16
1615  const Bool_t isCentralityOk = CentBin16 >= 0 && CentBin16 < mNCentrality ;
1616 
1617  // No centrality is defined
1618  if ( !isCentralityOk ) return CentBin9;
1619 
1620  // First handle the exceptions
1621  if( mRefMult_corr > mCentrality_bins[15][mParameterIndex] &&
1622  mRefMult_corr <= mCentrality_bins[16][mParameterIndex] ) {
1623  CentBin9 = 8; // most central 5%
1624  }
1625  else if( mRefMult_corr > mCentrality_bins[14][mParameterIndex] &&
1626  mRefMult_corr <= mCentrality_bins[15][mParameterIndex]) {
1627  CentBin9 = 7; // most central 5-10%
1628  }
1629  else {
1630  CentBin9 = (Int_t)(0.5*CentBin16);
1631  }
1632 
1633  // Vz dependent centrality definition for Vz<-27 and Vz>25
1634  // Run17 54.4 GeV, trigid = 580001
1635  // std::cout << mParameterIndex << std::endl;
1636  if ( mParameterIndex == 28 && ( mVz<-27 || mVz>25 ) ) {
1637  CentBin9 = getCentralityBin9VzDep();
1638  }
1639 
1640  if (mVerbose) {
1641  std::cout << "\tCentrality bin (out of 9): " << CentBin9 << "\t[DONE]" << std::endl;
1642  }
1643 
1644  return CentBin9;
1645 }
1646 
1647 //_________________
1648 Int_t StRefMultCorr::getVzWindowForVzDepCentDef() const {
1649  Int_t iBinVz = -1;
1650 
1651  if( mParameterIndex==28 ) { // 54.4 GeV, RefMult, 580001
1652  if( mVz>-30 && mVz<-29 ) iBinVz = 0;
1653  else if( mVz>-29 && mVz<-27 ) iBinVz = 1;
1654  else if( mVz>25 && mVz<27 ) iBinVz = 2;
1655  else if( mVz>27 && mVz<29 ) iBinVz = 3;
1656  else if( mVz>29 && mVz<30 ) iBinVz = 4;
1657  else iBinVz = -1;
1658  }
1659  else if( mParameterIndex>=30 && mParameterIndex<=35 ) { //Run18 27 GeV MB, 6 triggerIds
1660  if( mVz>=-70. && mVz<-60. ) iBinVz = 0;
1661  else if( mVz>=-60.&&mVz<-50.) iBinVz = 1;
1662  else if( mVz>=-50.&&mVz<-40.) iBinVz = 2;
1663  else if( mVz>=-40.&&mVz<-30.) iBinVz = 3;
1664  else if( mVz>=-30.&&mVz<-20.) iBinVz = 4;
1665  else if( mVz>=-20.&&mVz<-10.) iBinVz = 5;
1666  else if( mVz>=-10.&&mVz<0.0 ) iBinVz = 6;
1667  else if( mVz>=0.0 &&mVz<10. ) iBinVz = 7;
1668  else if( mVz>=10. &&mVz<20. ) iBinVz = 8;
1669  else if( mVz>=20. &&mVz<30. ) iBinVz = 9;
1670  else if( mVz>=30. &&mVz<40. ) iBinVz = 10;
1671  else if( mVz>=40. &&mVz<50. ) iBinVz = 11;
1672  else if( mVz>=50. &&mVz<60. ) iBinVz = 12;
1673  else if( mVz>=60. &&mVz<=70 ) iBinVz = 13;
1674  else iBinVz = -1;
1675  }
1676  else if( mParameterIndex>=36 && mParameterIndex<=37 ) { //Run18 200 GeV isobar
1677  Double_t VtxZBinDouble = mVz/2. + 17.;
1678  iBinVz = 0;
1679  if (mVz == 25.) {
1680  iBinVz = 29;
1681  }
1682  else if (mVz == -35.) {
1683  iBinVz = 0;
1684  }
1685  else {
1686  iBinVz = TMath::Nint(VtxZBinDouble);
1687  }
1688  }
1689  else if ( (mParameterIndex == 38) || (mParameterIndex == 39) ) { // Au+Au 19.6 GeV 2019 || Au+Au 14.6 GeV 2019
1690 
1691  for ( Int_t iVz=0; iVz<auau19_run19_nVzBins; iVz++ ) { // Utilize same Vz binning: (29 bins, -145., 145.)
1692  if ( auau19_run19_vzRangeLimits[iVz][0] <= mVz && mVz < auau19_run19_vzRangeLimits[iVz][1] ) {
1693  iBinVz = iVz;
1694  break;
1695  }
1696  } // for ( Int_t iVz=0; iVz<auau19_run19_nVzBins; iVz++ )
1697  } // else if ( mParameterIndex == 38 )
1698  else if ( mParameterIndex == 40 ) { // Au+Au 200 GeV 2019
1699  for ( Int_t iVz=0; iVz<auau200_run19_nVzBins; iVz++ ) { // Utilize Vz binning: (20 bins, -100., 100.)
1700  if ( auau200_run19_vzRangeLimits[iVz][0] <= mVz && mVz < auau200_run19_vzRangeLimits[iVz][1] ) {
1701  iBinVz = iVz;
1702  break;
1703  }
1704  } // for ( Int_t iVz=0; iVz<auau200_run19_nVzBins; iVz++ )
1705  } // else if ( mParameterIndex == 40 )
1706  else if ( mParameterIndex == 41 ) { // Au+Au 7.7 GeV 2020
1707  for ( Int_t iVz=0; iVz<auau7_run21_nVzBins; iVz++ ) { // Utilize Vz binning: (29 bins, -145., 145.)
1708  if ( auau7_run21_vzRangeLimits[iVz][0] <= mVz && mVz < auau7_run21_vzRangeLimits[iVz][1] ) {
1709  iBinVz = iVz;
1710  break;
1711  }
1712  } // for ( Int_t iVz=0; iVz<auau7_run21_nVzBins; iVz++ )
1713  } // else if ( mParameterIndex == 41 )
1714  else if ( mParameterIndex == 42 ) { // Au+Au 9.2 GeV Run 20 TriggerID = 780020
1715  for ( Int_t iVz=0; iVz<auau9_run20_nVzBins; iVz++ ) { // Utilize Vz binning: (29 bins, -145., 145.)
1716  if ( auau9_run20_vzRangeLimits[iVz][0] <= mVz && mVz < auau9_run20_vzRangeLimits[iVz][1] ) {
1717  iBinVz = iVz;
1718  break;
1719  }
1720  } // for ( Int_t iVz=0; iVz<auau9_run20_nVzBins; iVz++ )
1721  } // else if ( mParameterIndex == 42 )
1722  else if ( mParameterIndex == 43 ) { // Au+Au 17.3 GeV 2021
1723  for ( Int_t iVz=0; iVz<auau17_run21_nVzBins; iVz++ ) { // Utilize Vz binning: (29 bins, -145., 145.)
1724  if ( auau17_run21_vzRangeLimits[iVz][0] <= mVz && mVz < auau17_run21_vzRangeLimits[iVz][1] ) {
1725  iBinVz = iVz;
1726  break;
1727  }
1728  } // for ( Int_t iVz=0; iVz<auau17_run21_nVzBins; iVz++ )
1729  } // else if ( mParameterIndex == 43 )
1730  else if ( mParameterIndex == 44 ) { // Au+Au 11.5 GeV 2020
1731  for ( Int_t iVz=0; iVz<auau11_run20_nVzBins; iVz++ ) { // Utilize Vz binning: (29 bins, -145., 145.)
1732  if ( auau11_run20_vzRangeLimits[iVz][0] <= mVz && mVz < auau11_run20_vzRangeLimits[iVz][1] ) {
1733  iBinVz = iVz;
1734  break;
1735  }
1736  } // for ( Int_t iVz=0; iVz<auau17_run21_nVzBins; iVz++ )
1737  } // else if ( mParameterIndex == 44 )
1738  /*
1739  else if ( mRefX == 6 && mParameterIndex == 0 ) { // d+Au 200 GeV 2021
1740  for ( Int_t iVz=0; iVz<dau200_run21_nVzBins; iVz++ ) {
1741  if ( dau200_run21_vzRangeLimits[iVz][0] <= mVz && mVz < dau200_run21_vzRangeLimits[iVz][1] ) {
1742  iBinVz = iVz;
1743  break;
1744  }
1745  } // for ( Int_t iVz=0; iVz<dau200_run21_nVzBins; iVz++ )
1746  } // else if ( mRefX == 6 && mParameterIndex == 0 )
1747  */
1748  else {
1749  iBinVz = -1;
1750  }
1751 
1752  return iBinVz;
1753 }
1754 
1755 //_________________
1756 Int_t StRefMultCorr::getCentralityBin9VzDep() const {
1757  const Int_t vzid = getVzWindowForVzDepCentDef();
1758  Int_t iCent = 9999;
1759  for(Int_t i=0; i<9; i++) {
1760  if ( i==8 ) {
1761  if( mRefMult_corr>CentBin9_vzdep[vzid][i] && mRefMult_corr<50000 ) iCent = i;
1762  }
1763  else if( mRefMult_corr>CentBin9_vzdep[vzid][i] &&
1764  mRefMult_corr<CentBin9_vzdep[vzid][i+1] ) iCent = i;
1765  }
1766  // 80-100% for icent=-1
1767  if( mRefMult_corr>0 && mRefMult_corr<CentBin9_vzdep[vzid][0] ) iCent = -1;
1768  return ( iCent==9999 ) ? -1 : iCent;
1769 }
1770 
1771 //_________________
1772 Int_t StRefMultCorr::getCentralityBin16VzDep() const {
1773  const Int_t vzid = getVzWindowForVzDepCentDef();
1774  Int_t iCent = 9999;
1775  for(Int_t i=0; i<16; i++) {
1776  if (i == 15) {
1777  if (mRefMult_corr > CentBin16_vzdep[vzid][i] && mRefMult_corr < 50000) {
1778  iCent = i;
1779  }
1780  }
1781  else if (mRefMult_corr > CentBin16_vzdep[vzid][i] &&
1782  mRefMult_corr < CentBin16_vzdep[vzid][i + 1]) {
1783  iCent = i;
1784  }
1785  }
1786  // 80-100% for icent=-1
1787  if( mRefMult_corr>0 && mRefMult_corr<CentBin16_vzdep[vzid][0] ) {
1788  iCent = -1;
1789  }
1790  return (iCent==9999) ? -1 : iCent;
1791 }
1792 
1793 //_________________
1794 const Int_t StRefMultCorr::getRefX() const {
1795  if ( mName.CompareTo("grefmult", TString::kIgnoreCase) == 0 ) return 0;
1796  else if ( mName.CompareTo("refmult", TString::kIgnoreCase) == 0 ) return 1;
1797  else if ( mName.CompareTo("refmult2", TString::kIgnoreCase) == 0 ) return 2;
1798  else if ( mName.CompareTo("refmult3", TString::kIgnoreCase) == 0 ) return 3;
1799  else if ( mName.CompareTo("refmult4", TString::kIgnoreCase) == 0 ) return 4;
1800  else if ( mName.CompareTo("fxtmult", TString::kIgnoreCase) == 0 ) return 5;
1801  // else if ( mName.CompareTo("refmult6", TString::kIgnoreCase) == 0 ) return 6;
1802  else return 9999;
1803 }
1804 
1805 //_________________
1806 const Int_t StRefMultCorr::getNumberOfDatasets() const {
1807  if ( mName.CompareTo("grefmult", TString::kIgnoreCase) == 0 ) return nID_gref;
1808  else if ( mName.CompareTo("refmult", TString::kIgnoreCase) == 0 ) return nID_ref1;
1809  else if ( mName.CompareTo("refmult2", TString::kIgnoreCase) == 0 ) return nID_ref2;
1810  else if ( mName.CompareTo("refmult3", TString::kIgnoreCase) == 0 ) return nID_ref3;
1811  else if ( mName.CompareTo("refmult4", TString::kIgnoreCase) == 0 ) return nID_ref4;
1812  else if ( mName.CompareTo("fxtmult", TString::kIgnoreCase) == 0 ) return nID_ref5;
1813  // else if ( mName.CompareTo("refmult6", TString::kIgnoreCase) == 0 ) return nID_ref6;
1814  else return 9999;
1815 }
1816 
1817 //_________________
1818 void StRefMultCorr::readHeaderFile() {
1819 
1820  //std::vector<std::string> sParam_ShapeWeight = StringSplit(getParamX_ShapeWeight(1,1),',');
1821  //for(int ib=0;ib<sParam_ShapeWeight.size(); ib++) std::cout<<"sParam_ShapeWeight[i]: "<<sParam_ShapeWeight[ib]<<std::endl;
1822 
1823  // Read reference multiplicity to work with (grefmult, refMult, refMult2, refMult3, fxtMult)
1824  const Int_t refX = getRefX();
1825  mRefX = refX;
1826  // Retrieve number of datasets for the given mRefX
1827  const Int_t nID = getNumberOfDatasets();
1828 
1829  for(int iID=0; iID<nID; iID++) {
1830 
1831  //
1832  // Year, energy, run numbers, Vz cut
1833  //
1834  Int_t year; Double_t energy;
1835  std::vector<std::string> sParam = StringSplit(getParamX(refX,iID,0),':');
1836  year = stoi(sParam[0]);
1837  energy = stoi(sParam[1]);
1838  std::vector<std::string> sRuns = StringSplit(sParam[2],',');
1839 
1840  Int_t startRunId=0, stopRunId=0;
1841  startRunId = stoi(sRuns[0]);
1842  stopRunId = stoi(sRuns[1]);
1843 
1844  Double_t startZvertex=-9999., stopZvertex=-9999. ;
1845  std::vector<std::string> sVz = StringSplit(sParam[3],',');
1846  startZvertex = stod(sVz[0]);
1847  stopZvertex = stod(sVz[1]);
1848 
1849  mYear.push_back(year);
1850  mBeginRun.insert(std::make_pair(std::make_pair(energy, year), startRunId));
1851  mEndRun.insert( std::make_pair(std::make_pair(energy, year), stopRunId));
1852  mStart_runId.push_back( startRunId ) ;
1853  mStop_runId.push_back( stopRunId ) ;
1854  mStart_zvertex.push_back( startZvertex ) ;
1855  mStop_zvertex.push_back( stopZvertex ) ;
1856 
1857  //
1858  // Centrality definition
1859  //
1860  std::vector<std::string> sParamCent = StringSplit(getParamX(refX,iID,1),',');
1861  for(UInt_t i=0; i<sParamCent.size(); i++) {
1862  mCentrality_bins[i].push_back( stoi(sParamCent[i]) );
1863  }
1864  mCentrality_bins[mNCentrality].push_back( 5000 );
1865 
1866  //
1867  // Normalization range
1868  //
1869  Double_t normalize_stop=-1.0 ;
1870  normalize_stop = stod(getParamX(refX,iID,2)) ;
1871  mNormalize_stop.push_back( normalize_stop );
1872 
1873  //
1874  // Acceptance (vz) correction
1875  //
1876  std::vector<std::string> sParamVz = StringSplit(getParamX(refX,iID,3),',');
1877  for(UInt_t i=0; i<mNPar_z_vertex; i++) {
1878  Double_t val = -9999.;
1879  if(i<sParamVz.size()) val = stod(sParamVz[i]);
1880  else val = 0.0;
1881  mPar_z_vertex[i].push_back( val );
1882  }
1883 
1884  //
1885  // Trigger inefficiency correction
1886  //
1887  std::vector<std::string> sParamTrig = StringSplit(getParamX(refX,iID,4),',');
1888  for(UInt_t i=0; i<mNPar_weight; i++) {
1889  Double_t val = -9999.;
1890  if(i<sParamTrig.size()) val = stod(sParamTrig[i]);
1891  else val = 0.0;
1892  mPar_weight[i].push_back( val );
1893  }
1894 
1895  //
1896  // Luminosity correction
1897  //
1898  std::vector<std::string> sParamLumi = StringSplit(getParamX(refX,iID,5),',');
1899  for(UInt_t i=0; i<mNPar_luminosity; i++) {
1900  Double_t val = -9999.;
1901  if(i<sParamLumi.size()) val = stod(sParamLumi[i]);
1902  else val = 0.0;
1903  mPar_luminosity[i].push_back( val );
1904  }
1905 
1906  // std::cout << refX <<" "<< iID <<"/"<< nID << std::endl;
1907  }
1908 
1909  std::cout << "StRefMultCorr::readHeaderFile [" << mName
1910  << "] Correction parameters and centrality definitions have been successfully read for refX=" << refX
1911  << std::endl;
1912 }
1913 
1914 //_________________
1915 void StRefMultCorr::readBadRunsFromHeaderFile() {
1916 
1917  //
1918  // TODO: Modify to read only bad runs associated with the requested year
1919  //
1920 
1921  for(Int_t i=0; i<nBadRun_refmult_2010; i++) {
1922  mBadRun.push_back(badrun_refmult_2010[i]);
1923  }
1924  std::cout << "read in nBadRun_refmult_2010: " << nBadRun_refmult_2010 << std::endl;
1925 
1926  for(Int_t i=0; i<nBadRun_refmult_2011; i++) {
1927  mBadRun.push_back(badrun_refmult_2011[i]);
1928  }
1929  std::cout << "read in nBadRun_refmult_2011: " << nBadRun_refmult_2011 << std::endl;
1930 
1931  for(Int_t i=0; i<nBadRun_grefmult_2014; i++) {
1932  mBadRun.push_back(badrun_grefmult_2014[i]);
1933  }
1934  std::cout << "read in nBadRun_grefmult_2014: " << nBadRun_grefmult_2014 << std::endl;
1935 
1936  for(Int_t i=0; i<nBadRun_grefmult_2016; i++) {
1937  mBadRun.push_back(badrun_grefmult_2016[i]);
1938  }
1939  std::cout << "read in nBadRun_grefmult_2016: " << nBadRun_grefmult_2016 << std::endl;
1940 
1941  for(Int_t i=0; i<nBadRun_refmult_2017; i++) {
1942  mBadRun.push_back(badrun_refmult_2017[i]);
1943  }
1944  std::cout << "read in nBadRun_refmult_2017: " << nBadRun_refmult_2017 << std::endl;
1945 
1946  for (Int_t i = 0; i < nBadRun_refmult_2018; i++) {
1947  mBadRun.push_back(badrun_refmult_2018[i]);
1948  }
1949  std::cout << "read in nBadRun_refmult_2018: " << nBadRun_refmult_2018 << std::endl;
1950 
1951  for (Int_t i = 0; i < nBadRun_refmult_2019; i++) {
1952  mBadRun.push_back(badrun_refmult_2019[i]);
1953  }
1954  std::cout << "read in nBadRun_refmult_2019: " << nBadRun_refmult_2019 << std::endl;
1955 
1956  for (Int_t i = 0; i < nBadRun_refmult_2020; i++) {
1957  mBadRun.push_back(badrun_refmult_2020[i]);
1958  }
1959  std::cout << "read in nBadRun_refmult_2020: " << nBadRun_refmult_2020 << std::endl;
1960 
1961  for (Int_t i = 0; i < nBadRun_refmult_2021; i++) {
1962  mBadRun.push_back(badrun_refmult_2021[i]);
1963  }
1964  std::cout << "read in nBadRun_refmult_2021: " << nBadRun_refmult_2021 << std::endl;
1965 
1967  if ( mName.CompareTo("grefmult", TString::kIgnoreCase) == 0 ) {
1968  std::cout << "StRefMultCorr::readBadRunsFromHeaderFile Bad runs for year 2010, 2011, 2017, 2018 and 2019 have been read." << std::endl;
1969  }
1970 }
1971 
1972 //_________________
1973 void StRefMultCorr::print(const Option_t* option) const {
1974  std::cout << "StRefMultCorr::print Print input parameters for "
1975  << mName << " ========================================" << std::endl << std::endl;
1976  // Option switched off, can be used to specify parameters
1977  // const TString opt(option);
1978 
1979  // Int_t input_counter = 0;
1980  for(UInt_t id=0; id<mStart_runId.size(); id++) {
1981  //std::cout << "Data line = " << input_counter << ", Start_runId = " << Start_runId[input_counter] << ", Stop_runId = " << Stop_runId[input_counter] << std::endl;
1982  // const UInt_t id = mStart_runId.size()-1;
1983 
1984  // Removed line break
1985  std::cout << " Index=" << id;
1986  std::cout << Form(" Run=[%8d, %8d]", mStart_runId[id], mStop_runId[id]);
1987  std::cout << Form(" z-vertex=[%1.1f, %1.1f]", mStart_zvertex[id], mStop_zvertex[id]);
1988  std::cout << ", Normalize_stop=" << mNormalize_stop[id];
1989  std::cout << std::endl;
1990 
1991  // if(opt.IsWhitespace()){
1992  // continue ;
1993  // }
1994 
1995  std::cout << "Centrality: ";
1996 
1997  for(Int_t i=0;i<mNCentrality;i++) {
1998  std::cout << Form(" >%2d%%", 80-5*i);
1999  }
2000  std::cout << std::endl;
2001  std::cout << "RefMult: ";
2002 
2003  for(Int_t i=0;i<mNCentrality;i++) {
2004  // std::cout << Form("StRefMultCorr::read Centrality %3d-%3d %%, refmult > %4d", 75-5*i, 80-5*i, mCentrality_bins[i][id]) << std::endl;
2005  const TString tmp(">");
2006  const TString centrality = tmp + Form("%d", mCentrality_bins[i][id]);
2007  std::cout << Form("%6s", centrality.Data());
2008  }
2009  std::cout << std::endl;
2010 
2011  for(Int_t i=0;i<mNPar_z_vertex;i++) {
2012  std::cout << " mPar_z_vertex[" << i << "] = " << mPar_z_vertex[i][id];
2013  }
2014  std::cout << std::endl;
2015 
2016  for(Int_t i=0;i<mNPar_weight;i++) {
2017  std::cout << " mPar_weight[" << i << "] = " << mPar_weight[i][id];
2018  }
2019  std::cout << std::endl;
2020 
2021  for(Int_t i=0;i<mNPar_luminosity;i++) {
2022  std::cout << " mPar_luminosity[" << i << "] = " << mPar_luminosity[i][id];
2023  }
2024  std::cout << std::endl << std::endl;
2025  }
2026  std::cout << "=====================================================================================" << std::endl;
2027 }
2028 
2029 //_________________
2030 std::vector<std::string> StRefMultCorr::StringSplit( const std::string str, const char sep ) const {
2031  std::vector<std::string> vstr;
2032  std::stringstream ss(str);
2033  std::string buffer;
2034  while( getline(ss,buffer,sep) ) {
2035  vstr.push_back(buffer);
2036  }
2037  return vstr;
2038 }
2039 
2040 //________________
2041 Double_t StRefMultCorr::calcPileUpRefMult(Double_t ntofmatch, Double_t x0, Double_t x1,
2042  Double_t x2, Double_t x3, Double_t x4) const {
2043  return ( x0 + x1*(ntofmatch) + x2*pow(ntofmatch,2) + x3*pow(ntofmatch,3) + x4*pow(ntofmatch,4) );
2044 }
Double_t luminosityCorrection(Double_t zdcCoincidenceRate) const
Luminosity correction factor.
Definition: FJcore.h:367
Int_t getCentralityBin9() const
Get 9 centrality bins (10% increment except for 0-5 and 5-10)
Double_t getShapeWeight_SubVz2Center() const
Shape reweighting of refmult: ratio of refMult in each Vz bin to that in the center (|Vz|&lt;10cm) ...
Int_t getCentralityBin16() const
Get 16 centrality bins (5% increment, 0-5, 5-10, ..., 75-80)
Int_t getBeginRun(const Double_t energy, const Int_t year)
Return the first runId from energy and year.
Double_t sampleRefMult(Int_t refMult) const
Sample refMult -&gt; convert integer to double.
Double_t getRefMultCorr() const
Get corrected multiplicity, correction as a function of primary z-vertex.
void print(const Option_t *option="") const
Print all parameters.
Double_t getWeight() const
Total weighting factor: incorporates shape and trigger efficiency weights.
Double_t triggerWeight() const
Trigger efficiency: fit of the Glauber/Data.
Bool_t passnTofMatchRefmultCut(Double_t refmult, Double_t ntofmatch, Double_t vz=0.) const
Check if NOT pile-up event.
Double_t vzCorrection(Double_t z) const
Vz correction factor.
Bool_t isBadRun(const Int_t RunId)
Check if run is bad.
Int_t getEndRun(const Double_t energy, const Int_t year)
Return the last runId from energy and year.
StRefMultCorr(const TString name="refmult", const TString subname="Def", const TString libname="Def")
virtual ~StRefMultCorr()
Destructor.
void init(const Int_t RunId)
Initialization of centrality bins etc.