MpdFemtoModelBPLCMS3DCorrFctnKt.cxx 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452
  1. //
  2. // Three-dimensional Bertsch-Pratt correlation function in LCMS for the model estimations
  3. //
  4. // MpdFemtoMaker headers
  5. #include "MpdFemtoModelGausLCMSFreezeOutGenerator.h"
  6. #include "MpdFemtoModelHiddenInfo.h"
  7. #include "MpdFemtoModelManager.h"
  8. #include "MpdFemtoPair.h"
  9. // MpdFemtoMakerUser headers
  10. #include "MpdFemtoModelBPLCMS3DCorrFctnKt.h"
  11. // ROOT headers
  12. #include "TMath.h"
  13. #include "TString.h"
  14. #include "TH3F.h"
  15. // C++ headers
  16. #include <cmath>
  17. ClassImp(MpdFemtoModelBPLCMS3DCorrFctnKt)
  18. //_________________
  19. MpdFemtoModelBPLCMS3DCorrFctnKt::MpdFemtoModelBPLCMS3DCorrFctnKt(const char* title, const int& nBins,
  20. const double& qLo, const double& qHi,
  21. const int& ktBins, const double& ktLo,
  22. const double& ktHi, const bool isUseDenominator) :
  23. MpdFemtoBaseCorrFctn(),
  24. mManager(nullptr) {
  25. // Parametrized constructor
  26. // Set general parameters
  27. setKtRange( ktBins, ktLo, ktHi );
  28. setUseDenominator( isUseDenominator );
  29. // Define string parameters
  30. TString baseName(title);
  31. TString baseTitle(title);
  32. TString appendix = ";q_{out} (GeV/c);q_{side} (GeV/c);q_{long} (GeV/c)";
  33. // Loop over kTBins
  34. for (int iKtBin = 0; iKtBin < mNKtBins; iKtBin++) {
  35. // Set numerator
  36. baseName = title;
  37. baseName += "_num_";
  38. baseName += iKtBin;
  39. baseTitle = baseName;
  40. baseTitle += " ";
  41. baseTitle += Form("%3.2f",( mKtRange[0] + mKtStep * iKtBin) );
  42. baseTitle += "#leq k_{T} (GeV/c) #leq";
  43. baseTitle += Form("%3.2f",( mKtRange[0] + mKtStep * (iKtBin+1) ) );
  44. baseTitle += appendix;
  45. TH3F *histoNum = new TH3F(baseName.Data(), baseTitle.Data(),
  46. nBins, qLo, qHi,
  47. nBins, qLo, qHi,
  48. nBins, qLo, qHi);
  49. histoNum->Sumw2();
  50. mNumerator.push_back( histoNum );
  51. // Weighted numerator
  52. baseName = title;
  53. baseName += "_num_wei_";
  54. baseName += iKtBin;
  55. baseTitle = baseName;
  56. baseTitle += " ";
  57. baseTitle += Form("%3.2f",( mKtRange[0] + mKtStep * iKtBin) );
  58. baseTitle += "#leq kT (GeV/c) #leq";
  59. baseTitle += Form("%3.2f",( mKtRange[0] + mKtStep * (iKtBin+1) ) );
  60. baseTitle += appendix;
  61. TH3F *histoNumW = new TH3F(baseName.Data(), baseTitle.Data(),
  62. nBins, qLo, qHi,
  63. nBins, qLo, qHi,
  64. nBins, qLo, qHi);
  65. histoNumW->Sumw2();
  66. mNumeratorWeighted.push_back( histoNumW );
  67. // Qinv weighted numerator
  68. baseName = title;
  69. baseName += "_num_qinv_";
  70. baseName += iKtBin;
  71. baseTitle = baseName;
  72. baseTitle += " ";
  73. baseTitle += Form("%3.2f",( mKtRange[0] + mKtStep * iKtBin) );
  74. baseTitle += "#leq kT (GeV/c) #leq";
  75. baseTitle += Form("%3.2f",( mKtRange[0] + mKtStep * (iKtBin+1) ) );
  76. baseTitle += appendix;
  77. TH3F *histoNumQinvW = new TH3F(baseName.Data(), baseTitle.Data(),
  78. nBins, qLo, qHi,
  79. nBins, qLo, qHi,
  80. nBins, qLo, qHi);
  81. histoNumQinvW->Sumw2();
  82. mNumeratorQinvWeighted.push_back( histoNumQinvW );
  83. if ( mIsUseDenominator ) {
  84. // Denominator
  85. baseName = title;
  86. baseName += "_den_";
  87. baseName += iKtBin;
  88. baseTitle = baseName;
  89. baseTitle += " ";
  90. baseTitle += Form("%3.2f",( mKtRange[0] + mKtStep * iKtBin) );
  91. baseTitle += "#leq kT (GeV/c) #leq";
  92. baseTitle += Form("%3.2f",( mKtRange[0] + mKtStep * (iKtBin+1) ) );
  93. baseTitle += appendix;
  94. TH3F *histoDen = new TH3F(baseName.Data(), baseTitle.Data(),
  95. nBins, qLo, qHi,
  96. nBins, qLo, qHi,
  97. nBins, qLo, qHi);
  98. histoDen->Sumw2();
  99. mDenominator.push_back( histoDen );
  100. // Weighted denominator
  101. baseName = title;
  102. baseName += "_den_wei_";
  103. baseName += iKtBin;
  104. baseTitle = baseName;
  105. baseTitle += " ";
  106. baseTitle += Form("%3.2f",( mKtRange[0] + mKtStep * iKtBin) );
  107. baseTitle += "#leq kT (GeV/c) #leq";
  108. baseTitle += Form("%3.2f",( mKtRange[0] + mKtStep * (iKtBin+1) ) );
  109. baseTitle += appendix;
  110. TH3F *histoDenW = new TH3F(baseName.Data(), baseTitle.Data(),
  111. nBins, qLo, qHi,
  112. nBins, qLo, qHi,
  113. nBins, qLo, qHi);
  114. histoDenW->Sumw2();
  115. mDenominatorWeighted.push_back( histoDenW );
  116. // Qinv weighted denominator
  117. baseName = title;
  118. baseName += "_den_qinv_";
  119. baseName += iKtBin;
  120. baseTitle = baseName;
  121. baseTitle += " ";
  122. baseTitle += Form("%3.2f",( mKtRange[0] + mKtStep * iKtBin) );
  123. baseTitle += "#leq kT (GeV/c) #leq";
  124. baseTitle += Form("%3.2f",( mKtRange[0] + mKtStep * (iKtBin+1) ) );
  125. baseTitle += appendix;
  126. TH3F *histoDenQinvW = new TH3F(baseName.Data(), baseTitle.Data(),
  127. nBins, qLo, qHi,
  128. nBins, qLo, qHi,
  129. nBins, qLo, qHi);
  130. histoDenQinvW->Sumw2();
  131. mDenominatorQinvWeighted.push_back( histoDenQinvW );
  132. } // if ( mIsUseDenominator )
  133. } // for (int iKtBin=0; iKtBin<mNKtBins; iKtBin++ )
  134. }
  135. //_________________
  136. MpdFemtoModelBPLCMS3DCorrFctnKt::MpdFemtoModelBPLCMS3DCorrFctnKt(const MpdFemtoModelBPLCMS3DCorrFctnKt& copy) :
  137. MpdFemtoBaseCorrFctn(copy),
  138. mManager(copy.mManager),
  139. mNKtBins(copy.mNKtBins),
  140. mKtStep(copy.mKtStep),
  141. mIsUseDenominator(copy.mIsUseDenominator) {
  142. // Copy constructor
  143. for (int i = 0; i < 2; i++) {
  144. mKtRange[i] = copy.mKtRange[i];
  145. }
  146. // Loop over kT bins and make histogram copies
  147. for ( int iCF=0; iCF<mNKtBins; iCF++ ) {
  148. mNumerator.push_back( ( copy.mNumerator.at(iCF) ) ?
  149. new TH3F( *copy.mNumerator.at(iCF) ) : nullptr );
  150. mNumeratorWeighted.push_back( ( copy.mNumeratorWeighted.at(iCF) ) ?
  151. new TH3F( *copy.mNumeratorWeighted.at(iCF) ) : nullptr );
  152. mNumeratorQinvWeighted.push_back( ( copy.mNumeratorQinvWeighted.at(iCF) ) ?
  153. new TH3F( *copy.mNumeratorQinvWeighted.at(iCF) ) : nullptr );
  154. if ( mIsUseDenominator ) {
  155. mDenominator.push_back( ( copy.mDenominator.at(iCF) ) ?
  156. new TH3F( *copy.mDenominator.at(iCF) ) : nullptr );
  157. mDenominatorWeighted.push_back( ( copy.mDenominatorWeighted.at(iCF) ) ?
  158. new TH3F( *copy.mDenominatorWeighted.at(iCF) ) : nullptr );
  159. mDenominatorQinvWeighted.push_back( ( copy.mDenominatorQinvWeighted.at(iCF) ) ?
  160. new TH3F( *copy.mDenominatorQinvWeighted.at(iCF) ) : nullptr );
  161. } // if ( mIsUseDenominator )
  162. } // for ( int iCF=0; iCF<mNKtBins; iCF++ )
  163. }
  164. //____________________________
  165. MpdFemtoModelBPLCMS3DCorrFctnKt::~MpdFemtoModelBPLCMS3DCorrFctnKt() {
  166. // Destructor
  167. if ( !mNumerator.empty() ) {
  168. mNumerator.clear();
  169. }
  170. if ( !mNumeratorWeighted.empty() ) {
  171. mNumeratorWeighted.clear();
  172. }
  173. if ( !mNumeratorQinvWeighted.empty() ) {
  174. mNumeratorQinvWeighted.clear();
  175. }
  176. if ( !mDenominator.empty() ) {
  177. mDenominator.clear();
  178. }
  179. if ( !mDenominatorWeighted.empty() ) {
  180. mDenominatorWeighted.clear();
  181. }
  182. if ( !mDenominatorQinvWeighted.empty() ) {
  183. mDenominatorQinvWeighted.clear();
  184. }
  185. }
  186. //_________________________
  187. MpdFemtoModelBPLCMS3DCorrFctnKt& MpdFemtoModelBPLCMS3DCorrFctnKt::operator=(const MpdFemtoModelBPLCMS3DCorrFctnKt& copy) {
  188. // Assignment operator
  189. if (this != &copy) {
  190. mManager = copy.mManager;
  191. mNKtBins = copy.mNKtBins;
  192. mKtStep = copy.mKtStep;
  193. mKtRange[0] = copy.mKtRange[0];
  194. mKtRange[1] = copy.mKtRange[1];
  195. mIsUseDenominator = copy.mIsUseDenominator;
  196. // Loop over kT bins and make histogram copies
  197. for ( int iCF=0; iCF<mNKtBins; iCF++ ) {
  198. mNumerator.push_back( ( copy.mNumerator.at(iCF) ) ?
  199. new TH3F( *copy.mNumerator.at(iCF) ) : nullptr );
  200. mNumeratorWeighted.push_back( ( copy.mNumeratorWeighted.at(iCF) ) ?
  201. new TH3F( *copy.mNumeratorWeighted.at(iCF) ) : nullptr );
  202. mNumeratorQinvWeighted.push_back( ( copy.mNumeratorQinvWeighted.at(iCF) ) ?
  203. new TH3F( *copy.mNumeratorQinvWeighted.at(iCF) ) : nullptr );
  204. if ( mIsUseDenominator ) {
  205. mDenominator.push_back( ( copy.mDenominator.at(iCF) ) ?
  206. new TH3F( *copy.mDenominator.at(iCF) ) : nullptr );
  207. mDenominatorWeighted.push_back( ( copy.mDenominatorWeighted.at(iCF) ) ?
  208. new TH3F( *copy.mDenominatorWeighted.at(iCF) ) : nullptr );
  209. mDenominatorQinvWeighted.push_back( ( copy.mDenominatorQinvWeighted.at(iCF) ) ?
  210. new TH3F( *copy.mDenominatorQinvWeighted.at(iCF) ) : nullptr );
  211. } // if ( mIsUseDenominator )
  212. } // for ( int iCF=0; iCF<mNKtBins; iCF++ )
  213. } // if (this != &copy)
  214. return *this;
  215. }
  216. //_________________
  217. void MpdFemtoModelBPLCMS3DCorrFctnKt::connectToManager(MpdFemtoModelManager *manager) {
  218. mManager = manager;
  219. }
  220. //_________________
  221. void MpdFemtoModelBPLCMS3DCorrFctnKt::eventBegin(const MpdFemtoEvent* /* event */) {
  222. /* empty */
  223. }
  224. //_________________
  225. void MpdFemtoModelBPLCMS3DCorrFctnKt::eventEnd(const MpdFemtoEvent* /* event */) {
  226. /* empty */
  227. }
  228. //_________________________
  229. void MpdFemtoModelBPLCMS3DCorrFctnKt::writeOutHistos() {
  230. // Write out all histograms to file
  231. if ( !mNumerator.empty() ) {
  232. for (unsigned int i=0; i<mNumerator.size(); i++ ) {
  233. mNumerator.at(i)->Write();
  234. }
  235. }
  236. if ( !mNumeratorWeighted.empty() ) {
  237. for (unsigned int i=0; i<mNumeratorWeighted.size(); i++ ) {
  238. mNumeratorWeighted.at(i)->Write();
  239. }
  240. }
  241. if ( !mNumeratorQinvWeighted.empty() ) {
  242. for (unsigned int i=0; i<mNumeratorQinvWeighted.size(); i++ ) {
  243. mNumeratorQinvWeighted.at(i)->Write();
  244. }
  245. }
  246. if ( !mDenominator.empty() ) {
  247. for (unsigned int i=0; i<mDenominator.size(); i++ ) {
  248. mDenominator.at(i)->Write();
  249. }
  250. }
  251. if ( !mDenominatorWeighted.empty() ) {
  252. for (unsigned int i=0; i<mDenominatorWeighted.size(); i++ ) {
  253. mDenominatorWeighted.at(i)->Write();
  254. }
  255. }
  256. if ( !mDenominatorQinvWeighted.empty() ) {
  257. for (unsigned int i=0; i<mDenominatorQinvWeighted.size(); i++ ) {
  258. mDenominatorQinvWeighted.at(i)->Write();
  259. }
  260. }
  261. }
  262. //______________________________
  263. TList* MpdFemtoModelBPLCMS3DCorrFctnKt::getOutputList() {
  264. // Prepare the list of objects to be written to the output
  265. TList *tOutputList = new TList();
  266. for (int i = 0; i < mNKtBins; i++) {
  267. tOutputList->Add(mNumerator.at(i));
  268. tOutputList->Add(mNumeratorWeighted.at(i));
  269. tOutputList->Add(mNumeratorQinvWeighted.at(i));
  270. if ( mIsUseDenominator ) {
  271. tOutputList->Add(mDenominator.at(i));
  272. tOutputList->Add(mDenominatorWeighted.at(i));
  273. tOutputList->Add(mDenominatorQinvWeighted.at(i));
  274. }
  275. } // for (int i=0; i<mNKtBins; i++)
  276. return tOutputList;
  277. }
  278. //_________________________
  279. void MpdFemtoModelBPLCMS3DCorrFctnKt::finish() {
  280. /* empty */
  281. }
  282. //____________________________
  283. MpdFemtoString MpdFemtoModelBPLCMS3DCorrFctnKt::report() {
  284. // Construct the report
  285. TString report("LCMS frame Bertsch-Pratt 3D correlation functions with kT binning report with model weights:\n");
  286. for (int iKt = 0; iKt < mNKtBins; iKt++) {
  287. report += TString::Format("Number of entries in %d-th numerator :\t%E\n",
  288. iKt, mNumerator.at(iKt)->GetEntries());
  289. report += TString::Format("Number of entries in %d-th weighted numerator :\t%E\n",
  290. iKt, mNumeratorWeighted.at(iKt)->GetEntries());
  291. report += TString::Format("Number of entries in %d-th qInv weighted numerator :\t%E\n",
  292. iKt, mNumeratorWeighted.at(iKt)->GetEntries());
  293. if ( mIsUseDenominator ) {
  294. report += TString::Format("Number of entries in %d-th denominator :\t%E\n",
  295. iKt, mDenominator.at(iKt)->GetEntries());
  296. report += TString::Format("Number of entries in %d-th weighted denominator :\t%E\n",
  297. iKt, mDenominatorWeighted.at(iKt)->GetEntries());
  298. report += TString::Format("Number of entries in %d-th qInv weighted denominator:\t%E\n",
  299. iKt, mDenominatorWeighted.at(iKt)->GetEntries());
  300. }
  301. } // for (int iKt=0; iKt<mNKtBins; iKt++)
  302. if (mPairCut) {
  303. report += "Report from the front-loaded PairCut specific to this CorrFctn\n";
  304. report += mPairCut->report();
  305. } else {
  306. report += "No front-loaded PairCut specific to this CorrFctn\n";
  307. }
  308. return MpdFemtoString((const char *) report);
  309. }
  310. //_________________
  311. void MpdFemtoModelBPLCMS3DCorrFctnKt::addRealPair(MpdFemtoPair* pair) {
  312. // Perform operations on real pairs
  313. // Check front-loaded pair cut
  314. if (mPairCut && !mPairCut->pass(pair)) {
  315. return;
  316. }
  317. // Find index for the current kT value
  318. int mIndexKt = (int) ( ( pair->kT() - mKtRange[0] ) / mKtStep );
  319. // Check that index is in the requested kT range
  320. if ( ( mIndexKt>=0 ) && ( mIndexKt<mNKtBins ) ) {
  321. // Check that histrogram exists
  322. if ( mNumerator.at(mIndexKt) ) {
  323. mNumerator.at(mIndexKt)->Fill(pair->qOutCMS(),
  324. pair->qSideCMS(),
  325. pair->qLongCMS(),
  326. 1.);
  327. } // if ( mNumerator[mIndexKt] )
  328. // Check that histrogram exists
  329. if ( mNumeratorWeighted.at(mIndexKt) ) {
  330. double weight = mManager->weight(pair);
  331. mNumeratorWeighted.at(mIndexKt)->Fill(pair->qOutCMS(),
  332. pair->qSideCMS(),
  333. pair->qLongCMS(),
  334. weight);
  335. } // if ( mNumeratorWeighted[mIndexKt] )
  336. // Check that histogram exists
  337. if ( mNumeratorQinvWeighted.at(mIndexKt) ) {
  338. mNumeratorQinvWeighted.at(mIndexKt)->Fill(pair->qOutCMS(),
  339. pair->qSideCMS(),
  340. pair->qLongCMS(),
  341. pair->qInv());
  342. }
  343. } // if ( (mIndexKt>=0) && (mIndexKt<mNKtBins) )
  344. }
  345. //_________________
  346. void MpdFemtoModelBPLCMS3DCorrFctnKt::addMixedPair(MpdFemtoPair* pair) {
  347. if ( !mIsUseDenominator ) return;
  348. // Check front-loaded pair cut
  349. if (mPairCut && !mPairCut->pass(pair)) {
  350. return;
  351. }
  352. // Find index for the current kT value
  353. int mIndexKt = (int) ( (pair->kT() - mKtRange[0]) / mKtStep );
  354. // Check that index is in the requested kT range
  355. if ( ( mIndexKt>=0 ) && ( mIndexKt<mNKtBins ) ) {
  356. // Check that histrogram exists
  357. if ( mDenominator.at(mIndexKt) ) {
  358. mDenominator.at(mIndexKt)->Fill(pair->qOutCMS(),
  359. pair->qSideCMS(),
  360. pair->qLongCMS(),
  361. 1.);
  362. } // if ( mDenominator[mIndexKt] )
  363. // Check that histrogram exists
  364. if ( mDenominatorWeighted.at(mIndexKt) ) {
  365. double weight = mManager->weight(pair);
  366. mDenominatorWeighted.at(mIndexKt)->Fill(pair->qOutCMS(),
  367. pair->qSideCMS(),
  368. pair->qLongCMS(),
  369. weight);
  370. } // if ( mDenominatorWeighted[mIndexKt] )
  371. // Check that histrogram exists
  372. if ( mDenominatorQinvWeighted.at(mIndexKt) ) {
  373. mDenominatorQinvWeighted.at(mIndexKt)->Fill(pair->qOutCMS(),
  374. pair->qSideCMS(),
  375. pair->qLongCMS(),
  376. pair->qInv());
  377. } // if ( mDenominatorQinvWeighted[mIndexKt] )
  378. } // if ( (mIndexKt>=0) && (mIndexKt<mNKtBins) )
  379. }
  380. //_________________
  381. void MpdFemtoModelBPLCMS3DCorrFctnKt::setKtRange(const int& nKtBins, const float& kTLow, const float& kTHi) {
  382. if (nKtBins >= 1) {
  383. mNKtBins = nKtBins;
  384. } else {
  385. std::cout << "[WARNING} void MpdFemtoModelBPLCMS3DCorrFctnKt::setKtRange - "
  386. << "Number of kT bins must be positive" << std::endl;
  387. mNKtBins = 1;
  388. }
  389. if (kTLow < kTHi) {
  390. mKtRange[0] = kTLow;
  391. mKtRange[1] = kTHi;
  392. } else {
  393. std::cout << "[WARNING} void MpdFemtoModelBPLCMS3DCorrFctnKt::setKtRange - "
  394. << "kTLow < kTHi. Resetting to the integral" << std::endl;
  395. mKtRange[0] = 0.15;
  396. mKtRange[1] = 1.05;
  397. }
  398. mKtStep = (mKtRange[1] - mKtRange[0]) / mNKtBins;
  399. }