MpdFemtoModelQinvCorrFctnKt.cxx 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359
  1. //
  2. // One-dimensional correlation function with kT-binning
  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 "MpdFemtoModelQinvCorrFctnKt.h"
  11. // ROOT headers
  12. #include "TString.h"
  13. #include "TH1D.h"
  14. // C++ headers
  15. #include <iostream>
  16. ClassImp(MpdFemtoModelQinvCorrFctnKt)
  17. //_________________
  18. MpdFemtoModelQinvCorrFctnKt::MpdFemtoModelQinvCorrFctnKt(const char* title,
  19. const int &nQbins, const double &qLo, const double &qHi,
  20. const int& nKtBins, const double& ktLo, const double& ktHi,
  21. const bool useDenominator) :
  22. MpdFemtoBaseCorrFctn(),
  23. mManager(nullptr) {
  24. // Constructor
  25. // Set general parameters
  26. setKtRange(nKtBins, ktLo, ktHi);
  27. setUseDenominator( useDenominator );
  28. TString baseName(title);
  29. TString baseTitle(title);
  30. TString appendix = ";q_{inv} (GeV/c);Entries";
  31. // Histogram loop
  32. for (int iKtBin = 0; iKtBin < mNKtBins; iKtBin++) {
  33. // Numerator
  34. baseName = title;
  35. baseName += "_num_";
  36. baseName += iKtBin;
  37. baseTitle = baseName;
  38. baseTitle += " ";
  39. baseTitle += Form("%3.2f", (mKtRange[0] + mKtStep * iKtBin));
  40. baseTitle += "#leq k_{T} (GeV/c) #leq";
  41. baseTitle += Form("%3.2f", (mKtRange[0] + mKtStep * (iKtBin + 1)));
  42. baseTitle += appendix;
  43. TH1D *histoNum = new TH1D(baseName.Data(), baseTitle.Data(),
  44. nQbins, qLo, qHi);
  45. histoNum->Sumw2();
  46. mNumerator.push_back(histoNum);
  47. // Numerator weighted
  48. baseName = title;
  49. baseName += "_num_wei_";
  50. baseName += iKtBin;
  51. baseTitle = baseName;
  52. baseTitle += " ";
  53. baseTitle += Form("%3.2f", (mKtRange[0] + mKtStep * iKtBin));
  54. baseTitle += "#leq k_{T} (GeV/c) #leq";
  55. baseTitle += Form("%3.2f", (mKtRange[0] + mKtStep * (iKtBin + 1)));
  56. baseTitle += appendix;
  57. TH1D *histoNumW = new TH1D(baseName.Data(), baseTitle.Data(),
  58. nQbins, qLo, qHi);
  59. histoNumW->Sumw2();
  60. mNumeratorWeighted.push_back(histoNumW);
  61. if ( mIsUseDenominator ) {
  62. // Denominator
  63. baseName = title;
  64. baseName += "_den_";
  65. baseName += iKtBin;
  66. baseTitle = baseName;
  67. baseTitle += " ";
  68. baseTitle += Form("%3.2f", (mKtRange[0] + mKtStep * iKtBin));
  69. baseTitle += "#leq k_{T} (GeV/c) #leq";
  70. baseTitle += Form("%3.2f", (mKtRange[0] + mKtStep * (iKtBin + 1)));
  71. baseTitle += appendix;
  72. TH1D *histoDen = new TH1D(baseName.Data(), baseTitle.Data(),
  73. nQbins, qLo, qHi);
  74. histoDen->Sumw2();
  75. mDenominator.push_back(histoDen);
  76. // Denominator weighted
  77. baseName = title;
  78. baseName += "_den_wei_";
  79. baseName += iKtBin;
  80. baseTitle = baseName;
  81. baseTitle += " ";
  82. baseTitle += Form("%3.2f", (mKtRange[0] + mKtStep * iKtBin));
  83. baseTitle += "#leq k_{T} (GeV/c) #leq";
  84. baseTitle += Form("%3.2f", (mKtRange[0] + mKtStep * (iKtBin + 1)));
  85. baseTitle += appendix;
  86. TH1D *histoDenW = new TH1D(baseName.Data(), baseTitle.Data(),
  87. nQbins, qLo, qHi);
  88. histoDenW->Sumw2();
  89. mDenominatorWeighted.push_back(histoDen);
  90. } // if ( mIsUseDenominator )
  91. } // for(int i=0; i<mNumberKt; i++)
  92. }
  93. //_________________
  94. MpdFemtoModelQinvCorrFctnKt::MpdFemtoModelQinvCorrFctnKt(const MpdFemtoModelQinvCorrFctnKt& c) :
  95. MpdFemtoBaseCorrFctn(c),
  96. mManager(c.mManager),
  97. mNKtBins(c.mNKtBins),
  98. mKtStep(c.mKtStep),
  99. mIsUseDenominator(c.mIsUseDenominator) {
  100. // Copy constructor
  101. mKtRange[0] = c.mKtRange[0];
  102. mKtRange[1] = c.mKtRange[1];
  103. for (int iCF = 0; iCF < mNKtBins; iCF++) {
  104. mNumerator.push_back( (c.mNumerator.at(iCF)) ?
  105. new TH1D(*c.mNumerator.at(iCF)) : nullptr);
  106. mNumeratorWeighted.push_back( (c.mNumeratorWeighted.at(iCF)) ?
  107. new TH1D(*c.mNumeratorWeighted.at(iCF)) : nullptr);
  108. if ( mIsUseDenominator ) {
  109. mDenominator.push_back( (c.mDenominator.at(iCF)) ?
  110. new TH1D(*c.mDenominator.at(iCF)) : nullptr);
  111. mDenominatorWeighted.push_back( (c.mDenominatorWeighted.at(iCF)) ?
  112. new TH1D(*c.mDenominatorWeighted.at(iCF)) : nullptr);
  113. }
  114. } // for ( int iBin=0; iBin<mNumberKt; iBin++ )
  115. }
  116. //_________________
  117. MpdFemtoModelQinvCorrFctnKt& MpdFemtoModelQinvCorrFctnKt::operator=(const MpdFemtoModelQinvCorrFctnKt& c) {
  118. // Assignment operator
  119. if (this != &c) {
  120. mManager = c.mManager;
  121. mNKtBins = c.mNKtBins;
  122. mKtRange[0] = c.mKtRange[0];
  123. mKtRange[1] = c.mKtRange[1];
  124. mKtStep = c.mKtStep;
  125. mIsUseDenominator = c.mIsUseDenominator;
  126. // Loop over kT bins
  127. for (int iCF = 0; iCF < mNKtBins; iCF++) {
  128. mNumerator.push_back( (c.mNumerator.at(iCF)) ?
  129. new TH1D(*c.mNumerator.at(iCF)) : nullptr);
  130. mNumeratorWeighted.push_back( (c.mNumeratorWeighted.at(iCF)) ?
  131. new TH1D(*c.mNumeratorWeighted.at(iCF)) : nullptr);
  132. if ( mIsUseDenominator ) {
  133. mDenominator.push_back( (c.mDenominator.at(iCF)) ?
  134. new TH1D(*c.mDenominator.at(iCF)) : nullptr);
  135. mDenominatorWeighted.push_back( (c.mDenominatorWeighted.at(iCF)) ?
  136. new TH1D(*c.mDenominatorWeighted.at(iCF)) : nullptr);
  137. }
  138. } // for ( int iBin=0; iBin<mNumberKt; iBin++ )
  139. } // if (this != &c)
  140. return *this;
  141. }
  142. //_________________
  143. MpdFemtoModelQinvCorrFctnKt::~MpdFemtoModelQinvCorrFctnKt() {
  144. // Destructor
  145. if (!mNumerator.empty()) {
  146. mNumerator.clear();
  147. }
  148. if (!mNumeratorWeighted.empty()) {
  149. mNumeratorWeighted.clear();
  150. }
  151. if (!mDenominator.empty()) {
  152. mDenominator.clear();
  153. }
  154. if (!mDenominatorWeighted.empty()) {
  155. mDenominatorWeighted.clear();
  156. }
  157. delete mManager;
  158. }
  159. //_________________
  160. void MpdFemtoModelQinvCorrFctnKt::writeOutHistos() {
  161. // Write histograms to the file
  162. if (!mNumerator.empty()) {
  163. for (unsigned int i = 0; i < mNumerator.size(); i++) {
  164. mNumerator.at(i)->Write();
  165. }
  166. }
  167. if (!mNumeratorWeighted.empty()) {
  168. for (unsigned int i = 0; i < mNumeratorWeighted.size(); i++) {
  169. mNumeratorWeighted.at(i)->Write();
  170. }
  171. }
  172. if ( mIsUseDenominator ) {
  173. if (!mDenominator.empty()) {
  174. for (unsigned int i = 0; i < mDenominator.size(); i++) {
  175. mDenominator.at(i)->Write();
  176. }
  177. }
  178. if (!mDenominatorWeighted.empty()) {
  179. for (unsigned int i = 0; i < mDenominatorWeighted.size(); i++) {
  180. mDenominatorWeighted.at(i)->Write();
  181. }
  182. }
  183. }
  184. }
  185. //_________________
  186. TList* MpdFemtoModelQinvCorrFctnKt::getOutputList() {
  187. // Prepare the list of objects to be written to the output
  188. TList *outputList = new TList();
  189. for (int i = 0; i < mNKtBins; i++) {
  190. outputList->Add(mNumerator.at(i));
  191. outputList->Add(mNumeratorWeighted.at(i));
  192. if ( mIsUseDenominator ) {
  193. outputList->Add(mDenominator.at(i));
  194. outputList->Add(mDenominatorWeighted.at(i));
  195. }
  196. }
  197. return outputList;
  198. }
  199. //_________________
  200. void MpdFemtoModelQinvCorrFctnKt::eventBegin(const MpdFemtoEvent* /* event */) {
  201. /* empty */
  202. }
  203. //_________________
  204. void MpdFemtoModelQinvCorrFctnKt::eventEnd(const MpdFemtoEvent* /* event */) {
  205. /* empty */
  206. }
  207. //_________________
  208. void MpdFemtoModelQinvCorrFctnKt::finish() {
  209. /* empty */
  210. }
  211. //_________________
  212. void MpdFemtoModelQinvCorrFctnKt::setKtRange(const int& nbins, const double& ktLo, const double& ktHi) {
  213. mNKtBins = nbins;
  214. mKtRange[0] = ktLo;
  215. mKtRange[1] = ktHi;
  216. mKtStep = (mKtRange[1] - mKtRange[0]) / mNKtBins;
  217. }
  218. //_________________
  219. void MpdFemtoModelQinvCorrFctnKt::addRealPair(MpdFemtoPair* pair) {
  220. // Check if pair passes front-loaded cut if exists
  221. if (mPairCut && !mPairCut->pass(pair)) {
  222. return;
  223. }
  224. int mIndexKt = (int) ( (pair->kT() - mKtRange[0]) / mKtStep );
  225. if ((mIndexKt >= 0) && (mIndexKt < mNKtBins)) {
  226. if ( mNumerator.at(mIndexKt) ) {
  227. mNumerator.at(mIndexKt)->Fill( fabs(pair->qInv()), 1.);
  228. } else {
  229. std::cout << "[ERROR] void MpdFemtoModelQinvCorrFctnKt::addRealPair(MpdFemtoPair* pair) - "
  230. << "No correlation function with index: " << mIndexKt << " was found"
  231. << std::endl;
  232. }
  233. if ( mNumeratorWeighted.at(mIndexKt) ) {
  234. double weight = mManager->weight(pair);
  235. mNumeratorWeighted.at(mIndexKt)->Fill( fabs(pair->qInv()), weight);
  236. } else {
  237. std::cout << "[ERROR] void MpdFemtoModelQinvCorrFctnKt::addRealPair(MpdFemtoPair* pair) - "
  238. << "No weighted correlation function with index: " << mIndexKt << " was found"
  239. << std::endl;
  240. }
  241. } // if ( ( mIndexKt>=0 ) && ( mIndexKt<mNumberKt ) )
  242. }
  243. //_________________
  244. void MpdFemtoModelQinvCorrFctnKt::addMixedPair(MpdFemtoPair* pair) {
  245. if ( !mIsUseDenominator ) return;
  246. // Check if pair passes front-loaded cut if exists
  247. if (mPairCut && !mPairCut->pass(pair)) {
  248. return;
  249. }
  250. int mIndexKt = (int) ( (fabs(pair->kT()) - mKtRange[0]) / mKtStep );
  251. if ( (mIndexKt >= 0) && (mIndexKt < mNKtBins) ) {
  252. if (mDenominator.at(mIndexKt)) {
  253. mDenominator.at(mIndexKt)->Fill( fabs(pair->qInv()) );
  254. } else {
  255. std::cout << "[ERROR] void MpdFemtoModelQinvCorrFctnKt::addMixedPair(MpdFemtoPair* pair) - "
  256. << "No correlation function with index: " << mIndexKt << " was found"
  257. << std::endl;
  258. }
  259. if (mDenominatorWeighted.at(mIndexKt)) {
  260. double weight = mManager->weight(pair);
  261. mDenominatorWeighted.at(mIndexKt)->Fill( fabs(pair->qInv()), weight);
  262. } else {
  263. std::cout << "[ERROR] void MpdFemtoModelQinvCorrFctnKt::addMixedPair(MpdFemtoPair* pair) - "
  264. << "No weighted correlation function with index: " << mIndexKt << " was found"
  265. << std::endl;
  266. }
  267. } // if ( ( mIndexKt>=0 ) && ( mIndexKt<mNumberKt ) )
  268. }
  269. //_________________
  270. MpdFemtoString MpdFemtoModelQinvCorrFctnKt::report() {
  271. // Make a report
  272. MpdFemtoString tStr = "MpdFemtoModelQinvCorrFctnKt report";
  273. int mNumeratorEntries = 0;
  274. int mNumeratorWEntries = 0;
  275. int mDenominatorEntries = 0;
  276. int mDenominatorWEntries = 0;
  277. for (int i = 0; i < mNKtBins; i++) {
  278. mNumeratorEntries += (int) mNumerator.at(i)->GetEntries();
  279. mNumeratorWEntries += (int) mNumeratorWeighted.at(i)->GetEntries();
  280. if ( mIsUseDenominator ) {
  281. mDenominatorEntries += (int) mDenominator.at(i)->GetEntries();
  282. mDenominatorWEntries += (int) mDenominatorWeighted.at(i)->GetEntries();
  283. }
  284. } // for (int i=0; i<mNumberKt; i++)
  285. tStr += TString::Format("\nTotal number of pairs in the numerator :\t%d\n", mNumeratorEntries);
  286. tStr += TString::Format("\nTotal number of pairs in the weighted numerator:\t%d\n", mNumeratorEntries);
  287. if ( mIsUseDenominator ) {
  288. tStr += TString::Format("Total number of pairs in the denominator :\t%d\n", mDenominatorEntries);
  289. tStr += TString::Format("Total number of pairs in the weighted denominator:\t%d\n", mDenominatorEntries);
  290. }
  291. for (int i = 0; i < mNKtBins; i++) {
  292. tStr += TString::Format("Total number of pairs in %d-th numerator :\t%E\n",
  293. i, mNumerator.at(i)->GetEntries());
  294. tStr += TString::Format("Total number of pairs in %d-th weighted numerator:\t%E\n",
  295. i, mNumeratorWeighted.at(i)->GetEntries());
  296. if ( mIsUseDenominator ) {
  297. tStr += TString::Format("Total number of pairs in %d-th denominator :\t%E\n",
  298. i, mDenominator.at(i)->GetEntries());
  299. tStr += TString::Format("Total number of pairs in %d-th weighted denominator:\t%E\n",
  300. i, mDenominatorWeighted.at(i)->GetEntries());
  301. }
  302. }
  303. if (mPairCut) {
  304. tStr += "Here is the PairCut specific to this CorrFctn\n";
  305. tStr += mPairCut->report();
  306. } else {
  307. tStr += "No PairCut specific to this CorrFctn\n";
  308. }
  309. return tStr;
  310. }
  311. //_________________
  312. void MpdFemtoModelQinvCorrFctnKt::connectToManager(MpdFemtoModelManager *manager) {
  313. mManager = manager;
  314. }