StRefMultCorr.cxx 23 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712
  1. //------------------------------------------------------------------------------
  2. // $Id$
  3. // $Log$
  4. // Revision 1.14 2015/05/22 06:52:05 hmasui
  5. // Add grefmult for Run14 Au+Au 200 GeV
  6. //
  7. // Revision 1.13 2013/05/10 18:33:33 hmasui
  8. // Add TOF tray mult, preliminary update for Run12 U+U
  9. //
  10. // Revision 1.12 2012/05/19 00:48:20 hmasui
  11. // Update refmult3
  12. //
  13. // Revision 1.11 2012/05/14 00:40:25 hmasui
  14. // Exit code if no valid run number found. Fix the initialization of mParameterIndex
  15. //
  16. // Revision 1.10 2012/05/09 22:26:46 hmasui
  17. // Commented out option use in the print() function
  18. //
  19. // Revision 1.9 2012/05/08 03:19:49 hmasui
  20. // Move parameters to Centrality_def_refmult.txt
  21. //
  22. // Revision 1.8 2012/04/23 21:29:37 hmasui
  23. // Added isBadRun() function for outlier rejection, getBeginRun() and getEndRun() to obtain the run range for a given (energy,year)
  24. //
  25. // Revision 1.7 2011/11/30 00:25:07 hmasui
  26. // Additional check for ifstream to avoid reading one extra line from input file
  27. //
  28. // Revision 1.6 2011/11/08 19:11:05 hmasui
  29. // Add luminosity corrections for 200 GeV
  30. //
  31. // Revision 1.5 2011/10/11 19:35:20 hmasui
  32. // Fix typo. Add z-vertex check in getWeight() function
  33. //
  34. // Revision 1.4 2011/10/10 21:30:37 hmasui
  35. // Replaced hard coded parameters for z-vertex and weight corrections by input parameters from text file
  36. //
  37. //
  38. // Revision 1.3 2011/08/12 20:28:07 hmasui
  39. // Avoid varying corrected refmult in the same event by random number
  40. //
  41. // Revision 1.2 2011/08/11 23:51:10 hmasui
  42. // Suppress cout in the setParameterIndex function. Use TError for error messages.
  43. //
  44. // Revision 1.1 2011/08/11 18:38:28 hmasui
  45. // First version of Refmult correction class
  46. //
  47. //------------------------------------------------------------------------------
  48. #include <algorithm>
  49. #include <fstream>
  50. #include <iostream>
  51. #include <string>
  52. #include "StRefMultCorr.h"
  53. #include "TError.h"
  54. #include "TRandom.h"
  55. #include "TMath.h"
  56. ClassImp(StRefMultCorr)
  57. using namespace std ;
  58. namespace {
  59. typedef pair<Double_t, Int_t> keys;
  60. }
  61. //______________________________________________________________________________
  62. // Default constructor
  63. StRefMultCorr::StRefMultCorr(const TString name)
  64. : mName(name)
  65. {
  66. mRefMult = 0 ;
  67. mVz = -9999. ;
  68. mRefMult_corr = -1.0 ;
  69. // Clear all data members
  70. clear() ;
  71. // Read parameters
  72. read() ;
  73. readBadRuns() ;
  74. }
  75. //______________________________________________________________________________
  76. // Default destructor
  77. StRefMultCorr::~StRefMultCorr()
  78. {
  79. }
  80. //______________________________________________________________________________
  81. Int_t StRefMultCorr::getBeginRun(const Double_t energy, const Int_t year)
  82. {
  83. keys key(std::make_pair(energy, year));
  84. // Make sure key exists
  85. multimap<keys, Int_t>::iterator iterCheck = mBeginRun.find(key);
  86. if ( iterCheck == mBeginRun.end() ) {
  87. Error("StRefMultCorr::getBeginRun", "can't find energy = %1.1f, year = %d", energy, year);
  88. return -1;
  89. }
  90. pair<multimap<keys, Int_t>::iterator, multimap<keys, Int_t>::iterator> iterRange = mBeginRun.equal_range(key);
  91. return (*(iterRange.first)).second ;
  92. }
  93. //______________________________________________________________________________
  94. Int_t StRefMultCorr::getEndRun(const Double_t energy, const Int_t year)
  95. {
  96. keys key(std::make_pair(energy, year));
  97. // Make sure key exists
  98. multimap<keys, Int_t>::iterator iterCheck = mEndRun.find(key);
  99. if ( iterCheck == mEndRun.end() ) {
  100. Error("StRefMultCorr::getEndRun", "can't find energy = %1.1f, year = %d", energy, year);
  101. return -1;
  102. }
  103. pair<multimap<keys, Int_t>::iterator, multimap<keys, Int_t>::iterator> iterRange = mEndRun.equal_range(key);
  104. multimap<keys, Int_t>::iterator iter = iterRange.second ;
  105. iter--;
  106. return (*iter).second ;
  107. }
  108. //______________________________________________________________________________
  109. void StRefMultCorr::clear()
  110. {
  111. // Clear all arrays, and set parameter index = -1
  112. mYear.clear() ;
  113. mStart_runId.clear() ;
  114. mStop_runId.clear() ;
  115. mStart_zvertex.clear() ;
  116. mStop_zvertex.clear() ;
  117. mNormalize_stop.clear() ;
  118. for(Int_t i=0;i<mNCentrality;i++) {
  119. mCentrality_bins[i].clear() ;
  120. }
  121. mParameterIndex = -1 ;
  122. for(Int_t i=0;i<mNPar_z_vertex;i++) {
  123. mPar_z_vertex[i].clear() ;
  124. }
  125. for(Int_t i=0;i<mNPar_weight;i++) {
  126. mPar_weight[i].clear();
  127. }
  128. for(Int_t i=0;i<mNPar_luminosity;i++) {
  129. mPar_luminosity[i].clear();
  130. }
  131. mBeginRun.clear() ;
  132. mEndRun.clear() ;
  133. mBadRun.clear() ;
  134. mnVzBinForWeight = 0 ;
  135. mVzEdgeForWeight.clear();
  136. mgRefMultTriggerCorrDiffVzScaleRatio.clear() ;
  137. }
  138. //______________________________________________________________________________
  139. Bool_t StRefMultCorr::isBadRun(const Int_t RunId)
  140. {
  141. // Return true if a given run id is bad run
  142. vector<Int_t>::iterator iter = std::find(mBadRun.begin(), mBadRun.end(), RunId);
  143. #if 0
  144. if ( iter != mBadRun.end() ) {
  145. // QA
  146. cout << "StRefMultCorr::isBadRun Find bad run = " << (*iter) << endl;
  147. }
  148. #endif
  149. return ( iter != mBadRun.end() ) ;
  150. }
  151. //______________________________________________________________________________
  152. void StRefMultCorr::initEvent(const UShort_t RefMult, const Double_t z, const Double_t zdcCoincidenceRate)
  153. {
  154. // Set refmult, vz and corrected refmult if current (refmult,vz) are different from inputs
  155. // User must call this function event-by-event before
  156. // calling any other public functions
  157. if ( mRefMult != RefMult || mVz != z || mZdcCoincidenceRate != zdcCoincidenceRate ) {
  158. mRefMult = RefMult ;
  159. mVz = z ;
  160. mZdcCoincidenceRate = zdcCoincidenceRate ;
  161. mRefMult_corr = getRefMultCorr(mRefMult, mVz, mZdcCoincidenceRate) ;
  162. }
  163. }
  164. //______________________________________________________________________________
  165. Bool_t StRefMultCorr::isIndexOk() const
  166. {
  167. // mParameterIndex not initialized (-1)
  168. if ( mParameterIndex == -1 ) {
  169. Error("StRefMultCorr::isIndexOk", "mParameterIndex = -1. Call init(const Int_t RunId) function to initialize centrality bins, corrections");
  170. Error("StRefMultCorr::isIndexOk", "mParameterIndex = -1. or use valid run numbers defined in Centrality_def_%s.txt", mName.Data());
  171. Error("StRefMultCorr::isIndexOk", "mParameterIndex = -1. exit");
  172. cout << endl;
  173. // Stop the process if invalid run number found
  174. exit(0);
  175. }
  176. // Out of bounds
  177. if ( mParameterIndex >= (Int_t)mStart_runId.size() ) {
  178. Error("StRefMultCorr::isIndexOk",
  179. Form("mParameterIndex = %d > max number of parameter set = %d. Make sure you put correct index for this energy",
  180. mParameterIndex, mStart_runId.size()));
  181. return kFALSE ;
  182. }
  183. return kTRUE ;
  184. }
  185. //______________________________________________________________________________
  186. Bool_t StRefMultCorr::isZvertexOk() const
  187. {
  188. // Primary z-vertex check
  189. return ( mVz > mStart_zvertex[mParameterIndex] && mVz < mStop_zvertex[mParameterIndex] ) ;
  190. }
  191. //______________________________________________________________________________
  192. Bool_t StRefMultCorr::isRefMultOk() const
  193. {
  194. // Invalid index
  195. if ( !isIndexOk() ) return kFALSE ;
  196. // select 0-80%
  197. return (mRefMult_corr > mCentrality_bins[0][mParameterIndex] && mRefMult_corr < mCentrality_bins[mNCentrality][mParameterIndex]);
  198. }
  199. //______________________________________________________________________________
  200. Bool_t StRefMultCorr::isCentralityOk(const Int_t icent) const
  201. {
  202. // Invalid centrality id
  203. if ( icent < -1 || icent >= mNCentrality+1 ) return kFALSE ;
  204. // Invalid index
  205. if ( !isIndexOk() ) return kFALSE ;
  206. // Special case
  207. // 1. 80-100% for icent=-1
  208. if ( icent == -1 ) return (mRefMult_corr <= mCentrality_bins[0][mParameterIndex]);
  209. // 2. icent = mNCentrality
  210. if ( icent == mNCentrality ) return (mRefMult_corr <= mCentrality_bins[mNCentrality][mParameterIndex]);
  211. const Bool_t ok = (mRefMult_corr > mCentrality_bins[icent][mParameterIndex] && mRefMult_corr <= mCentrality_bins[icent+1][mParameterIndex]);
  212. return ok ;
  213. }
  214. //______________________________________________________________________________
  215. void StRefMultCorr::init(const Int_t RunId)
  216. {
  217. // Reset mParameterIndex
  218. mParameterIndex = -1 ;
  219. // call setParameterIndex
  220. setParameterIndex(RunId) ;
  221. }
  222. //______________________________________________________________________________
  223. Int_t StRefMultCorr::setParameterIndex(const Int_t RunId)
  224. {
  225. // Determine the corresponding parameter set for the input RunId
  226. for(UInt_t npar = 0; npar < mStart_runId.size(); npar++)
  227. {
  228. if(RunId >= mStart_runId[npar] && RunId <= mStop_runId[npar])
  229. {
  230. mParameterIndex = npar ;
  231. // cout << "StRefMultCorr::setParameterIndex Parameter set = " << mParameterIndex << " for RUN " << RunId << endl;
  232. break ;
  233. }
  234. }
  235. if(mParameterIndex == -1){
  236. Error("StRefMultCorr::setParameterIndex", "Parameter set does not exist for RUN %d", RunId);
  237. }
  238. //else cout << "Parameter set = " << npar_set << endl;
  239. return mParameterIndex ;
  240. }
  241. //______________________________________________________________________________
  242. Double_t StRefMultCorr::getRefMultCorr() const
  243. {
  244. // Call initEvent() first
  245. return mRefMult_corr ;
  246. }
  247. //______________________________________________________________________________
  248. Double_t StRefMultCorr::getRefMultCorr(const UShort_t RefMult, const Double_t z,
  249. const Double_t zdcCoincidenceRate, const UInt_t flag) const
  250. {
  251. // Apply correction if parameter index & z-vertex are ok
  252. if (!isIndexOk() || !isZvertexOk()) return RefMult ;
  253. // Correction function for RefMult, takes into account z_vertex dependence
  254. // Luminosity corrections
  255. // 200 GeV only. correction = 1 for all the other energies
  256. const Double_t par0l = mPar_luminosity[0][mParameterIndex] ;
  257. const Double_t par1l = mPar_luminosity[1][mParameterIndex] ;
  258. const Double_t correction_luminosity = (par0l==0.0) ? 1.0 : 1.0/(1.0 + par1l/par0l*zdcCoincidenceRate/1000.);
  259. // par0 to par5 define the parameters of a polynomial to parametrize z_vertex dependence of RefMult
  260. const Double_t par0 = mPar_z_vertex[0][mParameterIndex];
  261. const Double_t par1 = mPar_z_vertex[1][mParameterIndex];
  262. const Double_t par2 = mPar_z_vertex[2][mParameterIndex];
  263. const Double_t par3 = mPar_z_vertex[3][mParameterIndex];
  264. const Double_t par4 = mPar_z_vertex[4][mParameterIndex];
  265. const Double_t par5 = mPar_z_vertex[5][mParameterIndex];
  266. const Double_t par6 = mPar_z_vertex[6][mParameterIndex];
  267. const Double_t par7 = mPar_z_vertex[7][mParameterIndex]; // this parameter is usually 0, it takes care for an additional efficiency, usually difference between phase A and phase B parameter 0
  268. const Double_t RefMult_ref = par0; // Reference mean RefMult at z=0
  269. 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
  270. Double_t Hovno = 1.0; // Correction factor for RefMult, takes into account z_vertex dependence
  271. if(RefMult_z > 0.0)
  272. {
  273. Hovno = (RefMult_ref + par7)/RefMult_z;
  274. }
  275. Double_t RefMult_d = (Double_t)(RefMult)+gRandom->Rndm(); // random sampling over bin width -> avoid peak structures in corrected distribution
  276. Double_t RefMult_corr = -9999. ;
  277. switch ( flag ) {
  278. case 0: return RefMult_d*correction_luminosity;
  279. case 1: return RefMult_d*Hovno;
  280. case 2: return RefMult_d*Hovno*correction_luminosity;
  281. default:
  282. {
  283. Error("StRefMultCorr::getRefMultCorr", "invalid flag, flag=%d, should be 0,1 or 2", flag);
  284. return -9999.;
  285. }
  286. }
  287. // cout << "Input RefMult = " << RefMult << ", input z = " << z << ", RefMult_corr = " << RefMult_corr << endl;
  288. return RefMult_corr ;
  289. }
  290. //______________________________________________________________________________
  291. void StRefMultCorr::readScaleForWeight(const Char_t* input)
  292. {
  293. ifstream fin(input) ;
  294. if(!fin) {
  295. Error("StRefMultCorr::readScaleForWeight", "can't open %s", input);
  296. return;
  297. }
  298. // Users must set the vz bin size by setVzForWeight() (see below)
  299. if(mnVzBinForWeight==0) {
  300. Error("StRefMultCorr::readScaleForWeight",
  301. "Please call setVzForWeight() to set vz bin size");
  302. return;
  303. }
  304. // Do not allow multiple calls
  305. if(!mgRefMultTriggerCorrDiffVzScaleRatio.empty()) {
  306. Error("StRefMultCorr::readScaleForWeight",
  307. "scale factor has already set in the array");
  308. return;
  309. }
  310. cout << "StRefMultCorr::readScaleForWeight Read scale factor ..."
  311. << flush;
  312. while(fin) {
  313. Double_t scale[mnVzBinForWeight] ;
  314. for(Int_t i=0; i<mnVzBinForWeight; i++) {
  315. fin >> scale[i] ;
  316. }
  317. if(fin.eof()) break ;
  318. for(Int_t i=0; i<mnVzBinForWeight; i++) {
  319. mgRefMultTriggerCorrDiffVzScaleRatio.push_back(scale[i]);
  320. }
  321. }
  322. cout << " [OK]" << endl;
  323. }
  324. //______________________________________________________________________________
  325. void StRefMultCorr::setVzForWeight(const Int_t nbin, const Double_t min,
  326. const Double_t max)
  327. {
  328. // Do not allow multiple calls
  329. if(!mVzEdgeForWeight.empty()) {
  330. Error("StRefMultCorr::setVzForWeight",
  331. "z-vertex range for weight has already been defined");
  332. return;
  333. }
  334. mnVzBinForWeight = nbin ;
  335. // calculate increment size
  336. const Double_t step = (max-min)/(Double_t)nbin;
  337. for(Int_t i=0; i<mnVzBinForWeight+1; i++) {
  338. mVzEdgeForWeight.push_back( min + step*i );
  339. }
  340. // Debug
  341. for(Int_t i=0; i<mnVzBinForWeight; i++) {
  342. cout << i << " " << step << " " << mVzEdgeForWeight[i]
  343. << ", " << mVzEdgeForWeight[i+1] << endl;
  344. }
  345. }
  346. //______________________________________________________________________________
  347. Double_t StRefMultCorr::getScaleForWeight() const
  348. {
  349. // Special scale factor for global refmult in Run14
  350. // to account for the difference between
  351. // VPDMB-30 and VPDMB-5
  352. // return 1 if mgRefMultTriggerCorrDiffVzScaleRatio array is empty
  353. if(mgRefMultTriggerCorrDiffVzScaleRatio.empty()) return 1.0 ;
  354. Double_t VPD5weight=1.0;
  355. for(Int_t j=0;j<mnVzBinForWeight;j++) {
  356. if(mVz>mVzEdgeForWeight[j] && mVz<=mVzEdgeForWeight[j+1]) {
  357. const Int_t refMultbin=static_cast<Int_t>(mRefMult_corr);
  358. VPD5weight=mgRefMultTriggerCorrDiffVzScaleRatio[refMultbin*mnVzBinForWeight + j];
  359. const Double_t tmpContent=VPD5weight;
  360. if(tmpContent==0 || (mRefMult_corr>500 && tmpContent<=0.65)) VPD5weight=1.15;//Just because the value of the weight is around 1.15
  361. if(mRefMult_corr>500 && tmpContent>=1.35) VPD5weight=1.15;//Remove those Too large weight factor,gRefmult>500
  362. // this weight and reweight should be careful, after reweight(most peripheral),Then weight(whole range)
  363. }
  364. }
  365. return 1.0/VPD5weight;
  366. }
  367. //______________________________________________________________________________
  368. Double_t StRefMultCorr::getWeight() const
  369. {
  370. Double_t Weight = 1.0;
  371. // Invalid index
  372. if( !isIndexOk() ) return Weight ;
  373. // Invalid z-vertex
  374. if( !isZvertexOk() ) return Weight ;
  375. const Double_t par0 = mPar_weight[0][mParameterIndex];
  376. const Double_t par1 = mPar_weight[1][mParameterIndex];
  377. const Double_t par2 = mPar_weight[2][mParameterIndex];
  378. const Double_t par3 = mPar_weight[3][mParameterIndex];
  379. const Double_t par4 = mPar_weight[4][mParameterIndex];
  380. const Double_t A = mPar_weight[5][mParameterIndex];
  381. const Double_t par6 = mPar_weight[6][mParameterIndex];//Add by guannan for run14
  382. const Double_t par7 = mPar_weight[7][mParameterIndex];//Add by guannan for run14
  383. // Additional z-vetex dependent correction
  384. //const Double_t A = ((1.27/1.21))/(30.0*30.0); // Don't ask...
  385. //const Double_t A = (0.05/0.21)/(30.0*30.0); // Don't ask...
  386. if(isRefMultOk() // 0-80%
  387. && mRefMult_corr < mNormalize_stop[mParameterIndex] // re-weighting only apply up to normalization point
  388. && mRefMult_corr != -(par3/par2) // avoid denominator = 0
  389. )
  390. {
  391. Weight = par0 + par1/(par2*mRefMult_corr + par3) + par4*(par2*mRefMult_corr + par3) + par6/TMath::Power(par2*mRefMult_corr+par3 ,2) + par7*TMath::Power(par2*mRefMult_corr+par3 ,2); // Parametrization of MC/data RefMult ratio
  392. Weight = Weight + (Weight-1.0)*(A*mVz*mVz); // z-dependent weight correction
  393. }
  394. // Special scale factor for global refmult
  395. // for others, scale factor = 1
  396. const Double_t scale = getScaleForWeight() ;
  397. Weight *= scale ;
  398. return Weight ;
  399. }
  400. //______________________________________________________________________________
  401. Int_t StRefMultCorr::getCentralityBin16() const
  402. {
  403. Int_t CentBin16 = -1;
  404. // Invalid index
  405. if( !isIndexOk() ) return CentBin16 ;
  406. while(CentBin16 < mNCentrality && !isCentralityOk(CentBin16) )
  407. {
  408. CentBin16++;
  409. }
  410. // return -1 if CentBin16 = 16 (very large refmult, refmult>5000)
  411. return (CentBin16==16) ? -1 : CentBin16;
  412. }
  413. //______________________________________________________________________________
  414. Int_t StRefMultCorr::getCentralityBin9() const
  415. {
  416. Int_t CentBin9 = -1;
  417. // Invalid index
  418. if ( !isIndexOk() ) return CentBin9 ;
  419. const Int_t CentBin16 = getCentralityBin16(); // Centrality bin 16
  420. const Bool_t isCentralityOk = CentBin16 >= 0 && CentBin16 < mNCentrality ;
  421. // No centrality is defined
  422. if (!isCentralityOk) return CentBin9 ;
  423. // First handle the exceptions
  424. if(mRefMult_corr > mCentrality_bins[15][mParameterIndex] && mRefMult_corr <= mCentrality_bins[16][mParameterIndex])
  425. {
  426. CentBin9 = 8; // most central 5%
  427. }
  428. else if(mRefMult_corr > mCentrality_bins[14][mParameterIndex] && mRefMult_corr <= mCentrality_bins[15][mParameterIndex])
  429. {
  430. CentBin9 = 7; // most central 5-10%
  431. }
  432. else
  433. {
  434. CentBin9 = (Int_t)(0.5*CentBin16);
  435. }
  436. return CentBin9;
  437. }
  438. //______________________________________________________________________________
  439. const Char_t* StRefMultCorr::getTable() const
  440. {
  441. if ( mName.CompareTo("refmult", TString::kIgnoreCase) == 0 ) {
  442. return "StRoot/StRefMultCorr/Centrality_def_refmult.txt";
  443. }
  444. else if ( mName.CompareTo("refmult2", TString::kIgnoreCase) == 0 ) {
  445. return "StRoot/StRefMultCorr/Centrality_def_refmult2.txt";
  446. }
  447. else if ( mName.CompareTo("refmult3", TString::kIgnoreCase) == 0 ) {
  448. return "StRoot/StRefMultCorr/Centrality_def_refmult3.txt";
  449. }
  450. else if ( mName.CompareTo("refmult4", TString::kIgnoreCase) == 0 ) {
  451. return "StRoot/StRefMultCorr/Centrality_def_refmult4.txt";
  452. }
  453. else if ( mName.CompareTo("toftray", TString::kIgnoreCase) == 0 ) {
  454. return "StRoot/StRefMultCorr/Centrality_def_toftray.txt";
  455. }
  456. else if ( mName.CompareTo("grefmult", TString::kIgnoreCase) == 0 ) {
  457. return "StRoot/StRefMultCorr/Centrality_def_grefmult.txt";
  458. }
  459. else{
  460. Error("StRefMultCorr::getTable", "No implementation for %s", mName.Data());
  461. cout << "Current available option is refmult or refmult2 or refmult3 or toftray" << endl;
  462. return "";
  463. }
  464. }
  465. //______________________________________________________________________________
  466. void StRefMultCorr::read()
  467. {
  468. // Open the parameter file and read the data
  469. const Char_t* inputFileName(getTable());
  470. ifstream ParamFile(inputFileName);
  471. if(!ParamFile){
  472. Error("StRefMultCorr::read", "cannot open %s", inputFileName);
  473. return;
  474. }
  475. cout << "StRefMultCorr::read Open " << inputFileName << flush ;
  476. string line ;
  477. getline(ParamFile,line);
  478. if(line.find("Start_runId")!=string::npos)
  479. {
  480. while(ParamFile.good())
  481. {
  482. Int_t year;
  483. Double_t energy;
  484. ParamFile >> year >> energy ;
  485. Int_t startRunId=0, stopRunId=0 ;
  486. Double_t startZvertex=-9999., stopZvertex=-9999. ;
  487. ParamFile >> startRunId >> stopRunId >> startZvertex >> stopZvertex ;
  488. // Error check
  489. if(ParamFile.eof()) break;
  490. mYear.push_back(year) ;
  491. mBeginRun.insert(std::make_pair(std::make_pair(energy, year), startRunId));
  492. mEndRun.insert(std::make_pair(std::make_pair(energy, year), stopRunId));
  493. mStart_runId.push_back( startRunId ) ;
  494. mStop_runId.push_back( stopRunId ) ;
  495. mStart_zvertex.push_back( startZvertex ) ;
  496. mStop_zvertex.push_back( stopZvertex ) ;
  497. for(Int_t i=0;i<mNCentrality;i++) {
  498. Int_t centralitybins=-1;
  499. ParamFile >> centralitybins;
  500. mCentrality_bins[i].push_back( centralitybins );
  501. }
  502. Double_t normalize_stop=-1.0 ;
  503. ParamFile >> normalize_stop ;
  504. mNormalize_stop.push_back( normalize_stop );
  505. for(Int_t i=0;i<mNPar_z_vertex;i++) {
  506. Double_t param=-9999.;
  507. ParamFile >> param;
  508. mPar_z_vertex[i].push_back( param );
  509. }
  510. for(Int_t i=0;i<mNPar_weight;i++) {
  511. Double_t param=-9999.;
  512. ParamFile >> param;
  513. mPar_weight[i].push_back( param );
  514. }
  515. for(Int_t i=0;i<mNPar_luminosity;i++) {
  516. Double_t param=-9999.;
  517. ParamFile >> param;
  518. mPar_luminosity[i].push_back( param );
  519. }
  520. mCentrality_bins[mNCentrality].push_back( 5000 );
  521. }
  522. }
  523. else
  524. {
  525. cout << endl;
  526. Error("StRefMultCorr::read", "Input file is not correct! Wrong structure.");
  527. return;
  528. }
  529. ParamFile.close();
  530. cout << " [OK]" << endl;
  531. }
  532. //______________________________________________________________________________
  533. void StRefMultCorr::readBadRuns()
  534. {
  535. // Read bad run numbers
  536. // - From year 2010 - 2014
  537. // - If input file doesn't exist, skip to the next year without warning
  538. vector<int> years = { 2010, 2011, 2014 };
  539. for(int i = 0; i < years.size(); i++ ) {
  540. cout << "StRefMultCorr::readBadRuns For " << mName << ": open " << flush ;
  541. const Int_t year = years[ i ];
  542. const Char_t* inputFileName(Form("StRoot/StRefMultCorr/bad_runs_refmult_year%d.txt", year));
  543. ifstream fin(inputFileName);
  544. if(!fin){
  545. cout << "Cannot load " << inputFileName << endl;
  546. cout << endl;
  547. continue;
  548. }
  549. cout << " " << inputFileName << flush;
  550. Int_t runId = 0 ;
  551. while( fin >> runId ) {
  552. mBadRun.push_back(runId);
  553. }
  554. cout << " [OK]" << endl;
  555. }
  556. }
  557. //______________________________________________________________________________
  558. void StRefMultCorr::print(const Option_t* option) const
  559. {
  560. cout << "StRefMultCorr::print Print input parameters for " << mName << " ========================================" << endl << endl;
  561. // Option switched off, can be used to specify parameters
  562. // const TString opt(option);
  563. // Int_t input_counter = 0;
  564. for(UInt_t id=0; id<mStart_runId.size(); id++) {
  565. //cout << "Data line = " << input_counter << ", Start_runId = " << Start_runId[input_counter] << ", Stop_runId = " << Stop_runId[input_counter] << endl;
  566. // const UInt_t id = mStart_runId.size()-1;
  567. // Removed line break
  568. cout << " Index=" << id;
  569. cout << Form(" Run=[%8d, %8d]", mStart_runId[id], mStop_runId[id]);
  570. cout << Form(" z-vertex=[%1.1f, %1.1f]", mStart_zvertex[id], mStop_zvertex[id]);
  571. cout << ", Normalize_stop=" << mNormalize_stop[id];
  572. cout << endl;
  573. // if(opt.IsWhitespace()){
  574. // continue ;
  575. // }
  576. cout << "Centrality: ";
  577. for(Int_t i=0;i<mNCentrality;i++){
  578. cout << Form(" >%2d%%", 80-5*i);
  579. }
  580. cout << endl;
  581. cout << "RefMult: ";
  582. for(Int_t i=0;i<mNCentrality;i++){
  583. // cout << Form("StRefMultCorr::read Centrality %3d-%3d %%, refmult > %4d", 75-5*i, 80-5*i, mCentrality_bins[i][id]) << endl;
  584. const TString tmp(">");
  585. const TString centrality = tmp + Form("%d", mCentrality_bins[i][id]);
  586. cout << Form("%6s", centrality.Data());
  587. }
  588. cout << endl;
  589. for(Int_t i=0;i<mNPar_z_vertex;i++) {
  590. cout << " mPar_z_vertex[" << i << "] = " << mPar_z_vertex[i][id];
  591. }
  592. cout << endl;
  593. for(Int_t i=0;i<mNPar_weight;i++) {
  594. cout << " mPar_weight[" << i << "] = " << mPar_weight[i][id];
  595. }
  596. cout << endl;
  597. for(Int_t i=0;i<mNPar_luminosity;i++) {
  598. cout << " mPar_luminosity[" << i << "] = " << mPar_luminosity[i][id];
  599. }
  600. cout << endl << endl;
  601. }
  602. cout << "=====================================================================================" << endl;
  603. }