11 #include "StRefMultCorr.h"
22 typedef std::pair<Double_t, Int_t> keys;
28 const TString libname) :
29 mName(name), mSubName(subname), mLibName(libname),
30 mRefX(1), mVerbose(kFALSE) {
37 std::cout << mSubName.Data() <<
" "<< mLibName.Data() << std::endl;
43 readBadRunsFromHeaderFile() ;
54 keys key(std::make_pair(energy, year));
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);
63 std::pair<std::multimap<keys, Int_t>::iterator, std::multimap<keys, Int_t>::iterator> iterRange = mBeginRun.equal_range(key);
65 return (*(iterRange.first)).second ;
70 keys key(std::make_pair(energy, year));
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);
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;
83 return (*iter).second;
87 void StRefMultCorr::clear() {
91 mStart_runId.clear() ;
93 mStart_zvertex.clear() ;
94 mStop_zvertex.clear() ;
95 mNormalize_stop.clear() ;
97 for(Int_t i=0; i<mNCentrality; i++) {
98 mCentrality_bins[i].clear() ;
101 mParameterIndex = -1 ;
103 for(Int_t i=0;i<mNPar_z_vertex;i++) {
104 mPar_z_vertex[i].clear() ;
107 for(Int_t i=0;i<mNPar_weight;i++) {
108 mPar_weight[i].clear();
111 for(Int_t i=0;i<mNPar_luminosity;i++) {
112 mPar_luminosity[i].clear();
119 mnVzBinForWeight = 0 ;
120 mVzEdgeForWeight.clear();
121 mgRefMultTriggerCorrDiffVzScaleRatio.clear() ;
127 std::vector<Int_t>::iterator iter = std::find(mBadRun.begin(), mBadRun.end(), RunId);
129 if ( iter != mBadRun.end() ) {
131 std::cout <<
"StRefMultCorr::isBadRun Find bad run = " << (*iter) << std::endl;
135 return ( iter != mBadRun.end() ) ;
139 void StRefMultCorr::initEvent(
const UShort_t RefMult,
const Double_t z,
140 const Double_t zdcCoincidenceRate) {
142 std::cout <<
"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
143 std::cout <<
"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
144 std::cout <<
"+++++++++++++++++++++ New Event +++++++++++++++++++\n";
145 std::cout <<
"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
146 std::cout <<
"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
151 if ( mRefMult != RefMult || mVz != z || mZdcCoincidenceRate != zdcCoincidenceRate ) {
154 mZdcCoincidenceRate = zdcCoincidenceRate ;
155 mRefMult_corr =
getRefMultCorr(mRefMult, mVz, mZdcCoincidenceRate) ;
161 Double_t vz )
const {
165 std::cout <<
"Pile up check...";
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{};
173 Bool_t notPileUp = kFALSE;
181 if( mParameterIndex>=30 && mParameterIndex<=35 ) {
182 const Double_t min = 4.0;
183 const Double_t max = 5.0;
185 if(ntofmatch<=2)
return false;
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;
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);
207 refmultcutmin = refmultCenter - min*refmultLower;
208 refmultcutmax = refmultCenter + max*refmultUpper;
209 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
211 else if( mParameterIndex>=36 && mParameterIndex<=37 ) {
214 if(mParameterIndex==36) {
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;
227 else if(mParameterIndex==37) {
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;
239 refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
240 refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
241 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
243 else if ( mParameterIndex==38 ) {
245 if ( -145. <= vz && vz < -87. ) {
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;
257 else if ( -87. <= vz && vz < -29. ) {
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;
269 else if ( -29. <= vz && vz < 29. ) {
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;
281 else if ( 29. <= vz && vz < 87. ) {
284 b2=-0.00315093393202473;
285 b3=8.35823870487966e-06;
286 b4=-9.14822467807924e-09;
287 c0=-15.9411138444366;
289 c2=0.0011824174541949;
290 c3=-1.48902496972716e-06;
291 c4=-2.29371463231934e-10;
293 else if ( 87. <= vz && vz <= 145. ) {
296 b2=-0.00319461487211403;
297 b3=9.56612691680003e-06;
298 b4=-1.31049413530369e-08;;
299 c0=-22.4679773305418;
301 c2=-0.000106213494253586;
302 c3=4.93946486222714e-06;
303 c4=-1.1450298089717e-08;
306 refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
307 refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
308 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
310 else if ( mParameterIndex==39 ) {
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);
325 else if ( mParameterIndex==40 ) {
338 refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
339 refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
340 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
342 else if ( mParameterIndex==41 ) {
344 if ( -145. <= vz && vz < -87. ) {
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;
356 else if ( -87. <= vz && vz < -29. ) {
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;
368 else if ( -29. <= vz && vz < 29. ) {
371 b2=-0.00100184372825271;
372 b3=7.76378744751984e-07;
373 b4=-6.46469867000365e-09;
374 c0=-11.4340781454132;
376 c2=0.00121092416745035;
377 c3=1.17875404059176e-07;
378 c4=-9.81658682040738e-09;
380 else if ( 29. <= vz && vz < 87. ) {
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;
392 else if ( 87. <= vz && vz <= 145. ) {
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;
405 refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
406 refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
407 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
409 else if ( mParameterIndex==42 ) {
411 if ( -145. <= vz && vz < -87. ) {
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;
423 else if ( -87. <= vz && vz < -29. ) {
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;
435 else if ( -29. <= vz && vz < 29. ) {
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;
447 else if ( 29. <= vz && vz < 87. ) {
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;
459 else if ( 87. <= vz && vz <= 145. ) {
462 b2=-0.00195374948885512;
463 b3=-6.14244087431038e-06;
464 b4=1.99930095058841e-08;
465 c0=-15.6624325989392;
467 c2=0.00794996911844969;
468 c3=-4.09239155250494e-05;
469 c4=6.40163739983216e-08;
472 refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
473 refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
474 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
476 else if ( mParameterIndex==43 ) {
478 if ( -145. <= vz && vz < -87. ) {
481 b2=-0.00230107205687879;
482 b3=1.04069753338853e-05;
483 b4=-2.43265995270951e-08;
486 c2=-0.00285234327923795;
487 c3=1.68279361312683e-05;
488 c4=-2.89872992178789e-08;
490 else if ( -87. <= vz && vz < -29. ) {
493 b2=-0.000197781802002694;
494 b3=1.02666189094347e-06;
495 b4=-5.52762010064236e-09;
496 c0=-21.4352021999217;
498 c2=-0.00160328567162831;
499 c3=8.94486444751978e-06;
500 c4=-1.46093145276812e-08;
502 else if ( -29. <= vz && vz < 29. ) {
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;
514 else if ( 29. <= vz && vz < 87. ) {
517 b2=-5.45184310087876e-05;
518 b3=6.25643605701836e-07;
519 b4=-4.90542835006027e-09;
520 c0=-20.7158089395719;
522 c2=-0.00138806953636318;
523 c3=7.92595642206008e-06;
524 c4=-1.32107375325913e-08;
526 else if ( 87. <= vz && vz <= 145. ) {
529 b2=-0.000569887807630565;
530 b3=3.95821109316978e-06;
531 b4=-1.60367555403757e-08;
532 c0=-26.3129222166004;
534 c2=-0.00341644731702994;
535 c3=1.84782571448044e-05;
536 c4=-3.03333077890128e-08;
539 refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
540 refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
541 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
543 else if ( mParameterIndex==44 ) {
545 if ( -145. <= vz && vz < -87. ) {
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;
557 else if ( -87. <= vz && vz < -29. ) {
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;
569 else if ( -29. <= vz && vz < 29. ) {
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;
581 else if ( 29. <= vz && vz < 87. ) {
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;
593 else if ( 87. <= vz && vz <= 145. ) {
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;
606 refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
607 refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
608 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
615 else if ( mRefX == 5 ) {
616 if (mParameterIndex == 1) {
628 refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
629 refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
630 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
632 else if (mParameterIndex == 2) {
644 refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
645 refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
646 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
648 else if (mParameterIndex == 3) {
660 refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
661 refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
662 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
664 else if (mParameterIndex == 4) {
676 refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
677 refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
678 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
680 else if (mParameterIndex == 5) {
683 b2=-0.000336278299464254;
684 b3=-0.000549204114892259;
685 b4=1.89274000668251e-06;
686 c0=-24.3335976220474;
688 c2=0.0042539477677528;
689 c3=0.000725893349545147;
690 c4=-6.29726910263091e-06;
692 refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
693 refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
694 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
696 else if (mParameterIndex == 6) {
708 refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
709 refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
710 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
712 else if (mParameterIndex == 7) {
724 refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
725 refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
726 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
728 else if (mParameterIndex == 8) {
740 refmultcutmax = calcPileUpRefMult(ntofmatch, b0, b1, b2, b3, b4);
741 refmultcutmin = calcPileUpRefMult(ntofmatch, c0, c1, c2, c3, c4);
742 notPileUp = isInPileUpRefMultLimits(refmult, refmultcutmin, refmultcutmax);
773 std::cout <<
"\t notPileUp: ";
778 std::cout <<
"FALSE";
780 std::cout <<
"\t[DONE]\n";
787 Bool_t StRefMultCorr::isIndexOk()
const {
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;
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()));
811 Bool_t StRefMultCorr::isZvertexOk()
const {
813 return (mVz > mStart_zvertex[mParameterIndex] &&
814 mVz < mStop_zvertex[mParameterIndex]);
818 Bool_t StRefMultCorr::isRefMultOk()
const {
820 if ( !isIndexOk() )
return kFALSE ;
823 return (mRefMult_corr > mCentrality_bins[0][mParameterIndex] &&
824 mRefMult_corr < mCentrality_bins[mNCentrality][mParameterIndex]);
828 Bool_t StRefMultCorr::isCentralityOk(
const Int_t icent)
const {
830 if ( icent < -1 || icent >= mNCentrality+1 )
return kFALSE ;
833 if ( !isIndexOk() )
return kFALSE ;
837 if ( icent == -1 )
return ( mRefMult_corr <= mCentrality_bins[0][mParameterIndex] );
840 if ( icent == mNCentrality )
return ( mRefMult_corr <= mCentrality_bins[mNCentrality][mParameterIndex] );
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]
856 mParameterIndex = -1 ;
859 setParameterIndex(RunId) ;
863 Int_t StRefMultCorr::setParameterIndex(
const Int_t 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 ;
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 ) {
881 if (mVzEdgeForWeight.empty()) {
882 setVzForWeight(nWeightVzBin_Run14_P16id,
883 WeightVzEdgeMin_Run14_P16id,
884 WeightVzEdgeMax_Run14_P16id);
886 if (mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
887 readScaleForWeight(nWeightgRefmultBin_Run14_P16id,
888 weight_VpdMB5ToVpdMB30_Run14_P16id);
892 else if ( mSubName.CompareTo(
"Run16_AuAu200_VpdMB5", TString::kIgnoreCase) == 0 ) {
893 if(mLibName.CompareTo(
"P16ij", TString::kIgnoreCase) == 0) {
895 if (RunId / 1000 >= 17039 && RunId / 1000 <= 17130) {
897 if (mVzEdgeForWeight.empty()) {
898 setVzForWeight(nWeightVzBin_Run16_P16ij_prod1,
899 WeightVzEdgeMin_Run16_P16ij_prod1,
900 WeightVzEdgeMax_Run16_P16ij_prod1);
902 if (mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
903 readScaleForWeight(nWeightgRefmultBin_Run16_P16ij_prod1,
904 weight_VpdMB5ToVpdMBnoVtx_Run16_P16ij_prod1);
908 else if (RunId / 1000 >= 17169 && RunId / 1000 <= 17179) {
910 if (mVzEdgeForWeight.empty()) {
911 setVzForWeight(nWeightVzBin_Run16_P16ij_prod2,
912 WeightVzEdgeMin_Run16_P16ij_prod2,
913 WeightVzEdgeMax_Run16_P16ij_prod2);
915 if (mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
916 readScaleForWeight(nWeightgRefmultBin_Run16_P16ij_prod2,
917 weight_VpdMB5ToVpdMBnoVtx_Run16_P16ij_prod2);
922 else if ( mSubName.CompareTo(
"Run14_AuAu200_VpdMB30", TString::kIgnoreCase) == 0 ) {
923 if (mLibName.CompareTo(
"P16id", TString::kIgnoreCase) == 0) {
927 else if ( mSubName.CompareTo(
"Run14_AuAu200_VpdMBnoVtx_LowMid", TString::kIgnoreCase) == 0 ) {
928 if(mLibName.CompareTo(
"P16id", TString::kIgnoreCase) == 0) {
932 else if ( mSubName.CompareTo(
"Run14_AuAu200_VpdMBnoVtx_High", TString::kIgnoreCase) == 0 ) {
933 if(mLibName.CompareTo(
"P15ic", TString::kIgnoreCase) == 0) {
937 else if ( mSubName.CompareTo(
"Run16_AuAu200_VpdMBnoVtx", TString::kIgnoreCase) == 0 ) {
938 if(mLibName.CompareTo(
"P16ij", TString::kIgnoreCase) == 0) {
943 mParameterIndex = -1;
948 if ( mName.CompareTo(
"refmult", TString::kIgnoreCase) == 0 &&
949 mSubName.CompareTo(
"Isobar", TString::kIgnoreCase) == 0 ) {
952 mIsZr = mIsRu = kFALSE;
953 for (Int_t i = 0; i < nRunId_Zr; i++) {
954 if (RunId == IsobarRunId_Zr[i]) {
956 mParameterIndex = 36;
961 for (Int_t i = 0; i < nRunId_Ru; i++) {
962 if (RunId == IsobarRunId_Ru[i]) {
964 mParameterIndex = 37;
969 if(!mIsZr && !mIsRu) {
970 Error(
"StRefMultCorr::setParameterIndex",
"RUN %d is not isobaric data", RunId);
976 if(mParameterIndex == -1) {
977 Error(
"StRefMultCorr::setParameterIndex",
"Parameter set does not exist for RUN %d", RunId);
982 std::cout <<
"Parameter index set to: " << mParameterIndex << std::endl;
985 return mParameterIndex ;
991 return mRefMult_corr ;
997 Double_t lumiCorr = 1.;
1000 std::cout <<
"Estimation of luminosity correction factor..." << std::endl;
1001 std::cout <<
"\t ZDC coincidence rate: " << zdcCoincidenceRate <<
" BBC coincidence rate: " << 0 << std::endl;
1006 const Double_t par0_lum = mPar_luminosity[0][mParameterIndex] ;
1007 const Double_t par1_lum = mPar_luminosity[1][mParameterIndex] ;
1009 if( mParameterIndex==36 || mParameterIndex==37 || mParameterIndex==40 ||
1010 ( mRefX==5 && mParameterIndex==7 ) ||
1011 ( mRefX==5 && mParameterIndex==8 ) ) {
1014 Double_t b_prime = 1.;
1015 if(mParameterIndex==36) b_prime = 96.9914;
1016 if(mParameterIndex==37) b_prime = 97.9927;
1017 if(mParameterIndex==40) b_prime = 213.383;
1018 if(mParameterIndex==7 ) b_prime = 106.245;
1019 if(mParameterIndex==8 ) b_prime = 114.041;
1020 lumiCorr = (par0_lum<std::numeric_limits<double>::epsilon() ) ? 1.0 : b_prime/(par0_lum+zdcCoincidenceRate*par1_lum);
1023 lumiCorr = (par0_lum < std::numeric_limits<double>::epsilon() ) ? 1.0 : 1.0/(1.0 + par1_lum/par0_lum*zdcCoincidenceRate/1000.);
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 )
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;
1043 std::cout <<
"\tLuminosity correction factor: " << lumiCorr << std::endl;
1044 std::cout <<
"\t[DONE]\n";
1053 std::cout <<
"Estimation of acceptance correction factor..." << std::endl;
1054 std::cout <<
"\t Vertex z position: " << z <<
" mParameterIndex: " << mParameterIndex << std::endl;
1057 Double_t vzCorr = 1.;
1058 if ( mParameterIndex < 38 ) {
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];
1072 const Double_t par7 = mPar_z_vertex[7][mParameterIndex];
1074 const Double_t RefMult_ref = par0;
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;
1076 if(RefMult_z > 0.0) {
1077 vzCorr = (RefMult_ref + par7)/RefMult_z;
1080 else if ( mParameterIndex == 38 ) {
1082 vzCorr = auau19_run19_vzCorr[ getVzWindowForVzDepCentDef() ];
1084 else if ( mParameterIndex == 39 ) {
1086 vzCorr = auau14_run19_vzCorr[ getVzWindowForVzDepCentDef() ];
1088 else if ( mParameterIndex == 40 ) {
1090 vzCorr = auau200_run19_vzCorr[ getVzWindowForVzDepCentDef() ];
1092 else if ( mParameterIndex == 41 ) {
1094 vzCorr = auau7_run21_vzCorr[ getVzWindowForVzDepCentDef() ];
1096 else if ( mParameterIndex == 42 ) {
1098 vzCorr = auau9_trig2_run20_vzCorr[ getVzWindowForVzDepCentDef() ];
1100 else if ( mParameterIndex == 43 ) {
1102 vzCorr = auau17_run21_vzCorr[ getVzWindowForVzDepCentDef() ];
1104 else if ( mParameterIndex == 44 ) {
1106 vzCorr = auau11_run20_vzCorr[ getVzWindowForVzDepCentDef() ];
1110 std::cout <<
"\t Acceptance correction factor: " << vzCorr << std::endl;
1111 std::cout <<
"\t[DONE]\n";
1120 std::cout <<
"Sampling refMult value: " << refMult << std::endl;
1123 Double_t refMult_d = -9999.;
1124 if( mParameterIndex>=30 && mParameterIndex<=44 ) {
1125 refMult_d = (Double_t)refMult - 0.5 + gRandom->Rndm();
1128 refMult_d = (Double_t)refMult + gRandom->Rndm();
1132 std::cout <<
"\tSampled refMult value: " << refMult_d << std::endl;
1133 std::cout <<
"\t[DONE]\n";
1141 const Double_t zdcCoincidenceRate,
1142 const UInt_t flag)
const {
1145 std::cout <<
"Start refMultCorr calculations" << std::endl
1146 <<
"Initial values refMult / vz / zdcX / flag: "
1147 << refMult <<
" / " << z <<
" / " << zdcCoincidenceRate <<
" / "
1148 << flag << std::endl;
1151 if (!isIndexOk() || !isZvertexOk())
return refMult ;
1158 Double_t refMultCorr = -9999. ;
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;
1164 Error(
"StRefMultCorr::getRefMultCorr",
"invalid flag, flag=%d, should be 0, 1 or 2", flag);
1165 refMultCorr = -9999.;
1170 std::cout <<
"Final refMultCorr value: " << refMultCorr << std::endl;
1177 void StRefMultCorr::readScaleForWeight(
const Char_t* input) {
1178 std::ifstream fin(input) ;
1180 Error(
"StRefMultCorr::readScaleForWeight",
"can't open %s", input);
1185 if(mnVzBinForWeight==0) {
1186 Error(
"StRefMultCorr::readScaleForWeight",
1187 "Please call setVzForWeight() to set vz bin size");
1192 if(!mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
1193 Error(
"StRefMultCorr::readScaleForWeight",
1194 "scale factor has already set in the array");
1198 std::cout <<
"StRefMultCorr::readScaleForWeight Read scale factor ..." << std::flush;
1201 Double_t scale[mnVzBinForWeight] ;
1202 for(Int_t i=0; i<mnVzBinForWeight; i++) {
1206 if(fin.eof())
break ;
1208 for(Int_t i=0; i<mnVzBinForWeight; i++) {
1209 mgRefMultTriggerCorrDiffVzScaleRatio.push_back(scale[i]);
1212 std::cout <<
" [OK]" << std::endl;
1218 void StRefMultCorr::readScaleForWeight(
const Int_t nRefmultbin,
const Double_t *weight) {
1221 if( mnVzBinForWeight==0 ) {
1222 Error(
"StRefMultCorr::readScaleForWeight",
1223 "Please call setVzForWeight() to set vz bin size");
1228 if(!mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
1229 Error(
"StRefMultCorr::readScaleForWeight",
1230 "scale factor has already set in the array");
1234 std::cout <<
"StRefMultCorr::readScaleForWeight Read scale factor ..." << std::flush;
1236 for(Int_t i=0; i<nRefmultbin*mnVzBinForWeight; i++) {
1237 mgRefMultTriggerCorrDiffVzScaleRatio.push_back(weight[i]);
1240 std::cout <<
" [OK]" << std::endl;
1247 void StRefMultCorr::setVzForWeight(
const Int_t nbin,
const Double_t min,
const Double_t max) {
1249 if ( !mVzEdgeForWeight.empty() ) {
1250 Error(
"StRefMultCorr::setVzForWeight",
1251 "z-vertex range for weight has already been defined");
1255 mnVzBinForWeight = nbin ;
1257 const Double_t step = (max-min)/(Double_t)nbin;
1259 for(Int_t i=0; i<mnVzBinForWeight+1; i++) {
1260 mVzEdgeForWeight.push_back( min + step*i );
1264 for(Int_t i=0; i<mnVzBinForWeight; i++) {
1265 std::cout << i <<
" " << step <<
" " << mVzEdgeForWeight[i]
1266 <<
", " << mVzEdgeForWeight[i+1] << std::endl;
1271 Double_t StRefMultCorr::getScaleForWeight()
const {
1276 if(mgRefMultTriggerCorrDiffVzScaleRatio.empty())
return 1.0 ;
1281 Double_t VPD5weight=1.0;
1282 for(Int_t j=0;j<mnVzBinForWeight;j++) {
1283 if(mVz>mVzEdgeForWeight[j] && mVz<=mVzEdgeForWeight[j+1]) {
1296 const Int_t refMultbin=
static_cast<Int_t
>(mRefMult_corr);
1297 VPD5weight=mgRefMultTriggerCorrDiffVzScaleRatio[refMultbin*mnVzBinForWeight + j];
1298 const Double_t tmpContent=VPD5weight;
1301 if(mRefMult_corr>550 && (tmpContent>3.0||tmpContent<0.3)) VPD5weight=1.0;
1306 return 1.0/VPD5weight;
1314 std::cout <<
"Estimating shape correction factor" << std::endl;
1317 Double_t weight = 1.;
1318 Int_t iVzBinIndex = getVzWindowForVzDepCentDef();
1320 if ( mParameterIndex>=30 && mParameterIndex<=35 ) {
1322 if(mRefMult_corr>=500.)
return 1.0;
1324 Int_t iRunPdIndex = mParameterIndex-30;
1325 Int_t iRefmultBin = (Int_t)(mRefMult_corr/2.);
1328 if(iRunPdIndex<0 || iRunPdIndex>5)
return 1.0;
1329 if(iVzBinIndex<0 || iVzBinIndex>13)
return 1.0;
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;
1343 Double_t ShapeReweight = 1.0;
1345 Double_t tem_ShapeReweight = std::stod( sParam_ShapeWeight[iRefmultBin] );
1347 if( tem_ShapeReweight<0.1 ) {
1348 ShapeReweight = 0.1;
1350 else if(tem_ShapeReweight>10) {
1351 ShapeReweight = 10.;
1353 else if(tem_ShapeReweight>=0.1 && tem_ShapeReweight<=10.) {
1354 ShapeReweight = tem_ShapeReweight;
1357 weight = 1./ShapeReweight;
1360 else if ( mParameterIndex>=36 && mParameterIndex<=37 ) {
1362 if (mVz >= -9 && mVz <= 9) {
1366 Int_t mShapeIndex = 0;
1372 if (iVzBinIndex >= 22) {
1373 weight = ShapeWeightArray[mShapeIndex][iVzBinIndex - 9][TMath::Nint(mRefMult_corr)];
1376 weight = ShapeWeightArray[mShapeIndex][iVzBinIndex][TMath::Nint(mRefMult_corr)];
1379 if (weight == 0 || TMath::IsNaN(weight)) {
1383 else if (mParameterIndex == 38) {
1385 if (iVzBinIndex < 0 || iVzBinIndex > auau19_run19_nVzBins)
return 1.0;
1387 weight = auau19_run19_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1389 if (weight == 0 || TMath::IsNaN(weight)) {
1393 else if (mParameterIndex == 39) {
1395 if (iVzBinIndex < 0 || iVzBinIndex > auau14_run19_nVzBins)
return 1.0;
1397 weight = auau14_run19_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1399 if (weight == 0 || TMath::IsNaN(weight)) {
1403 else if (mParameterIndex == 40) {
1405 if (iVzBinIndex < 0 || iVzBinIndex > auau200_run19_nVzBins)
return 1.0;
1407 weight = auau200_run19_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1409 if (weight == 0 || TMath::IsNaN(weight)) {
1413 else if (mParameterIndex == 41) {
1415 if (iVzBinIndex < 0 || iVzBinIndex > auau7_run21_nVzBins)
return 1.0;
1417 weight = auau7_run21_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1419 if (weight == 0 || TMath::IsNaN(weight)) {
1423 else if (mParameterIndex == 42) {
1425 if (iVzBinIndex < 0 || iVzBinIndex > auau9_run20_nVzBins)
return 1.0;
1427 weight = auau9_trig2_run20_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1429 if (weight == 0 || TMath::IsNaN(weight)) {
1433 else if (mParameterIndex == 43) {
1435 if (iVzBinIndex < 0 || iVzBinIndex > auau17_run21_nVzBins)
return 1.0;
1437 weight = auau17_run21_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1439 if (weight == 0 || TMath::IsNaN(weight)) {
1443 else if (mParameterIndex == 44) {
1445 if (iVzBinIndex < 0 || iVzBinIndex > auau11_run20_nVzBins)
return 1.0;
1447 weight = auau11_run20_shapeWeightArray[iVzBinIndex][TMath::Nint(mRefMult_corr)];
1449 if (weight == 0 || TMath::IsNaN(weight)) {
1472 std::cout <<
"\tShape weight: " << weight << std::endl;
1473 std::cout <<
"\t[DONE]\n";
1483 std::cout <<
"Estimating the trigger weight" << std::endl;
1486 Double_t weight = 1.;
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];
1494 const Double_t par6 = mPar_weight[6][mParameterIndex];
1495 const Double_t par7 = mPar_weight[7][mParameterIndex];
1498 std::cout <<
"Trigger efficiency parameters: "
1506 <<
" [7]: " << par7 << std::endl;
1515 && mRefMult_corr < mNormalize_stop[mParameterIndex]
1516 && mRefMult_corr != -(par3/par2) ) {
1519 if (mRefX == 5 && mParameterIndex == 0) {
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) ) );
1537 weight = weight + (weight - 1.0) * (A * mVz * mVz);
1540 std::cout <<
"\tTrigger weight: " << weight << std::endl;
1541 std::cout <<
"\t[DONE]\n";
1549 Double_t Weight = 1.0;
1552 if( !isIndexOk() )
return Weight ;
1554 if( !isZvertexOk() )
return Weight ;
1562 Weight *= getScaleForWeight();
1574 std::cout <<
"Finding centrality bin (out of 16)\n";
1576 Int_t CentBin16 = -1;
1579 if( !isIndexOk() )
return CentBin16 ;
1581 while( CentBin16 < mNCentrality && !isCentralityOk(CentBin16) ) {
1587 if( mParameterIndex==28&&(mVz<-27||mVz>25) ) {
1588 CentBin16 = getCentralityBin16VzDep();
1595 std::cout <<
"\tCentrality bin (out of 16): " << CentBin16 <<
"\t[DONE]" << std::endl;
1599 return ( CentBin16==16 ) ? -1 : CentBin16;
1606 std::cout <<
"Finding centrality bin (out of 9)\n";
1609 Int_t CentBin9 = -1;
1612 if ( !isIndexOk() )
return CentBin9 ;
1615 const Bool_t isCentralityOk = CentBin16 >= 0 && CentBin16 < mNCentrality ;
1618 if ( !isCentralityOk )
return CentBin9;
1621 if( mRefMult_corr > mCentrality_bins[15][mParameterIndex] &&
1622 mRefMult_corr <= mCentrality_bins[16][mParameterIndex] ) {
1625 else if( mRefMult_corr > mCentrality_bins[14][mParameterIndex] &&
1626 mRefMult_corr <= mCentrality_bins[15][mParameterIndex]) {
1630 CentBin9 = (Int_t)(0.5*CentBin16);
1636 if ( mParameterIndex == 28 && ( mVz<-27 || mVz>25 ) ) {
1637 CentBin9 = getCentralityBin9VzDep();
1641 std::cout <<
"\tCentrality bin (out of 9): " << CentBin9 <<
"\t[DONE]" << std::endl;
1648 Int_t StRefMultCorr::getVzWindowForVzDepCentDef()
const {
1651 if( mParameterIndex==28 ) {
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;
1659 else if( mParameterIndex>=30 && mParameterIndex<=35 ) {
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;
1676 else if( mParameterIndex>=36 && mParameterIndex<=37 ) {
1677 Double_t VtxZBinDouble = mVz/2. + 17.;
1682 else if (mVz == -35.) {
1686 iBinVz = TMath::Nint(VtxZBinDouble);
1689 else if ( (mParameterIndex == 38) || (mParameterIndex == 39) ) {
1691 for ( Int_t iVz=0; iVz<auau19_run19_nVzBins; iVz++ ) {
1692 if ( auau19_run19_vzRangeLimits[iVz][0] <= mVz && mVz < auau19_run19_vzRangeLimits[iVz][1] ) {
1698 else if ( mParameterIndex == 40 ) {
1699 for ( Int_t iVz=0; iVz<auau200_run19_nVzBins; iVz++ ) {
1700 if ( auau200_run19_vzRangeLimits[iVz][0] <= mVz && mVz < auau200_run19_vzRangeLimits[iVz][1] ) {
1706 else if ( mParameterIndex == 41 ) {
1707 for ( Int_t iVz=0; iVz<auau7_run21_nVzBins; iVz++ ) {
1708 if ( auau7_run21_vzRangeLimits[iVz][0] <= mVz && mVz < auau7_run21_vzRangeLimits[iVz][1] ) {
1714 else if ( mParameterIndex == 42 ) {
1715 for ( Int_t iVz=0; iVz<auau9_run20_nVzBins; iVz++ ) {
1716 if ( auau9_run20_vzRangeLimits[iVz][0] <= mVz && mVz < auau9_run20_vzRangeLimits[iVz][1] ) {
1722 else if ( mParameterIndex == 43 ) {
1723 for ( Int_t iVz=0; iVz<auau17_run21_nVzBins; iVz++ ) {
1724 if ( auau17_run21_vzRangeLimits[iVz][0] <= mVz && mVz < auau17_run21_vzRangeLimits[iVz][1] ) {
1730 else if ( mParameterIndex == 44 ) {
1731 for ( Int_t iVz=0; iVz<auau11_run20_nVzBins; iVz++ ) {
1732 if ( auau11_run20_vzRangeLimits[iVz][0] <= mVz && mVz < auau11_run20_vzRangeLimits[iVz][1] ) {
1756 Int_t StRefMultCorr::getCentralityBin9VzDep()
const {
1757 const Int_t vzid = getVzWindowForVzDepCentDef();
1759 for(Int_t i=0; i<9; i++) {
1761 if( mRefMult_corr>CentBin9_vzdep[vzid][i] && mRefMult_corr<50000 ) iCent = i;
1763 else if( mRefMult_corr>CentBin9_vzdep[vzid][i] &&
1764 mRefMult_corr<CentBin9_vzdep[vzid][i+1] ) iCent = i;
1767 if( mRefMult_corr>0 && mRefMult_corr<CentBin9_vzdep[vzid][0] ) iCent = -1;
1768 return ( iCent==9999 ) ? -1 : iCent;
1772 Int_t StRefMultCorr::getCentralityBin16VzDep()
const {
1773 const Int_t vzid = getVzWindowForVzDepCentDef();
1775 for(Int_t i=0; i<16; i++) {
1777 if (mRefMult_corr > CentBin16_vzdep[vzid][i] && mRefMult_corr < 50000) {
1781 else if (mRefMult_corr > CentBin16_vzdep[vzid][i] &&
1782 mRefMult_corr < CentBin16_vzdep[vzid][i + 1]) {
1787 if( mRefMult_corr>0 && mRefMult_corr<CentBin16_vzdep[vzid][0] ) {
1790 return (iCent==9999) ? -1 : iCent;
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;
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;
1818 void StRefMultCorr::readHeaderFile() {
1824 const Int_t refX = getRefX();
1827 const Int_t nID = getNumberOfDatasets();
1829 for(
int iID=0; iID<nID; iID++) {
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],
',');
1840 Int_t startRunId=0, stopRunId=0;
1841 startRunId = stoi(sRuns[0]);
1842 stopRunId = stoi(sRuns[1]);
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]);
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 ) ;
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]) );
1864 mCentrality_bins[mNCentrality].push_back( 5000 );
1869 Double_t normalize_stop=-1.0 ;
1870 normalize_stop = stod(getParamX(refX,iID,2)) ;
1871 mNormalize_stop.push_back( normalize_stop );
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]);
1881 mPar_z_vertex[i].push_back( val );
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]);
1892 mPar_weight[i].push_back( val );
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]);
1903 mPar_luminosity[i].push_back( val );
1909 std::cout <<
"StRefMultCorr::readHeaderFile [" << mName
1910 <<
"] Correction parameters and centrality definitions have been successfully read for refX=" << refX
1915 void StRefMultCorr::readBadRunsFromHeaderFile() {
1921 for(Int_t i=0; i<nBadRun_refmult_2010; i++) {
1922 mBadRun.push_back(badrun_refmult_2010[i]);
1924 std::cout <<
"read in nBadRun_refmult_2010: " << nBadRun_refmult_2010 << std::endl;
1926 for(Int_t i=0; i<nBadRun_refmult_2011; i++) {
1927 mBadRun.push_back(badrun_refmult_2011[i]);
1929 std::cout <<
"read in nBadRun_refmult_2011: " << nBadRun_refmult_2011 << std::endl;
1931 for(Int_t i=0; i<nBadRun_grefmult_2014; i++) {
1932 mBadRun.push_back(badrun_grefmult_2014[i]);
1934 std::cout <<
"read in nBadRun_grefmult_2014: " << nBadRun_grefmult_2014 << std::endl;
1936 for(Int_t i=0; i<nBadRun_grefmult_2016; i++) {
1937 mBadRun.push_back(badrun_grefmult_2016[i]);
1939 std::cout <<
"read in nBadRun_grefmult_2016: " << nBadRun_grefmult_2016 << std::endl;
1941 for(Int_t i=0; i<nBadRun_refmult_2017; i++) {
1942 mBadRun.push_back(badrun_refmult_2017[i]);
1944 std::cout <<
"read in nBadRun_refmult_2017: " << nBadRun_refmult_2017 << std::endl;
1946 for (Int_t i = 0; i < nBadRun_refmult_2018; i++) {
1947 mBadRun.push_back(badrun_refmult_2018[i]);
1949 std::cout <<
"read in nBadRun_refmult_2018: " << nBadRun_refmult_2018 << std::endl;
1951 for (Int_t i = 0; i < nBadRun_refmult_2019; i++) {
1952 mBadRun.push_back(badrun_refmult_2019[i]);
1954 std::cout <<
"read in nBadRun_refmult_2019: " << nBadRun_refmult_2019 << std::endl;
1956 for (Int_t i = 0; i < nBadRun_refmult_2020; i++) {
1957 mBadRun.push_back(badrun_refmult_2020[i]);
1959 std::cout <<
"read in nBadRun_refmult_2020: " << nBadRun_refmult_2020 << std::endl;
1961 for (Int_t i = 0; i < nBadRun_refmult_2021; i++) {
1962 mBadRun.push_back(badrun_refmult_2021[i]);
1964 std::cout <<
"read in nBadRun_refmult_2021: " << nBadRun_refmult_2021 << std::endl;
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;
1974 std::cout <<
"StRefMultCorr::print Print input parameters for "
1975 << mName <<
" ========================================" << std::endl << std::endl;
1980 for(UInt_t
id=0;
id<mStart_runId.size();
id++) {
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;
1995 std::cout <<
"Centrality: ";
1997 for(Int_t i=0;i<mNCentrality;i++) {
1998 std::cout << Form(
" >%2d%%", 80-5*i);
2000 std::cout << std::endl;
2001 std::cout <<
"RefMult: ";
2003 for(Int_t i=0;i<mNCentrality;i++) {
2005 const TString tmp(
">");
2006 const TString centrality = tmp + Form(
"%d", mCentrality_bins[i][
id]);
2007 std::cout << Form(
"%6s", centrality.Data());
2009 std::cout << std::endl;
2011 for(Int_t i=0;i<mNPar_z_vertex;i++) {
2012 std::cout <<
" mPar_z_vertex[" << i <<
"] = " << mPar_z_vertex[i][id];
2014 std::cout << std::endl;
2016 for(Int_t i=0;i<mNPar_weight;i++) {
2017 std::cout <<
" mPar_weight[" << i <<
"] = " << mPar_weight[i][id];
2019 std::cout << std::endl;
2021 for(Int_t i=0;i<mNPar_luminosity;i++) {
2022 std::cout <<
" mPar_luminosity[" << i <<
"] = " << mPar_luminosity[i][id];
2024 std::cout << std::endl << std::endl;
2026 std::cout <<
"=====================================================================================" << std::endl;
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);
2034 while( getline(ss,buffer,sep) ) {
2035 vstr.push_back(buffer);
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) );
Double_t luminosityCorrection(Double_t zdcCoincidenceRate) const
Luminosity correction factor.
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|<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 -> 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.