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
48 StRefMultCorr::~StRefMultCorr() {
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 if (mParameterIndex == 9) { // Run 20 Au+Au 9.8 GeV (sqrt(s_NN)=4.5 GeV) TriggerID = 740000
745  b0=24.7124572851806;
746  b1=5.71868597816542;
747  b2=-0.00312345552143256;
748  b3=-0.00017256714660496;
749  b4=-6.34006022264486e-07;
750  c0=-32.7658746560952;
751  c1=2.70684190541053;
752  c2=0.00128464121679533;
753  c3=0.000470785124435377;
754  c4=-4.47307728635742e-06;
755 
756  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
757  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
758  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
759  }
760  else if (mParameterIndex == 10) { // Run 20 Au+Au 9.8 GeV (sqrt(s_NN)=4.5 GeV) TriggerID = 740010
761  b0=35.02;
762  b1=3.586;
763  b2=0.1368;
764  b3=-2.578e-3;
765  b4=1.189e-5;
766  c0=-24.84;
767  c1=3.289;
768  c2=-0.05722;
769  c3=1.491e-3;
770  c4=-9.678e-6;
771 
772  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
773  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
774  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
775  if (ntofmatch >= 95 && refmult > 367) notPileUp = kFALSE;
776  }
777  else {
778  notPileUp = kTRUE;
779  }
780  } // else if ( mRefX == 5 ) { // fxtMult
781  else if ( mRefX == 6 ) { // refMult6
782  if ( mParameterIndex==0 ) { // O+O 200 GeV 2021
783  b0 = 10.805229670974072;
784  b1 = 4.222122255649001;
785  b2 = -0.059023457533258925;
786  b3 = 0.0007297279166792341;
787  b4 = -3.730248246986408e-06;
788  c0 = -22.14994423125864;
789  c1 = 2.295827717755825;
790  c2 = -0.007304309390513328;
791  c3 = 9.397986288422607e-05;
792  c4 = -9.080150075680224e-07;
793  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
794  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
795  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
796  } // if ( mParameterIndex==0 ) { // O+O 200 GeV 2021
797  else if ( mParameterIndex==1 ) { // d+Au 200 GeV 2021
798  b0 = 2.4412914662443033;
799  b1 = 5.523540420923605;
800  b2 = -0.16458436958697667;
801  b3 = 0.002805908341435613;
802  b4 = -1.6300934820294975e-05;
803  c0 = -0.86595124167792;
804  c1 = 0.44263208748354943;
805  c2 = 0.06024976895762696;
806  c3 = -0.0013523620327006189;
807  c4 = 1.0553696607739253e-05;
808  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
809  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
810  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
811  } // else if ( mParameterIndex==1 ) { // d+Au 200 GeV 2021
812  else {
813  notPileUp = kTRUE;
814  }
815  } // else if ( mRefX == 6 ) { // refMult6
816  else if ( mRefX == 7 ) { // TotnMIP
817  if ( mParameterIndex==0 ) { // O+O 200 GeV 2021
818  b0 = 10.805229670974072;
819  b1 = 4.222122255649001;
820  b2 = -0.059023457533258925;
821  b3 = 0.0007297279166792341;
822  b4 = -3.730248246986408e-06;
823  c0 = -22.14994423125864;
824  c1 = 2.295827717755825;
825  c2 = -0.007304309390513328;
826  c3 = 9.397986288422607e-05;
827  c4 = -9.080150075680224e-07;
828  refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
829  refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
830  notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
831  } // if ( mParameterIndex==0 ) { // O+O 200 GeV 2021
832  else {
833  notPileUp = kTRUE;
834  }
835  } // else if ( mRefX == 7 ) { // TotnMIP
836 
837 
838  if (mVerbose) {
839  std::cout << "\t notPileUp: ";
840  if (notPileUp) {
841  std::cout << "TRUE";
842  }
843  else {
844  std::cout << "FALSE";
845  }
846  std::cout << "\t[DONE]\n";
847  }
848 
849  return notPileUp;
850 }
851 
852 //_________________
853 Bool_t StRefMultCorr::passnTofMatchTotnMIPCut(Double_t totnMIP, Double_t ntofmatch,
854  Double_t vz /* default = 0*/) const {
855 
856 
857  if (mVerbose) {
858  std::cout << "Pile up check using nBTOFMatch vs. totnMIP...";
859  }
860  Double_t b0{}, b1{}, b2{}, b3{}, b4{};
861  Double_t c0{}, c1{}, c2{}, c3{}, c4{};
862  Double_t totnMIPcutmax{};
863  Double_t totnMIPcutmin{};
864 
865  Bool_t notPileUp = kFALSE;
866 
867  if ( mRefX == 6 ) { // refMult6
868  if ( mParameterIndex==0 ) { // O+O 200 GeV 2021
869  b0 = 49.92964577941192;
870  b1 = 30.084027095279428;
871  b2 = -0.6414804471204509;
872  b3 = 0.006675174653594674;
873  b4 = -2.690799009087484e-05;
874  c0 = -114.51733250850462;
875  c1 = 2.2552394816896664;
876  c2 = 0.17302732482464722;
877  c3 = -0.002572787709722221;
878  c4 = 9.276963258195168e-06;
879  totnMIPcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
880  totnMIPcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
881  notPileUp = isInPileUpRefMultLimits(totnMIP, totnMIPcutmin, totnMIPcutmax);
882  } // if ( mParameterIndex==0 ) { // O+O 200 GeV 2021
883  else if ( mParameterIndex==1 ) { // d+Au 200 GeV 2021
884  b0 = 13.854308990154239;
885  b1 = 48.38889240643103;
886  b2 = -1.820427172052759;
887  b3 = 0.0313802142314363;
888  b4 = -0.00018680122215220316;
889  c0 = 7.022386992461804;
890  c1 = -3.747248801082249;
891  c2 = 0.36559511738864753;
892  c3 = -0.00836993562719528;
893  c4 = 7.050862326468557e-05;
894  totnMIPcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
895  totnMIPcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
896  notPileUp = isInPileUpRefMultLimits(totnMIP, totnMIPcutmin, totnMIPcutmax);
897  } // else if ( mParameterIndex==1 ) { // d+Au 200 GeV 2021
898  else {
899  notPileUp = kTRUE;
900  }
901  } // if ( mRefX == 6 ) { // refMult6
902  else if ( mRefX == 7 ) { // TotnMIP
903  if ( mParameterIndex==0 ) { // O+O 200 GeV 2021
904  b0 = 49.92964577941192;
905  b1 = 30.084027095279428;
906  b2 = -0.6414804471204509;
907  b3 = 0.006675174653594674;
908  b4 = -2.690799009087484e-05;
909  c0 = -114.51733250850462;
910  c1 = 2.2552394816896664;
911  c2 = 0.17302732482464722;
912  c3 = -0.002572787709722221;
913  c4 = 9.276963258195168e-06;
914  totnMIPcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
915  totnMIPcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
916  notPileUp = isInPileUpRefMultLimits(totnMIP, totnMIPcutmin, totnMIPcutmax);
917  } // if ( mParameterIndex==0 ) { // O+O 200 GeV 2021
918  else {
919  notPileUp = kTRUE;
920  }
921  } // else if ( mRefX == 7 ) { // TotnMIP
922 
923  if (mVerbose) {
924  std::cout << "\t notPileUp using : ";
925  if (notPileUp) {
926  std::cout << "TRUE";
927  }
928  else {
929  std::cout << "FALSE";
930  }
931  std::cout << "\t[DONE]\n";
932  }
933 
934  return notPileUp;
935 }
936 
937 //_________________
938 Bool_t StRefMultCorr::isIndexOk() const {
939  // mParameterIndex not initialized (-1)
940  if ( mParameterIndex == -1 ) {
941  Error("StRefMultCorr::isIndexOk", "mParameterIndex = -1. Call init(const Int_t RunId) function to initialize centrality bins, corrections");
942  Error("StRefMultCorr::isIndexOk", "mParameterIndex = -1. or use valid run numbers defined in Centrality_def_%s.txt", mName.Data());
943  Error("StRefMultCorr::isIndexOk", "mParameterIndex = -1. exit");
944  std::cout << std::endl;
945  // Stop the process if invalid run number found
946  // exit(0);
947  return kFALSE;
948  }
949 
950  // Out of bounds
951  if ( mParameterIndex >= (Int_t)mStart_runId.size() ) {
952  Error("StRefMultCorr::isIndexOk",
953  Form("mParameterIndex = %d > max number of parameter set = %d. Make sure you put correct index for this energy",
954  mParameterIndex, mStart_runId.size()));
955  return kFALSE ;
956  }
957 
958  return kTRUE ;
959 }
960 
961 //_________________
962 Bool_t StRefMultCorr::isZvertexOk() const {
963  // Primary z-vertex check
964  return (mVz > mStart_zvertex[mParameterIndex] &&
965  mVz < mStop_zvertex[mParameterIndex]);
966 }
967 
968 //_________________
969 Bool_t StRefMultCorr::isRefMultOk() const {
970  // Invalid index
971  if ( !isIndexOk() ) return kFALSE ;
972 
973  // select 0-80%
974  return (mRefMult_corr > mCentrality_bins[0][mParameterIndex] &&
975  mRefMult_corr < mCentrality_bins[mNCentrality][mParameterIndex]);
976 }
977 
978 //_________________
979 Bool_t StRefMultCorr::isCentralityOk(const Int_t icent) const {
980  // Invalid centrality id
981  if ( icent < -1 || icent >= mNCentrality+1 ) return kFALSE ;
982 
983  // Invalid index
984  if ( !isIndexOk() ) return kFALSE ;
985 
986  // Special case
987  // 1. 80-100% for icent=-1
988  if ( icent == -1 ) return ( mRefMult_corr <= mCentrality_bins[0][mParameterIndex] );
989 
990  // 2. icent = mNCentrality
991  if ( icent == mNCentrality ) return ( mRefMult_corr <= mCentrality_bins[mNCentrality][mParameterIndex] );
992 
993  const Bool_t isOK = ( mRefMult_corr > mCentrality_bins[icent][mParameterIndex] &&
994  mRefMult_corr <= mCentrality_bins[icent+1][mParameterIndex] );
995  if(mVerbose && isOK) {
996  std::cout << "StRefMultCorr::isCentralityOk refmultcorr = " << mRefMult_corr
997  << " min. bin = " << mCentrality_bins[icent][mParameterIndex]
998  << " max. bin = " << mCentrality_bins[icent+1][mParameterIndex]
999  << std::endl;
1000  }
1001  return isOK ;
1002 }
1003 
1004 //_________________
1005 void StRefMultCorr::init(const Int_t RunId) {
1006  // Reset mParameterIndex
1007  mParameterIndex = -1 ;
1008 
1009  // call setParameterIndex
1010  setParameterIndex(RunId) ;
1011 }
1012 
1013 //_________________
1014 Int_t StRefMultCorr::setParameterIndex(const Int_t RunId) {
1015  // Determine the corresponding parameter set for the input RunId
1016  for(UInt_t npar = 0; npar < mStart_runId.size(); npar++) {
1017  if(RunId >= mStart_runId[npar] && RunId <= mStop_runId[npar]) {
1018  mParameterIndex = npar ;
1019  // std::cout << "StRefMultCorr::setParameterIndex Parameter set = " << mParameterIndex << " for RUN " << RunId << std::endl;
1020  break ;
1021  }
1022  }
1023 
1024  // Multiple parameters/definitions for Run14/16 data
1025  // Set mParameterIndex by hand
1026  // For Run14 P16id production
1027  // For Run16 P16ij production
1028  if ( mName.CompareTo("grefmult", TString::kIgnoreCase) == 0 ) {
1029  if ( mSubName.CompareTo("Run14_AuAu200_VpdMB5", TString::kIgnoreCase) == 0 ) {
1030  if(RunId/1000000==15 && mLibName.CompareTo("P16id", TString::kIgnoreCase) == 0 ) {
1031  mParameterIndex = 0;
1032  if (mVzEdgeForWeight.empty()) {
1033  setVzForWeight(nWeightVzBin_Run14_P16id,
1034  WeightVzEdgeMin_Run14_P16id,
1035  WeightVzEdgeMax_Run14_P16id);
1036  }
1037  if (mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
1038  readScaleForWeight(nWeightgRefmultBin_Run14_P16id,
1039  weight_VpdMB5ToVpdMB30_Run14_P16id);
1040  }
1041  } // if(RunId/1000000==15 && mLibName.CompareTo("P16id", TString::kIgnoreCase) == 0 )
1042  } // if ( mSubName.CompareTo("Run14_AuAu200_VpdMB5", TString::kIgnoreCase) == 0 )
1043  else if ( mSubName.CompareTo("Run16_AuAu200_VpdMB5", TString::kIgnoreCase) == 0 ) {
1044  if(mLibName.CompareTo("P16ij", TString::kIgnoreCase) == 0) {
1045  // prod.1
1046  if (RunId / 1000 >= 17039 && RunId / 1000 <= 17130) {
1047  mParameterIndex = 4;
1048  if (mVzEdgeForWeight.empty()) {
1049  setVzForWeight(nWeightVzBin_Run16_P16ij_prod1,
1050  WeightVzEdgeMin_Run16_P16ij_prod1,
1051  WeightVzEdgeMax_Run16_P16ij_prod1);
1052  }
1053  if (mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
1054  readScaleForWeight(nWeightgRefmultBin_Run16_P16ij_prod1,
1055  weight_VpdMB5ToVpdMBnoVtx_Run16_P16ij_prod1);
1056  }
1057  } // if (RunId / 1000 >= 17039 && RunId / 1000 <= 17130)
1058  // prod.2
1059  else if (RunId / 1000 >= 17169 && RunId / 1000 <= 17179) {
1060  mParameterIndex = 5;
1061  if (mVzEdgeForWeight.empty()) {
1062  setVzForWeight(nWeightVzBin_Run16_P16ij_prod2,
1063  WeightVzEdgeMin_Run16_P16ij_prod2,
1064  WeightVzEdgeMax_Run16_P16ij_prod2);
1065  }
1066  if (mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
1067  readScaleForWeight(nWeightgRefmultBin_Run16_P16ij_prod2,
1068  weight_VpdMB5ToVpdMBnoVtx_Run16_P16ij_prod2);
1069  }
1070  } // else if (RunId / 1000 >= 17169 && RunId / 1000 <= 17179)
1071  } // if(mLibName.CompareTo("P16ij", TString::kIgnoreCase) == 0)
1072  } // else if ( mSubName.CompareTo("Run16_AuAu200_VpdMB5", TString::kIgnoreCase) == 0 )
1073  else if ( mSubName.CompareTo("Run14_AuAu200_VpdMB30", TString::kIgnoreCase) == 0 ) {
1074  if (mLibName.CompareTo("P16id", TString::kIgnoreCase) == 0) {
1075  mParameterIndex = 1;
1076  }
1077  }
1078  else if ( mSubName.CompareTo("Run14_AuAu200_VpdMBnoVtx_LowMid", TString::kIgnoreCase) == 0 ) {
1079  if(mLibName.CompareTo("P16id", TString::kIgnoreCase) == 0) {
1080  mParameterIndex = 2;
1081  }
1082  }
1083  else if ( mSubName.CompareTo("Run14_AuAu200_VpdMBnoVtx_High", TString::kIgnoreCase) == 0 ) {
1084  if(mLibName.CompareTo("P15ic", TString::kIgnoreCase) == 0) {
1085  mParameterIndex = 3;
1086  }
1087  }
1088  else if ( mSubName.CompareTo("Run16_AuAu200_VpdMBnoVtx", TString::kIgnoreCase) == 0 ) {
1089  if(mLibName.CompareTo("P16ij", TString::kIgnoreCase) == 0) {
1090  mParameterIndex = 6;
1091  }
1092  }
1093  else{
1094  mParameterIndex = -1;
1095  }
1096  }
1097 
1098  // Run numbers for isobar data are not successive
1099  if ( mName.CompareTo("refmult", TString::kIgnoreCase) == 0 &&
1100  mSubName.CompareTo("Isobar", TString::kIgnoreCase) == 0 ) {
1101 
1102  // std::cout <<"This is isobaric runs"<< std::endl;
1103  mIsZr = mIsRu = kFALSE;
1104  for (Int_t i = 0; i < nRunId_Zr; i++) {
1105  if (RunId == IsobarRunId_Zr[i]) {
1106  mIsZr = kTRUE;
1107  mParameterIndex = 36;
1108  // std::cout <<"--> Zr+Zr"<< std::endl;
1109  break;
1110  }
1111  }
1112  for (Int_t i = 0; i < nRunId_Ru; i++) {
1113  if (RunId == IsobarRunId_Ru[i]) {
1114  mIsRu = kTRUE;
1115  mParameterIndex = 37;
1116  // std::cout <<"--> Ru+Ru"<< std::endl;
1117  break;
1118  }
1119  }
1120  if(!mIsZr && !mIsRu) {
1121  Error("StRefMultCorr::setParameterIndex", "RUN %d is not isobaric data", RunId);
1122  }
1123  }
1124 
1125  // std::cout <<"mParameterIndex = "<< mParameterIndex << std::endl;
1126 
1127  if(mParameterIndex == -1) {
1128  Error("StRefMultCorr::setParameterIndex", "Parameter set does not exist for RUN %d", RunId);
1129  }
1130  //else std::cout << "Parameter set = " << npar_set << std::endl;
1131 
1132  if (mVerbose) {
1133  std::cout << "Parameter index set to: " << mParameterIndex << std::endl;
1134  }
1135 
1136  return mParameterIndex ;
1137 }
1138 
1139 //_________________
1140 Double_t StRefMultCorr::getRefMultCorr() const {
1141  // Call initEvent() first
1142  return mRefMult_corr ;
1143 }
1144 
1145 //________________
1146 Double_t StRefMultCorr::luminosityCorrection(Double_t zdcCoincidenceRate) const {
1147 
1148  Double_t lumiCorr = 1.;
1149 
1150  if ( mVerbose ) {
1151  std::cout << "Estimation of luminosity correction factor..." << std::endl;
1152  std::cout << "\t ZDC coincidence rate: " << zdcCoincidenceRate << " BBC coincidence rate: " << 0 << std::endl;
1153  }
1154  // 200 GeV only. correction = 1 for all the other energies for BES-I
1155  // the above statement may not true for BES-II, since the luminosity is much higher than BES-I, add by Zaochen
1156  // better to check the <Refmult> vs ZDCX to see whether they are flat or not, add by Zaochen
1157  const Double_t par0_lum = mPar_luminosity[0][mParameterIndex] ;
1158  const Double_t par1_lum = mPar_luminosity[1][mParameterIndex] ;
1159 
1160  if( mParameterIndex==36 || mParameterIndex==37 || mParameterIndex==40 ||
1161  ( mRefX==5 && mParameterIndex==7 ) ||
1162  ( mRefX==5 && mParameterIndex==8 ) ) {
1163  // if(mYear[mParameterIndex] == 2018 && mIsZr) zdcmean = 96.9914;
1164  // if(mYear[mParameterIndex] == 2018 && mIsRu) zdcmean = 97.9927;
1165  Double_t b_prime = 1.;
1166  if(mParameterIndex==36) b_prime = 96.9914; // Zr
1167  if(mParameterIndex==37) b_prime = 97.9927; // Ru
1168  if(mParameterIndex==40) b_prime = 213.383; // AuAu 200GeV Run19
1169  if(mParameterIndex==7 ) b_prime = 106.245; // AuAu 5.2GeV FXT Run20
1170  if(mParameterIndex==8 ) b_prime = 114.041; // AuAu 6.2GeV FXT Run20
1171  lumiCorr = (par0_lum<std::numeric_limits<double>::epsilon() ) ? 1.0 : b_prime/(par0_lum+zdcCoincidenceRate*par1_lum);
1172  }
1173  else {
1174  lumiCorr = (par0_lum < std::numeric_limits<double>::epsilon() ) ? 1.0 : 1.0/(1.0 + par1_lum/par0_lum*zdcCoincidenceRate/1000.);
1175  }
1176 
1177  // from Run14, P16id, for VpdMB5/VPDMB30/VPDMB-noVtx, use refMult at ZdcX=30, other is at ZdcX=0;
1178  // -->changed by xlchen@lbl.gov, Run16 ~ 50kHz
1179  if(
1180  ( mSubName.CompareTo("Run14_AuAu200_VpdMB5", TString::kIgnoreCase) == 0 && mLibName.CompareTo("P16id", TString::kIgnoreCase) == 0 )
1181  || ( mSubName.CompareTo("Run14_AuAu200_VpdMB30", TString::kIgnoreCase) == 0 && mLibName.CompareTo("P16id", TString::kIgnoreCase) == 0 )
1182  || ( mSubName.CompareTo("Run14_AuAu200_VpdMBnoVtx_LowMid", TString::kIgnoreCase) == 0 && mLibName.CompareTo("P16id", TString::kIgnoreCase) == 0 )
1183  || ( mSubName.CompareTo("Run14_AuAu200_VpdMBnoVtx_High", TString::kIgnoreCase) == 0 && mLibName.CompareTo("P15ic", TString::kIgnoreCase) == 0 )
1184  || ( mSubName.CompareTo("Run16_AuAu200_VpdMB5", TString::kIgnoreCase) == 0 && mLibName.CompareTo("P16ij", TString::kIgnoreCase) == 0 )
1185  || ( mSubName.CompareTo("Run16_AuAu200_VpdMBnoVtx", TString::kIgnoreCase) == 0 && mLibName.CompareTo("P16ij", TString::kIgnoreCase) == 0 )
1186  ) {
1187  float zdcmean = 0;
1188  if(mYear[mParameterIndex] == 2014) zdcmean = 30.;
1189  if(mYear[mParameterIndex] == 2016) zdcmean = 50.;
1190  lumiCorr = (par0_lum==0.0) ? lumiCorr : lumiCorr*(par0_lum+par1_lum*zdcmean)/par0_lum;
1191  }
1192 
1193  if (mVerbose) {
1194  std::cout << "\tLuminosity correction factor: " << lumiCorr << std::endl;
1195  std::cout << "\t[DONE]\n";
1196  }
1197  return lumiCorr;
1198 }
1199 
1200 //________________
1201 Double_t StRefMultCorr::vzCorrection(Double_t z) const {
1202 
1203  if ( mVerbose ) {
1204  std::cout << "Estimation of acceptance correction factor..." << std::endl;
1205  std::cout << "\t Vertex z position: " << z << " mParameterIndex: " << mParameterIndex << std::endl;
1206 
1207  }
1208  Double_t vzCorr = 1.;
1209  if ( mParameterIndex < 38 ) {
1210  if (mRefX == 6) { // refMult6
1211  if ( mParameterIndex == 0 ) { // O+O 200 GeV 2021
1212  // New Vz correction. All vz bins bins are normalize to that at the center
1213  vzCorr = oo200_run21_vzCorr_refMult[ getVzWindowForVzDepCentDef() ];
1214  } // if ( mParameterIndex == 0 ) { // O+O 200 GeV 2021
1215 
1216  else if ( mParameterIndex == 1 ) { // d+Au 200 GeV 2021
1217  // New Vz correction. All vz bins bins are normalize to that at the center
1218  vzCorr = dau200_run21_vzCorr_refMult[ getVzWindowForVzDepCentDef() ];
1219  } // else if ( mParameterIndex == 1 ) { // d+Au 200 GeV 2021
1220 
1221  } // if (mRefX == 6)
1222  else if (mRefX == 7) { // totnMIP
1223  if ( mParameterIndex == 0 ) { // O+O 200 GeV 2021
1224  // New Vz correction. All vz bins bins are normalize to that at the center
1225  vzCorr = oo200_run21_vzCorr_totnMIP[ getVzWindowForVzDepCentDef() ];
1226  } // if ( mParameterIndex == 0 ) { // O+O 200 GeV 2021
1227  }
1228  else {
1229  // Old correction based on the 6th-order polynomial fit of the high-end point
1230  // fit of refMult for the given Vz range
1231 
1232  // par0 to par5 define the parameters of a polynomial to parametrize z_vertex dependence of RefMult
1233  const Double_t par0 = mPar_z_vertex[0][mParameterIndex];
1234  const Double_t par1 = mPar_z_vertex[1][mParameterIndex];
1235  const Double_t par2 = mPar_z_vertex[2][mParameterIndex];
1236  const Double_t par3 = mPar_z_vertex[3][mParameterIndex];
1237  const Double_t par4 = mPar_z_vertex[4][mParameterIndex];
1238  const Double_t par5 = mPar_z_vertex[5][mParameterIndex];
1239  const Double_t par6 = mPar_z_vertex[6][mParameterIndex];
1240  // This parameter is usually 0, it takes care for an additional efficiency,
1241  // usually difference between phase A and phase B parameter 0
1242  const Double_t par7 = mPar_z_vertex[7][mParameterIndex];
1243 
1244  const Double_t RefMult_ref = par0; // Reference mean RefMult at z=0
1245  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
1246  if(RefMult_z > 0.0) {
1247  vzCorr = (RefMult_ref + par7)/RefMult_z;
1248  }
1249  }
1250  }
1251  else if ( mParameterIndex == 38 ) { // Au+Au 19.6 GeV Run 19
1252  // New Vz correction. All vz bins bins are normalize to the one at the center
1253  vzCorr = auau19_run19_vzCorr[ getVzWindowForVzDepCentDef() ];
1254  } // else if ( mParameterIndex == 38 )
1255  else if ( mParameterIndex == 39 ) {
1256  // New Vz correction. All vz bins bins are normalize to that at the center
1257  vzCorr = auau14_run19_vzCorr[ getVzWindowForVzDepCentDef() ];
1258  }
1259  else if ( mParameterIndex == 40 ) { // Au+Au 200 GeV Run 19
1260  // New Vz correction. All vz bins bins are normalize to that at the center
1261  vzCorr = auau200_run19_vzCorr[ getVzWindowForVzDepCentDef() ];
1262  }
1263  else if ( mParameterIndex == 41 ) { // Au+Au 7.7 GeV Run 21
1264  // New Vz correction. All vz bins bins are normalize to that at the center
1265  vzCorr = auau7_run21_vzCorr[ getVzWindowForVzDepCentDef() ];
1266  }
1267  else if ( mParameterIndex == 42 ) { // Au+Au 9.2 GeV Run 20 TriggerID = 780020
1268  // New Vz correction. All vz bins bins are normalize to that at the center
1269  vzCorr = auau9_trig2_run20_vzCorr[ getVzWindowForVzDepCentDef() ];
1270  }
1271  else if ( mParameterIndex == 43 ) { // Au+Au 17.3 GeV Run 21
1272  // New Vz correction. All vz bins bins are normalize to that at the center
1273  vzCorr = auau17_run21_vzCorr[ getVzWindowForVzDepCentDef() ];
1274  }
1275  else if ( mParameterIndex == 44 ) { // Au+Au 11.5 GeV Run 20
1276  // New Vz correction. All vz bins bins are normalize to that at the center
1277  vzCorr = auau11_run20_vzCorr[ getVzWindowForVzDepCentDef() ];
1278  }
1279 
1280  if (mVerbose) {
1281  std::cout << "\t Acceptance correction factor: " << vzCorr << std::endl;
1282  std::cout << "\t[DONE]\n";
1283  }
1284  return vzCorr;
1285 }
1286 
1287 //________________
1288 Double_t StRefMultCorr::sampleRefMult(Int_t refMult) const {
1289 
1290  if (mVerbose) {
1291  std::cout << "Sampling refMult value: " << refMult << std::endl;
1292  }
1293 
1294  Double_t refMult_d = -9999.;
1295  if( mParameterIndex>=30 && mParameterIndex<=44 ) {
1296  refMult_d = (Double_t)refMult - 0.5 + gRandom->Rndm();
1297  }
1298  else {
1299  refMult_d = (Double_t)refMult + gRandom->Rndm();
1300  }
1301 
1302  if (mVerbose) {
1303  std::cout << "\tSampled refMult value: " << refMult_d << std::endl;
1304  std::cout << "\t[DONE]\n";
1305  }
1306  return refMult_d;
1307 }
1308 
1309 //________________
1310 Double_t StRefMultCorr::getRefMultCorr(const UShort_t refMult,
1311  const Double_t z,
1312  const Double_t zdcCoincidenceRate,
1313  const UInt_t flag) const {
1314 
1315  if (mVerbose) {
1316  std::cout << "Start refMultCorr calculations" << std::endl
1317  << "Initial values refMult / vz / zdcX / flag: "
1318  << refMult << " / " << z << " / " << zdcCoincidenceRate << " / "
1319  << flag << std::endl;
1320  }
1321  // Apply correction if parameter index & z-vertex are ok
1322  if (!isIndexOk() || !isZvertexOk()) return refMult ;
1323 
1324  // Correction function for RefMult, takes into account z_vertex dependence
1325 
1326  Double_t lumiCorr = luminosityCorrection(zdcCoincidenceRate);
1327  Double_t vzCorr = vzCorrection(z);
1328  Double_t refMult_d = sampleRefMult( refMult );
1329  Double_t refMultCorr = -9999. ;
1330  switch (flag) {
1331  case 0: refMultCorr = refMult_d * lumiCorr; break;
1332  case 1: refMultCorr = refMult_d * vzCorr; break;
1333  case 2: refMultCorr = refMult_d * vzCorr * lumiCorr; break;
1334  default: {
1335  Error("StRefMultCorr::getRefMultCorr", "invalid flag, flag=%d, should be 0, 1 or 2", flag);
1336  refMultCorr = -9999.;
1337  }
1338  } // switch ( flag )
1339 
1340  if (mVerbose) {
1341  std::cout << "Final refMultCorr value: " << refMultCorr << std::endl;
1342  }
1343 
1344  return refMultCorr;
1345 }
1346 
1347 //_________________
1348 void StRefMultCorr::readScaleForWeight(const Char_t* input) {
1349  std::ifstream fin(input) ;
1350  if(!fin) {
1351  Error("StRefMultCorr::readScaleForWeight", "can't open %s", input);
1352  return;
1353  }
1354 
1355  // Users must set the vz bin size by setVzForWeight() (see below)
1356  if(mnVzBinForWeight==0) {
1357  Error("StRefMultCorr::readScaleForWeight",
1358  "Please call setVzForWeight() to set vz bin size");
1359  return;
1360  }
1361 
1362  // Do not allow multiple calls
1363  if(!mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
1364  Error("StRefMultCorr::readScaleForWeight",
1365  "scale factor has already set in the array");
1366  return;
1367  }
1368 
1369  std::cout << "StRefMultCorr::readScaleForWeight Read scale factor ..." << std::flush;
1370 
1371  while(fin) {
1372  Double_t scale[mnVzBinForWeight] ;
1373  for(Int_t i=0; i<mnVzBinForWeight; i++) {
1374  fin >> scale[i] ;
1375  }
1376 
1377  if(fin.eof()) break ;
1378 
1379  for(Int_t i=0; i<mnVzBinForWeight; i++) {
1380  mgRefMultTriggerCorrDiffVzScaleRatio.push_back(scale[i]);
1381  }
1382  }
1383  std::cout << " [OK]" << std::endl;
1384 }
1385 
1386 // NEW version to read Vz dependent weights from header
1387 // Implemented inside StRefMultCorr::setParameterIndex(RunId)
1388 //_________________
1389 void StRefMultCorr::readScaleForWeight(const Int_t nRefmultbin, const Double_t *weight) {
1390 
1391  // Users must set the vz bin size by setVzForWeight() (see below)
1392  if( mnVzBinForWeight==0 ) {
1393  Error("StRefMultCorr::readScaleForWeight",
1394  "Please call setVzForWeight() to set vz bin size");
1395  return;
1396  }
1397 
1398  // Do not allow multiple calls
1399  if(!mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
1400  Error("StRefMultCorr::readScaleForWeight",
1401  "scale factor has already set in the array");
1402  return;
1403  }
1404 
1405  std::cout << "StRefMultCorr::readScaleForWeight Read scale factor ..." << std::flush;
1406 
1407  for(Int_t i=0; i<nRefmultbin*mnVzBinForWeight; i++) {
1408  mgRefMultTriggerCorrDiffVzScaleRatio.push_back(weight[i]);
1409  }
1410 
1411  std::cout << " [OK]" << std::endl;
1412 }
1413 
1414 
1415 // In NEW version, setVzForWeight() is implemented inside StRefMultCorr::setParameterIndex(RunId)
1416 // It does not need to be called by users.
1417 //_________________
1418 void StRefMultCorr::setVzForWeight(const Int_t nbin, const Double_t min, const Double_t max) {
1419  // Do not allow multiple calls
1420  if ( !mVzEdgeForWeight.empty() ) {
1421  Error("StRefMultCorr::setVzForWeight",
1422  "z-vertex range for weight has already been defined");
1423  return;
1424  }
1425 
1426  mnVzBinForWeight = nbin ;
1427  // calculate increment size
1428  const Double_t step = (max-min)/(Double_t)nbin;
1429 
1430  for(Int_t i=0; i<mnVzBinForWeight+1; i++) {
1431  mVzEdgeForWeight.push_back( min + step*i );
1432  }
1433 
1434  // Debug
1435  for(Int_t i=0; i<mnVzBinForWeight; i++) {
1436  std::cout << i << " " << step << " " << mVzEdgeForWeight[i]
1437  << ", " << mVzEdgeForWeight[i+1] << std::endl;
1438  }
1439 }
1440 
1441 //_________________
1442 Double_t StRefMultCorr::getScaleForWeight() const {
1443  // Special scale factor for global refmult in Run14 (Run16)
1444  // to account for the relative difference of VPDMB5 w.r.t VPDMB30 (VPDMBnoVtx)
1445 
1446  // return 1 if mgRefMultTriggerCorrDiffVzScaleRatio array is empty
1447  if(mgRefMultTriggerCorrDiffVzScaleRatio.empty()) return 1.0 ;
1448 
1449  // const Int_t nVzBins =6;
1450  // Double_t VzEdge[nVzBins+1]={-6., -4., -2., 0., 2., 4., 6.};
1451 
1452  Double_t VPD5weight=1.0;
1453  for(Int_t j=0;j<mnVzBinForWeight;j++) {
1454  if(mVz>mVzEdgeForWeight[j] && mVz<=mVzEdgeForWeight[j+1]) {
1455  /*
1456  //refMultbin=mgRefMultTriggerCorrDiffVzScaleRatio_2[j]->FindBin(mRefMult_corr+1e-6);
1457  //VPD5weight=mgRefMultTriggerCorrDiffVzScaleRatio_2[j]->GetBinContent(refMultbin);
1458  const Int_t refMultbin=static_cast<Int_t>(mRefMult_corr);
1459  //VPD5weight=mgRefMultTriggerCorrDiffVzScaleRatio[j][refMultbin];
1460  VPD5weight=mgRefMultTriggerCorrDiffVzScaleRatio[refMultbin*mnVzBinForWeight + j];
1461  const Double_t tmpContent=VPD5weight;
1462  if(tmpContent==0 || (mRefMult_corr>500 && tmpContent<=0.65)) VPD5weight=1.15;//Just because the value of the weight is around 1.15
1463  if(mRefMult_corr>500 && tmpContent>=1.35) VPD5weight=1.15;//Remove those Too large weight factor,gRefmult>500
1464  // this weight and reweight should be careful, after reweight(most peripheral),Then weight(whole range)
1465  */
1466 
1467  const Int_t refMultbin=static_cast<Int_t>(mRefMult_corr);
1468  VPD5weight=mgRefMultTriggerCorrDiffVzScaleRatio[refMultbin*mnVzBinForWeight + j];
1469  const Double_t tmpContent=VPD5weight;
1470  // 1) Ratios fluctuate too much at very high gRefmult due to low statistics
1471  // 2) Avoid some events with too high weight
1472  if(mRefMult_corr>550 && (tmpContent>3.0||tmpContent<0.3)) VPD5weight=1.0;
1473  // this weight and reweight should be careful, after reweight(most peripheral),Then weight(whole range)
1474  }
1475  }
1476 
1477  return 1.0/VPD5weight;
1478 }
1479 
1480 //_________________
1481 // For Run18 27 GeV
1482 Double_t StRefMultCorr::getShapeWeight_SubVz2Center() const {
1483 
1484  if (mVerbose) {
1485  std::cout << "Estimating shape correction factor" << std::endl;
1486  }
1487 
1488  Double_t weight = 1.;
1489  Int_t iVzBinIndex = getVzWindowForVzDepCentDef();
1490 
1491  if ( mParameterIndex>=30 && mParameterIndex<=35 ) { // Run18 27 GeV MB
1492 
1493  if(mRefMult_corr>=500.) return 1.0; // almost no Refmult>500 for this collision energy
1494 
1495  Int_t iRunPdIndex = mParameterIndex-30;
1496  Int_t iRefmultBin = (Int_t)(mRefMult_corr/2.); //find the refmult bin matching to the Parameter bin, if binWidth=2
1497  //Int_t iRefmultBin = (Int_t)(mRefMult_corr); //find the refmult bin matching to the Parameter bin, if binWidth=1
1498 
1499  if(iRunPdIndex<0 || iRunPdIndex>5) return 1.0;
1500  if(iVzBinIndex<0 || iVzBinIndex>13) return 1.0;
1501 
1502  //----------------------------------------------------------
1503  //load the ShapeReweight factors in given RunPd and VzBin
1504  //----------------------------------------------------------
1505  std::vector<std::string> sParam_ShapeWeight = StringSplit(getParamX_ShapeWeight(iRunPdIndex,iVzBinIndex),',');
1506  if( iRefmultBin>=(Int_t)sParam_ShapeWeight.size() ) {
1507  std::cout<<"ERROR: sParam_ShapeWeight is out of ranges!!!!!"<<std::endl;
1508  return 1.0;
1509  }
1510 
1511  //std::cout<<"sParam_ShapeWeight.size(): "<<sParam_ShapeWeight.size()<<std::endl;
1512  //for(UInt_t is=0; is<sParam_ShapeWeight.size(); is++) std::cout<<"sParam_ShapeWeight[is]: "<<sParam_ShapeWeight[is]<<std::endl;
1513 
1514  Double_t ShapeReweight = 1.0;
1515 
1516  Double_t tem_ShapeReweight = std::stod( sParam_ShapeWeight[iRefmultBin] );
1517  //prevent the crazy numbers for the large fluctuations
1518  if( tem_ShapeReweight<0.1 ) {
1519  ShapeReweight = 0.1;
1520  }
1521  else if(tem_ShapeReweight>10) {
1522  ShapeReweight = 10.;
1523  }
1524  else if(tem_ShapeReweight>=0.1 && tem_ShapeReweight<=10.) {
1525  ShapeReweight = tem_ShapeReweight;
1526  }
1527 
1528  weight = 1./ShapeReweight;
1529 
1530  } // if ( mParameterIndex>=30 && mParameterIndex<=35 ) { // Run18 27 GeV MB
1531  else if ( mParameterIndex>=36 && mParameterIndex<=37 ) { // Isobar collision 200 GeV 2018
1532 
1533  if (mVz >= -9 && mVz <= 9) {
1534  return 1.;
1535  }
1536 
1537  Int_t mShapeIndex = 0;
1538  if (mIsZr) {
1539  mShapeIndex = 1;
1540  }
1541 
1542  //retrive shape weight
1543  if (iVzBinIndex >= 22) {
1544  weight = ShapeWeightArray[mShapeIndex][iVzBinIndex - 9][TMath::Nint(mRefMult_corr)];
1545  }
1546  else {
1547  weight = ShapeWeightArray[mShapeIndex][iVzBinIndex][TMath::Nint(mRefMult_corr)];
1548  }
1549  //handle bad weight
1550  if (weight == 0 || TMath::IsNaN(weight)) {
1551  weight = 1.;
1552  }
1553  } // else if ( mParameterIndex>=36 && mParameterIndex<=37 ) { // Isobar collision 200 GeV 2018
1554  else if (mParameterIndex == 38) { // Au+Au 19.6 GeV 2019
1555 
1556  if (iVzBinIndex < 0 || iVzBinIndex > auau19_run19_nVzBins) return 1.0;
1557 
1558  weight = auau19_run19_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1559  // Handle bad weight
1560  if (weight == 0 || TMath::IsNaN(weight)) {
1561  weight = 1.;
1562  }
1563  } // else if ( mParameterIndex == 38 ) { // Au+Au 19.6 GeV 2019
1564  else if (mParameterIndex == 39) { // Au+Au 14.6 GeV 2019
1565 
1566  if (iVzBinIndex < 0 || iVzBinIndex > auau14_run19_nVzBins) return 1.0;
1567 
1568  weight = auau14_run19_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1569  // Handle bad weight
1570  if (weight == 0 || TMath::IsNaN(weight)) {
1571  weight = 1.;
1572  }
1573  } // else if ( mParameterIndex == 39 ) { // Au+Au 14.6 GeV 2019
1574  else if (mParameterIndex == 40) { // Au+Au 200 GeV 2019
1575 
1576  if (iVzBinIndex < 0 || iVzBinIndex > auau200_run19_nVzBins) return 1.0;
1577 
1578  weight = auau200_run19_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1579  // Handle bad weight
1580  if (weight == 0 || TMath::IsNaN(weight)) {
1581  weight = 1.;
1582  }
1583  }
1584  else if (mParameterIndex == 41) { // Au+Au 7.7 GeV 2020
1585 
1586  if (iVzBinIndex < 0 || iVzBinIndex > auau7_run21_nVzBins) return 1.0;
1587 
1588  weight = auau7_run21_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1589  // Handle bad weight
1590  if (weight == 0 || TMath::IsNaN(weight)) {
1591  weight = 1.;
1592  }
1593  }
1594  else if (mParameterIndex == 42) { // Au+Au 9.2 GeV 2020 TrigerID = 780020
1595 
1596  if (iVzBinIndex < 0 || iVzBinIndex > auau9_run20_nVzBins) return 1.0;
1597 
1598  weight = auau9_trig2_run20_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1599  // Handle bad weight
1600  if (weight == 0 || TMath::IsNaN(weight)) {
1601  weight = 1.;
1602  }
1603  }
1604  else if (mParameterIndex == 43) { // Au+Au 17.3 GeV 2021
1605 
1606  if (iVzBinIndex < 0 || iVzBinIndex > auau17_run21_nVzBins) return 1.0;
1607 
1608  weight = auau17_run21_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1609  // Handle bad weight
1610  if (weight == 0 || TMath::IsNaN(weight)) {
1611  weight = 1.;
1612  }
1613  }
1614  else if (mParameterIndex == 44) { // Au+Au 11.5 GeV 2020
1615 
1616  if (iVzBinIndex < 0 || iVzBinIndex > auau11_run20_nVzBins) return 1.0;
1617 
1618  weight = auau11_run20_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1619  // Handle bad weight
1620  if (weight == 0 || TMath::IsNaN(weight)) {
1621  weight = 1.;
1622  }
1623  }
1624 
1625  else if (mRefX == 6 && mParameterIndex == 0) { // O+O 200 GeV 2021
1626 
1627  if (iVzBinIndex < 0 || iVzBinIndex > oo200_run21_nVzBins) return 1.0;
1628 
1629  weight = oo200_run21_shapeWeightArray_refMult[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1630  // Handle bad weight
1631  if (weight == 0 || TMath::IsNaN(weight)) {
1632  weight = 1.;
1633  }
1634  } // else if (mRefX == 6 && mParameterIndex == 0) { // O+O 200 GeV 2021
1635 
1636  else if (mRefX == 7 && mParameterIndex == 0) { // O+O 200 GeV 2021
1637 
1638  if (iVzBinIndex < 0 || iVzBinIndex > oo200_run21_nVzBins) return 1.0;
1639 
1640  weight = oo200_run21_shapeWeightArray_totnMIP[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1641  // Handle bad weight
1642  if (weight == 0 || TMath::IsNaN(weight)) {
1643  weight = 1.;
1644  }
1645  } // else if (mRefX == 7 && mParameterIndex == 0) { // O+O 200 GeV 2021
1646 
1647  else if (mRefX == 6 && mParameterIndex == 1) { // d+Au 200 GeV 2021
1648 
1649  if (iVzBinIndex < 0 || iVzBinIndex > dau200_run21_nVzBins) return 1.0;
1650 
1651  weight = dau200_run21_shapeWeightArray_refMult[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1652  // Handle bad weight
1653  if (weight == 0 || TMath::IsNaN(weight)) {
1654  weight = 1.;
1655  }
1656  } // else if (mRefX == 6 && mParameterIndex == 1) { // d+Au 200 GeV 2021
1657 
1658  else {
1659  weight = 1.0;
1660  }
1661 
1662  if (mVerbose) {
1663  std::cout << "\tShape weight: " << weight << std::endl;
1664  std::cout << "\t[DONE]\n";
1665  }
1666 
1667  return weight;
1668 }
1669 
1670 //________________
1671 Double_t StRefMultCorr::triggerWeight() const {
1672 
1673  if (mVerbose) {
1674  std::cout << "Estimating the trigger weight" << std::endl;
1675  }
1676 
1677  Double_t weight = 1.;
1678 
1679  const Double_t par0 = mPar_weight[0][mParameterIndex];
1680  const Double_t par1 = mPar_weight[1][mParameterIndex];
1681  const Double_t par2 = mPar_weight[2][mParameterIndex];
1682  const Double_t par3 = mPar_weight[3][mParameterIndex];
1683  const Double_t par4 = mPar_weight[4][mParameterIndex];
1684  const Double_t A = mPar_weight[5][mParameterIndex]; // Set to 0 if no z-dependent trigger weight is set
1685  const Double_t par6 = mPar_weight[6][mParameterIndex];
1686  const Double_t par7 = mPar_weight[7][mParameterIndex];
1687 
1688  if (mVerbose) {
1689  std::cout << "Trigger efficiency parameters: "
1690  << " [0]: " << par0
1691  << " [1]: " << par1
1692  << " [2]: " << par2
1693  << " [3]: " << par3
1694  << " [4]: " << par4
1695  << " A: " << A
1696  << " [6]: " << par6
1697  << " [7]: " << par7 << std::endl;
1698  }
1699 
1700  // Additional z-vetex dependent correction
1701  //const Double_t A = ((1.27/1.21))/(30.0*30.0); // Don't ask...
1702  //const Double_t A = (0.05/0.21)/(30.0*30.0); // Don't ask...
1703 
1704 
1705  if ( isRefMultOk() // 0-80%
1706  && mRefMult_corr < mNormalize_stop[mParameterIndex] // reweighting only apply up to normalization point
1707  && mRefMult_corr != -(par3/par2) ) { // avoid denominator = 0
1708 
1709  // Parametrization of MC/data RefMult ratio
1710  if (mRefX == 5 && mParameterIndex == 0) { // Run 18 Au+Au 3.85 GeV (sqrt(s_NN)=3 GeV)
1711  // Trigger efficiency correction does not exist. Temporarily set weight to 1
1712  weight = 1.;
1713  } // else if (mRefX == 5 && mParameterIndex == 0)
1714  else if ((mRefX == 6 || mRefX == 7) && (mParameterIndex == 0 || mParameterIndex == 1)) { // O+O and d+Au 200 GeV 2021
1715  weight = 1. + par0 * TMath::Exp(par1*mRefMult_corr + par2);
1716  } // else if ((mRefX == 6 || mRefX == 7) && (mParameterIndex == 0 || mParameterIndex == 1)) { // O+O and d+Au 200 GeV 2021
1717  else {
1718  weight = ( par0 +
1719  par1 / (par2 * mRefMult_corr + par3) +
1720  par4 * (par2 * mRefMult_corr + par3) +
1721  par6 / ( (par2 * mRefMult_corr + par3) * (par2 * mRefMult_corr + par3) ) +
1722  par7 * ( (par2 * mRefMult_corr + par3) * (par2 * mRefMult_corr + par3) ) );
1723  }
1724  /*
1725  std::cout << "par0: " << par0 << " par1: " << par1 << " par2: " << par2
1726  << " par3: " << par3 << " par4: " << par4 << " A: " << A
1727  << " par6: " << par6 << " par7: " << par7 << "\n"
1728  << "refMultCorr: " << mRefMult_corr << " weight: " << weight << std::endl;
1729  */
1730 
1731  weight = weight + (weight - 1.0) * (A * mVz * mVz); // z-dependent weight correction
1732  }
1733  if (mVerbose) {
1734  std::cout << "\tTrigger weight: " << weight << std::endl;
1735  std::cout << "\t[DONE]\n";
1736  }
1737  return weight;
1738 }
1739 
1740 //_________________
1741 Double_t StRefMultCorr::getWeight() const {
1742 
1743  Double_t Weight = 1.0;
1744 
1745  // Invalid index
1746  if( !isIndexOk() ) return Weight ;
1747  // Invalid z-vertex
1748  if( !isZvertexOk() ) return Weight ;
1749 
1750  // Trigger efficiency
1751  Weight = triggerWeight();
1752 
1753  //------------for Run14 and Run16----------------
1754  // Special scale factor for global refmult depending on Vz window
1755  // for others, scale factor = 1
1756  Weight *= getScaleForWeight();
1757 
1758  // Shape correction
1759  Weight *= getShapeWeight_SubVz2Center();
1760 
1761  return Weight;
1762 }
1763 
1764 //_________________
1765 Int_t StRefMultCorr::getCentralityBin16() const {
1766 
1767  if (mVerbose) {
1768  std::cout << "Finding centrality bin (out of 16)\n";
1769  }
1770  Int_t CentBin16 = -1;
1771 
1772  // Invalid index
1773  if( !isIndexOk() ) return CentBin16 ;
1774 
1775  while( CentBin16 < mNCentrality && !isCentralityOk(CentBin16) ) {
1776  CentBin16++;
1777  }
1778 
1779  // Vz dependent centrality definition for Vz<-27 and Vz>25
1780  // Run17 54.4 GeV, trigid = 580001
1781  if( mParameterIndex==28&&(mVz<-27||mVz>25) ) {
1782  CentBin16 = getCentralityBin16VzDep();
1783  // if(CentBin16==-1) std::cout <<"Vz="<< mVz <<" RefMult="<< mRefMult <<" RefMultCorr="<< mRefMult_corr << std::endl;
1784  // if(CentBin16==9999) std::cout << "Invalide number"<< std::endl;
1785  // std::cout <<"Vz dependent centrality definition for Vz<-27 and Vz>25 ... Vz="<< mVz <<" RefMult="<< mRefMult <<" RefMultCorr="<< mRefMult_corr <<" iCent="<< CentBin16 << std::endl;
1786  }
1787 
1788  if (mVerbose) {
1789  std::cout << "\tCentrality bin (out of 16): " << CentBin16 << "\t[DONE]" << std::endl;
1790  }
1791 
1792  // return -1 if CentBin16 = 16 (very large refmult, refmult>5000)
1793  return ( CentBin16==16 ) ? -1 : CentBin16;
1794 }
1795 
1796 //_________________
1797 Int_t StRefMultCorr::getCentralityBin9() const {
1798 
1799  if (mVerbose) {
1800  std::cout << "Finding centrality bin (out of 9)\n";
1801  }
1802 
1803  Int_t CentBin9 = -1;
1804 
1805  // Invalid index
1806  if ( !isIndexOk() ) return CentBin9 ;
1807 
1808  const Int_t CentBin16 = getCentralityBin16(); // Centrality bin 16
1809  const Bool_t isCentralityOk = CentBin16 >= 0 && CentBin16 < mNCentrality ;
1810 
1811  // No centrality is defined
1812  if ( !isCentralityOk ) return CentBin9;
1813 
1814  // First handle the exceptions
1815  if( mRefMult_corr > mCentrality_bins[15][mParameterIndex] &&
1816  mRefMult_corr <= mCentrality_bins[16][mParameterIndex] ) {
1817  CentBin9 = 8; // most central 5%
1818  }
1819  else if( mRefMult_corr > mCentrality_bins[14][mParameterIndex] &&
1820  mRefMult_corr <= mCentrality_bins[15][mParameterIndex]) {
1821  CentBin9 = 7; // most central 5-10%
1822  }
1823  else {
1824  CentBin9 = (Int_t)(0.5*CentBin16);
1825  }
1826 
1827  // Vz dependent centrality definition for Vz<-27 and Vz>25
1828  // Run17 54.4 GeV, trigid = 580001
1829  // std::cout << mParameterIndex << std::endl;
1830  if ( mParameterIndex == 28 && ( mVz<-27 || mVz>25 ) ) {
1831  CentBin9 = getCentralityBin9VzDep();
1832  }
1833 
1834  if (mVerbose) {
1835  std::cout << "\tCentrality bin (out of 9): " << CentBin9 << "\t[DONE]" << std::endl;
1836  }
1837 
1838  return CentBin9;
1839 }
1840 
1841 //_________________
1842 Int_t StRefMultCorr::getVzWindowForVzDepCentDef() const {
1843  Int_t iBinVz = -1;
1844 
1845  if( mParameterIndex==28 ) { // 54.4 GeV, RefMult, 580001
1846  if( mVz>-30 && mVz<-29 ) iBinVz = 0;
1847  else if( mVz>-29 && mVz<-27 ) iBinVz = 1;
1848  else if( mVz>25 && mVz<27 ) iBinVz = 2;
1849  else if( mVz>27 && mVz<29 ) iBinVz = 3;
1850  else if( mVz>29 && mVz<30 ) iBinVz = 4;
1851  else iBinVz = -1;
1852  }
1853  else if( mParameterIndex>=30 && mParameterIndex<=35 ) { //Run18 27 GeV MB, 6 triggerIds
1854  if( mVz>=-70. && mVz<-60. ) iBinVz = 0;
1855  else if( mVz>=-60.&&mVz<-50.) iBinVz = 1;
1856  else if( mVz>=-50.&&mVz<-40.) iBinVz = 2;
1857  else if( mVz>=-40.&&mVz<-30.) iBinVz = 3;
1858  else if( mVz>=-30.&&mVz<-20.) iBinVz = 4;
1859  else if( mVz>=-20.&&mVz<-10.) iBinVz = 5;
1860  else if( mVz>=-10.&&mVz<0.0 ) iBinVz = 6;
1861  else if( mVz>=0.0 &&mVz<10. ) iBinVz = 7;
1862  else if( mVz>=10. &&mVz<20. ) iBinVz = 8;
1863  else if( mVz>=20. &&mVz<30. ) iBinVz = 9;
1864  else if( mVz>=30. &&mVz<40. ) iBinVz = 10;
1865  else if( mVz>=40. &&mVz<50. ) iBinVz = 11;
1866  else if( mVz>=50. &&mVz<60. ) iBinVz = 12;
1867  else if( mVz>=60. &&mVz<=70 ) iBinVz = 13;
1868  else iBinVz = -1;
1869  }
1870  else if( mParameterIndex>=36 && mParameterIndex<=37 ) { //Run18 200 GeV isobar
1871  Double_t VtxZBinDouble = mVz/2. + 17.;
1872  iBinVz = 0;
1873  if (mVz == 25.) {
1874  iBinVz = 29;
1875  }
1876  else if (mVz == -35.) {
1877  iBinVz = 0;
1878  }
1879  else {
1880  iBinVz = TMath::Nint(VtxZBinDouble);
1881  }
1882  }
1883  else if ( (mParameterIndex == 38) || (mParameterIndex == 39) ) { // Au+Au 19.6 GeV 2019 || Au+Au 14.6 GeV 2019
1884 
1885  for ( Int_t iVz=0; iVz<auau19_run19_nVzBins; iVz++ ) { // Utilize same Vz binning: (29 bins, -145., 145.)
1886  if ( auau19_run19_vzRangeLimits[iVz][0] <= mVz && mVz < auau19_run19_vzRangeLimits[iVz][1] ) {
1887  iBinVz = iVz;
1888  break;
1889  }
1890  } // for ( Int_t iVz=0; iVz<auau19_run19_nVzBins; iVz++ )
1891  } // else if ( mParameterIndex == 38 )
1892  else if ( mParameterIndex == 40 ) { // Au+Au 200 GeV 2019
1893  for ( Int_t iVz=0; iVz<auau200_run19_nVzBins; iVz++ ) { // Utilize Vz binning: (20 bins, -100., 100.)
1894  if ( auau200_run19_vzRangeLimits[iVz][0] <= mVz && mVz < auau200_run19_vzRangeLimits[iVz][1] ) {
1895  iBinVz = iVz;
1896  break;
1897  }
1898  } // for ( Int_t iVz=0; iVz<auau200_run19_nVzBins; iVz++ )
1899  } // else if ( mParameterIndex == 40 )
1900  else if ( mParameterIndex == 41 ) { // Au+Au 7.7 GeV 2020
1901  for ( Int_t iVz=0; iVz<auau7_run21_nVzBins; iVz++ ) { // Utilize Vz binning: (29 bins, -145., 145.)
1902  if ( auau7_run21_vzRangeLimits[iVz][0] <= mVz && mVz < auau7_run21_vzRangeLimits[iVz][1] ) {
1903  iBinVz = iVz;
1904  break;
1905  }
1906  } // for ( Int_t iVz=0; iVz<auau7_run21_nVzBins; iVz++ )
1907  } // else if ( mParameterIndex == 41 )
1908  else if ( mParameterIndex == 42 ) { // Au+Au 9.2 GeV Run 20 TriggerID = 780020
1909  for ( Int_t iVz=0; iVz<auau9_run20_nVzBins; iVz++ ) { // Utilize Vz binning: (29 bins, -145., 145.)
1910  if ( auau9_run20_vzRangeLimits[iVz][0] <= mVz && mVz < auau9_run20_vzRangeLimits[iVz][1] ) {
1911  iBinVz = iVz;
1912  break;
1913  }
1914  } // for ( Int_t iVz=0; iVz<auau9_run20_nVzBins; iVz++ )
1915  } // else if ( mParameterIndex == 42 )
1916  else if ( mParameterIndex == 43 ) { // Au+Au 17.3 GeV 2021
1917  for ( Int_t iVz=0; iVz<auau17_run21_nVzBins; iVz++ ) { // Utilize Vz binning: (29 bins, -145., 145.)
1918  if ( auau17_run21_vzRangeLimits[iVz][0] <= mVz && mVz < auau17_run21_vzRangeLimits[iVz][1] ) {
1919  iBinVz = iVz;
1920  break;
1921  }
1922  } // for ( Int_t iVz=0; iVz<auau17_run21_nVzBins; iVz++ )
1923  } // else if ( mParameterIndex == 43 )
1924  else if ( mParameterIndex == 44 ) { // Au+Au 11.5 GeV 2020
1925  for ( Int_t iVz=0; iVz<auau11_run20_nVzBins; iVz++ ) { // Utilize Vz binning: (29 bins, -145., 145.)
1926  if ( auau11_run20_vzRangeLimits[iVz][0] <= mVz && mVz < auau11_run20_vzRangeLimits[iVz][1] ) {
1927  iBinVz = iVz;
1928  break;
1929  }
1930  } // for ( Int_t iVz=0; iVz<auau17_run21_nVzBins; iVz++ )
1931  } // else if ( mParameterIndex == 44 )
1932 
1933  else if ( (mRefX == 6 || mRefX == 7) && mParameterIndex == 0 ) { // O+O 200 GeV 2021
1934  for ( Int_t iVz=0; iVz<oo200_run21_nVzBins; iVz++ ) {
1935  if ( oo200_run21_vzRangeLimits[iVz][0] <= mVz && mVz < oo200_run21_vzRangeLimits[iVz][1] ) {
1936  iBinVz = iVz;
1937  break;
1938  }
1939  } // for ( Int_t iVz=0; iVz<oo200_run21_nVzBins; iVz++ )
1940  } // else if ( (mRefX == 6 || mRefX == 7) && mParameterIndex == 0 ) { // O+O 200 GeV 2021
1941 
1942  else if ( mRefX == 6 && mParameterIndex == 1 ) { // d+Au 200 GeV 2021
1943  for ( Int_t iVz=0; iVz<dau200_run21_nVzBins; iVz++ ) {
1944  if ( dau200_run21_vzRangeLimits[iVz][0] <= mVz && mVz < dau200_run21_vzRangeLimits[iVz][1] ) {
1945  iBinVz = iVz;
1946  break;
1947  }
1948  } // for ( Int_t iVz=0; iVz<dau200_run21_nVzBins; iVz++ )
1949  } // else if ( mRefX == 6 && mParameterIndex == 1 ) { // d+Au 200 GeV 2021
1950 
1951  else {
1952  iBinVz = -1;
1953  }
1954 
1955  return iBinVz;
1956 }
1957 
1958 //_________________
1959 Int_t StRefMultCorr::getCentralityBin9VzDep() const {
1960  const Int_t vzid = getVzWindowForVzDepCentDef();
1961  Int_t iCent = 9999;
1962  for(Int_t i=0; i<9; i++) {
1963  if ( i==8 ) {
1964  if( mRefMult_corr>CentBin9_vzdep[vzid][i] && mRefMult_corr<50000 ) iCent = i;
1965  }
1966  else if( mRefMult_corr>CentBin9_vzdep[vzid][i] &&
1967  mRefMult_corr<CentBin9_vzdep[vzid][i+1] ) iCent = i;
1968  }
1969  // 80-100% for icent=-1
1970  if( mRefMult_corr>0 && mRefMult_corr<CentBin9_vzdep[vzid][0] ) iCent = -1;
1971  return ( iCent==9999 ) ? -1 : iCent;
1972 }
1973 
1974 //_________________
1975 Int_t StRefMultCorr::getCentralityBin16VzDep() const {
1976  const Int_t vzid = getVzWindowForVzDepCentDef();
1977  Int_t iCent = 9999;
1978  for(Int_t i=0; i<16; i++) {
1979  if (i == 15) {
1980  if (mRefMult_corr > CentBin16_vzdep[vzid][i] && mRefMult_corr < 50000) {
1981  iCent = i;
1982  }
1983  }
1984  else if (mRefMult_corr > CentBin16_vzdep[vzid][i] &&
1985  mRefMult_corr < CentBin16_vzdep[vzid][i + 1]) {
1986  iCent = i;
1987  }
1988  }
1989  // 80-100% for icent=-1
1990  if( mRefMult_corr>0 && mRefMult_corr<CentBin16_vzdep[vzid][0] ) {
1991  iCent = -1;
1992  }
1993  return (iCent==9999) ? -1 : iCent;
1994 }
1995 
1996 //_________________
1997 const Int_t StRefMultCorr::getRefX() const {
1998  if ( mName.CompareTo("grefmult", TString::kIgnoreCase) == 0 ) return 0;
1999  else if ( mName.CompareTo("refmult", TString::kIgnoreCase) == 0 ) return 1;
2000  else if ( mName.CompareTo("refmult2", TString::kIgnoreCase) == 0 ) return 2;
2001  else if ( mName.CompareTo("refmult3", TString::kIgnoreCase) == 0 ) return 3;
2002  else if ( mName.CompareTo("refmult4", TString::kIgnoreCase) == 0 ) return 4;
2003  else if ( mName.CompareTo("fxtmult", TString::kIgnoreCase) == 0 ) return 5;
2004  else if ( mName.CompareTo("refmult6", TString::kIgnoreCase) == 0 ) return 6;
2005  else if ( mName.CompareTo("totnmip", TString::kIgnoreCase) == 0 ) return 7;
2006  else return 9999;
2007 }
2008 
2009 //_________________
2010 const Int_t StRefMultCorr::getNumberOfDatasets() const {
2011  if ( mName.CompareTo("grefmult", TString::kIgnoreCase) == 0 ) return nID_gref;
2012  else if ( mName.CompareTo("refmult", TString::kIgnoreCase) == 0 ) return nID_ref1;
2013  else if ( mName.CompareTo("refmult2", TString::kIgnoreCase) == 0 ) return nID_ref2;
2014  else if ( mName.CompareTo("refmult3", TString::kIgnoreCase) == 0 ) return nID_ref3;
2015  else if ( mName.CompareTo("refmult4", TString::kIgnoreCase) == 0 ) return nID_ref4;
2016  else if ( mName.CompareTo("fxtmult", TString::kIgnoreCase) == 0 ) return nID_ref5;
2017  else if ( mName.CompareTo("refmult6", TString::kIgnoreCase) == 0 ) return nID_ref6;
2018  else if ( mName.CompareTo("totnmip", TString::kIgnoreCase) == 0 ) return nID_ref7;
2019  else return 9999;
2020 }
2021 
2022 //_________________
2023 void StRefMultCorr::readHeaderFile() {
2024 
2025  //std::vector<std::string> sParam_ShapeWeight = StringSplit(getParamX_ShapeWeight(1,1),',');
2026  //for(int ib=0;ib<sParam_ShapeWeight.size(); ib++) std::cout<<"sParam_ShapeWeight[i]: "<<sParam_ShapeWeight[ib]<<std::endl;
2027 
2028  // Read reference multiplicity to work with (grefmult, refMult, refMult2, refMult3, fxtMult)
2029  const Int_t refX = getRefX();
2030  mRefX = refX;
2031  // Retrieve number of datasets for the given mRefX
2032  const Int_t nID = getNumberOfDatasets();
2033 
2034  for(int iID=0; iID<nID; iID++) {
2035 
2036  //
2037  // Year, energy, run numbers, Vz cut
2038  //
2039  Int_t year; Double_t energy;
2040  std::vector<std::string> sParam = StringSplit(getParamX(refX,iID,0),':');
2041  year = stoi(sParam[0]);
2042  energy = stoi(sParam[1]);
2043  std::vector<std::string> sRuns = StringSplit(sParam[2],',');
2044 
2045  Int_t startRunId=0, stopRunId=0;
2046  startRunId = stoi(sRuns[0]);
2047  stopRunId = stoi(sRuns[1]);
2048 
2049  Double_t startZvertex=-9999., stopZvertex=-9999. ;
2050  std::vector<std::string> sVz = StringSplit(sParam[3],',');
2051  startZvertex = stod(sVz[0]);
2052  stopZvertex = stod(sVz[1]);
2053 
2054  mYear.push_back(year);
2055  mBeginRun.insert(std::make_pair(std::make_pair(energy, year), startRunId));
2056  mEndRun.insert( std::make_pair(std::make_pair(energy, year), stopRunId));
2057  mStart_runId.push_back( startRunId ) ;
2058  mStop_runId.push_back( stopRunId ) ;
2059  mStart_zvertex.push_back( startZvertex ) ;
2060  mStop_zvertex.push_back( stopZvertex ) ;
2061 
2062  //
2063  // Centrality definition
2064  //
2065  std::vector<std::string> sParamCent = StringSplit(getParamX(refX,iID,1),',');
2066  for(UInt_t i=0; i<sParamCent.size(); i++) {
2067  mCentrality_bins[i].push_back( stoi(sParamCent[i]) );
2068  }
2069  mCentrality_bins[mNCentrality].push_back( 5000 );
2070 
2071  //
2072  // Normalization range
2073  //
2074  Double_t normalize_stop=-1.0 ;
2075  normalize_stop = stod(getParamX(refX,iID,2)) ;
2076  mNormalize_stop.push_back( normalize_stop );
2077 
2078  //
2079  // Acceptance (vz) correction
2080  //
2081  std::vector<std::string> sParamVz = StringSplit(getParamX(refX,iID,3),',');
2082  for(UInt_t i=0; i<mNPar_z_vertex; i++) {
2083  Double_t val = -9999.;
2084  if(i<sParamVz.size()) val = stod(sParamVz[i]);
2085  else val = 0.0;
2086  mPar_z_vertex[i].push_back( val );
2087  }
2088 
2089  //
2090  // Trigger inefficiency correction
2091  //
2092  std::vector<std::string> sParamTrig = StringSplit(getParamX(refX,iID,4),',');
2093  for(UInt_t i=0; i<mNPar_weight; i++) {
2094  Double_t val = -9999.;
2095  if(i<sParamTrig.size()) val = stod(sParamTrig[i]);
2096  else val = 0.0;
2097  mPar_weight[i].push_back( val );
2098  }
2099 
2100  //
2101  // Luminosity correction
2102  //
2103  std::vector<std::string> sParamLumi = StringSplit(getParamX(refX,iID,5),',');
2104  for(UInt_t i=0; i<mNPar_luminosity; i++) {
2105  Double_t val = -9999.;
2106  if(i<sParamLumi.size()) val = stod(sParamLumi[i]);
2107  else val = 0.0;
2108  mPar_luminosity[i].push_back( val );
2109  }
2110 
2111  // std::cout << refX <<" "<< iID <<"/"<< nID << std::endl;
2112  }
2113 
2114  std::cout << "StRefMultCorr::readHeaderFile [" << mName
2115  << "] Correction parameters and centrality definitions have been successfully read for refX=" << refX
2116  << std::endl;
2117 }
2118 
2119 //_________________
2120 void StRefMultCorr::readBadRunsFromHeaderFile() {
2121 
2122  //
2123  // TODO: Modify to read only bad runs associated with the requested year
2124  //
2125 
2126  for(Int_t i=0; i<nBadRun_refmult_2010; i++) {
2127  mBadRun.push_back(badrun_refmult_2010[i]);
2128  }
2129  std::cout << "read in nBadRun_refmult_2010: " << nBadRun_refmult_2010 << std::endl;
2130 
2131  for(Int_t i=0; i<nBadRun_refmult_2011; i++) {
2132  mBadRun.push_back(badrun_refmult_2011[i]);
2133  }
2134  std::cout << "read in nBadRun_refmult_2011: " << nBadRun_refmult_2011 << std::endl;
2135 
2136  for(Int_t i=0; i<nBadRun_grefmult_2014; i++) {
2137  mBadRun.push_back(badrun_grefmult_2014[i]);
2138  }
2139  std::cout << "read in nBadRun_grefmult_2014: " << nBadRun_grefmult_2014 << std::endl;
2140 
2141  for(Int_t i=0; i<nBadRun_grefmult_2016; i++) {
2142  mBadRun.push_back(badrun_grefmult_2016[i]);
2143  }
2144  std::cout << "read in nBadRun_grefmult_2016: " << nBadRun_grefmult_2016 << std::endl;
2145 
2146  for(Int_t i=0; i<nBadRun_refmult_2017; i++) {
2147  mBadRun.push_back(badrun_refmult_2017[i]);
2148  }
2149  std::cout << "read in nBadRun_refmult_2017: " << nBadRun_refmult_2017 << std::endl;
2150 
2151  for (Int_t i = 0; i < nBadRun_refmult_2018; i++) {
2152  mBadRun.push_back(badrun_refmult_2018[i]);
2153  }
2154  std::cout << "read in nBadRun_refmult_2018: " << nBadRun_refmult_2018 << std::endl;
2155 
2156  for (Int_t i = 0; i < nBadRun_refmult_2019; i++) {
2157  mBadRun.push_back(badrun_refmult_2019[i]);
2158  }
2159  std::cout << "read in nBadRun_refmult_2019: " << nBadRun_refmult_2019 << std::endl;
2160 
2161  for (Int_t i = 0; i < nBadRun_refmult_2020; i++) {
2162  mBadRun.push_back(badrun_refmult_2020[i]);
2163  }
2164  std::cout << "read in nBadRun_refmult_2020: " << nBadRun_refmult_2020 << std::endl;
2165 
2166  for (Int_t i = 0; i < nBadRun_refmult_2021; i++) {
2167  mBadRun.push_back(badrun_refmult_2021[i]);
2168  }
2169  std::cout << "read in nBadRun_refmult_2021: " << nBadRun_refmult_2021 << std::endl;
2170 
2172  if ( mName.CompareTo("grefmult", TString::kIgnoreCase) == 0 ) {
2173  std::cout << "StRefMultCorr::readBadRunsFromHeaderFile Bad runs for year 2010, 2011, 2017, 2018 and 2019 have been read." << std::endl;
2174  }
2175 }
2176 
2177 //_________________
2178 void StRefMultCorr::print(const Option_t* option) const {
2179  std::cout << "StRefMultCorr::print Print input parameters for "
2180  << mName << " ========================================" << std::endl << std::endl;
2181  // Option switched off, can be used to specify parameters
2182  // const TString opt(option);
2183 
2184  // Int_t input_counter = 0;
2185  for(UInt_t id=0; id<mStart_runId.size(); id++) {
2186  //std::cout << "Data line = " << input_counter << ", Start_runId = " << Start_runId[input_counter] << ", Stop_runId = " << Stop_runId[input_counter] << std::endl;
2187  // const UInt_t id = mStart_runId.size()-1;
2188 
2189  // Removed line break
2190  std::cout << " Index=" << id;
2191  std::cout << Form(" Run=[%8d, %8d]", mStart_runId[id], mStop_runId[id]);
2192  std::cout << Form(" z-vertex=[%1.1f, %1.1f]", mStart_zvertex[id], mStop_zvertex[id]);
2193  std::cout << ", Normalize_stop=" << mNormalize_stop[id];
2194  std::cout << std::endl;
2195 
2196  // if(opt.IsWhitespace()){
2197  // continue ;
2198  // }
2199 
2200  std::cout << "Centrality: ";
2201 
2202  for(Int_t i=0;i<mNCentrality;i++) {
2203  std::cout << Form(" >%2d%%", 80-5*i);
2204  }
2205  std::cout << std::endl;
2206  std::cout << "RefMult: ";
2207 
2208  for(Int_t i=0;i<mNCentrality;i++) {
2209  // std::cout << Form("StRefMultCorr::read Centrality %3d-%3d %%, refmult > %4d", 75-5*i, 80-5*i, mCentrality_bins[i][id]) << std::endl;
2210  const TString tmp(">");
2211  const TString centrality = tmp + Form("%d", mCentrality_bins[i][id]);
2212  std::cout << Form("%6s", centrality.Data());
2213  }
2214  std::cout << std::endl;
2215 
2216  for(Int_t i=0;i<mNPar_z_vertex;i++) {
2217  std::cout << " mPar_z_vertex[" << i << "] = " << mPar_z_vertex[i][id];
2218  }
2219  std::cout << std::endl;
2220 
2221  for(Int_t i=0;i<mNPar_weight;i++) {
2222  std::cout << " mPar_weight[" << i << "] = " << mPar_weight[i][id];
2223  }
2224  std::cout << std::endl;
2225 
2226  for(Int_t i=0;i<mNPar_luminosity;i++) {
2227  std::cout << " mPar_luminosity[" << i << "] = " << mPar_luminosity[i][id];
2228  }
2229  std::cout << std::endl << std::endl;
2230  }
2231  std::cout << "=====================================================================================" << std::endl;
2232 }
2233 
2234 //_________________
2235 std::vector<std::string> StRefMultCorr::StringSplit( const std::string str, const char sep ) const {
2236  std::vector<std::string> vstr;
2237  std::stringstream ss(str);
2238  std::string buffer;
2239  while( getline(ss,buffer,sep) ) {
2240  vstr.push_back(buffer);
2241  }
2242  return vstr;
2243 }
2244 
2245 //________________
2246 Double_t StRefMultCorr::calcPileUpRefMult(Double_t ntofmatch, Double_t x0, Double_t x1,
2247  Double_t x2, Double_t x3, Double_t x4) const {
2248  return ( x0 + x1*(ntofmatch) + x2*pow(ntofmatch,2) + x3*pow(ntofmatch,3) + x4*pow(ntofmatch,4) );
2249 }
2250 
2251 //________________
2252 Bool_t StRefMultCorr::isPileUpEvent(Double_t refmult, Double_t ntofmatch, Double_t vz, Double_t totnMIP) const {
2253  if ((mRefX==6) || (mRefX==7)) {
2254  // refMult6 and totnMIP require both refMult vs. nBTOFMatch and totnMIP vs. nBTOFMatch for pileup rejection
2255  if (totnMIP < 0.) {
2256  Error("StRefMultCorr::isPileUpEvent", "totnMIP<0");
2257  return kTRUE;
2258  } // if (totnMIP < 0.)
2259  return !( passnTofMatchRefmultCut(refmult, ntofmatch, vz) && passnTofMatchTotnMIPCut(totnMIP, ntofmatch, vz) );
2260  } // if ((mRefX==6) || (mRefX==7))
2261  else { // other refMult
2262  return !passnTofMatchRefmultCut(refmult, ntofmatch, vz);
2263  }
2264 }
StRefMultCorr::luminosityCorrection
Double_t luminosityCorrection(Double_t zdcCoincidenceRate) const
Luminosity correction factor.
Definition: StRefMultCorr.cxx:1146
Error
Definition: FJcore.h:367
StRefMultCorr::getCentralityBin9
Int_t getCentralityBin9() const
Get 9 centrality bins (10% increment except for 0-5 and 5-10)
Definition: StRefMultCorr.cxx:1797
StRefMultCorr::getShapeWeight_SubVz2Center
Double_t getShapeWeight_SubVz2Center() const
Shape reweighting of refmult: ratio of refMult in each Vz bin to that in the center (|Vz|&lt;10cm) ...
Definition: StRefMultCorr.cxx:1482
StRefMultCorr::getCentralityBin16
Int_t getCentralityBin16() const
Get 16 centrality bins (5% increment, 0-5, 5-10, ..., 75-80)
Definition: StRefMultCorr.cxx:1765
StRefMultCorr::getBeginRun
Int_t getBeginRun(const Double_t energy, const Int_t year)
Return the first runId from energy and year.
Definition: StRefMultCorr.cxx:53
StRefMultCorr::sampleRefMult
Double_t sampleRefMult(Int_t refMult) const
Sample refMult -&gt; convert integer to double.
Definition: StRefMultCorr.cxx:1288
StRefMultCorr::getRefMultCorr
Double_t getRefMultCorr() const
Get corrected multiplicity, correction as a function of primary z-vertex.
Definition: StRefMultCorr.cxx:1140
StRefMultCorr::print
void print(const Option_t *option="") const
Print all parameters.
Definition: StRefMultCorr.cxx:2178
StRefMultCorr::getWeight
Double_t getWeight() const
Total weighting factor: incorporates shape and trigger efficiency weights.
Definition: StRefMultCorr.cxx:1741
StRefMultCorr::triggerWeight
Double_t triggerWeight() const
Trigger efficiency: fit of the Glauber/Data.
Definition: StRefMultCorr.cxx:1671
StRefMultCorr::passnTofMatchRefmultCut
Bool_t passnTofMatchRefmultCut(Double_t refmult, Double_t ntofmatch, Double_t vz=0.) const
Check if NOT pile-up event.
Definition: StRefMultCorr.cxx:160
StRefMultCorr::vzCorrection
Double_t vzCorrection(Double_t z) const
Vz correction factor.
Definition: StRefMultCorr.cxx:1201
StRefMultCorr::isBadRun
Bool_t isBadRun(const Int_t RunId)
Check if run is bad.
Definition: StRefMultCorr.cxx:125
StRefMultCorr::getEndRun
Int_t getEndRun(const Double_t energy, const Int_t year)
Return the last runId from energy and year.
Definition: StRefMultCorr.cxx:69
StRefMultCorr::StRefMultCorr
StRefMultCorr(const TString name="refmult", const TString subname="Def", const TString libname="Def")
Definition: StRefMultCorr.cxx:27
StRefMultCorr::~StRefMultCorr
virtual ~StRefMultCorr()
Destructor.
Definition: StRefMultCorr.cxx:48
StRefMultCorr::isPileUpEvent
Bool_t isPileUpEvent(Double_t refmult, Double_t ntofmatch, Double_t vz=0., Double_t totnMIP=-999.) const
Check if pile-up event.
Definition: StRefMultCorr.cxx:2252
StRefMultCorr::init
void init(const Int_t RunId)
Initialization of centrality bins etc.
Definition: StRefMultCorr.cxx:1005
StRefMultCorr::passnTofMatchTotnMIPCut
Bool_t passnTofMatchTotnMIPCut(Double_t totnMIP, Double_t ntofmatch, Double_t vz=0.) const
Check if NOT pile-up event using EPD&#39;s total nMIP versus nBTOFMatch.
Definition: StRefMultCorr.cxx:853
Generated by   doxygen 1.8.5