MpdFemtoDeltaEtaDeltaPhiStarMinKt.cxx 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349
  1. //
  2. // Two-dimensional $\Delta\eta$ vs. $\Delta\phi^{*}_{min}$
  3. //
  4. // MpdFemtoMaker headers
  5. #include "MpdFemtoPair.h"
  6. // MpdFemtoMakerUser headers
  7. #include "MpdFemtoDeltaEtaDeltaPhiStarMinKt.h"
  8. // ROOT headers
  9. #include "TMath.h"
  10. #include "TString.h"
  11. ClassImp(MpdFemtoDeltaEtaDeltaPhiStarMinKt)
  12. //_________________
  13. MpdFemtoDeltaEtaDeltaPhiStarMinKt::MpdFemtoDeltaEtaDeltaPhiStarMinKt(const char* title,
  14. const int& nBinsEta, const double& etaLo,
  15. const double& etaHi,
  16. const int& nBinsPhi, const double& phiLo,
  17. const double& phiHi,
  18. const int& ktBins, const double& ktLo,
  19. const double& ktHi) :
  20. MpdFemtoBaseCorrFctn() {
  21. // Parametrized constructor
  22. // Set initial parameters
  23. setCharges(1, 1);
  24. setBField(-0.5);
  25. mRadiusStep = 0.35;
  26. mRadii[0] = 0.6;
  27. mRadii[1] = 2.;
  28. setHistoParameters(nBinsEta, etaLo, etaHi, nBinsPhi, phiLo, phiHi);
  29. setKtRange(ktBins, ktLo, ktHi);
  30. // Define string parameters
  31. TString baseName(title);
  32. TString baseTitle(title);
  33. TString appendix = ";#Delta#phi^{*}_{min} (rad);#Delta#eta";
  34. // Loop over kTBins
  35. for (int iKtBin = 0; iKtBin < mNKtBins; iKtBin++) {
  36. // Set numerator
  37. baseName = title;
  38. baseName += "_num_";
  39. baseName += iKtBin;
  40. baseTitle = baseName;
  41. baseTitle += " ";
  42. baseTitle += Form("%3.2f", (mKtRange[0] + mKtStep * iKtBin));
  43. baseTitle += "#leq k_{T} (GeV/c) #leq ";
  44. baseTitle += Form("%3.2f", (mKtRange[0] + mKtStep * (iKtBin + 1)));
  45. baseTitle += appendix;
  46. TH2F *hNumer = new TH2F(baseName.Data(), baseTitle.Data(),
  47. mPhiBins, mPhiRange[0], mPhiRange[1],
  48. mEtaBins, mEtaRange[0], mEtaRange[1]);
  49. hNumer->Sumw2();
  50. mNumerator.push_back(hNumer);
  51. // Denominator
  52. baseName = title;
  53. baseName += "_den_";
  54. baseName += iKtBin;
  55. baseTitle = baseName;
  56. baseTitle += " ";
  57. baseTitle += Form("%3.2f", (mKtRange[0] + mKtStep * iKtBin));
  58. baseTitle += "#leq k_{T} (GeV/c) #leq ";
  59. baseTitle += Form("%3.2f", (mKtRange[0] + mKtStep * (iKtBin + 1)));
  60. baseTitle += appendix;
  61. TH2F *hDenom = new TH2F(baseName.Data(), baseTitle.Data(),
  62. mPhiBins, mPhiRange[0], mPhiRange[1],
  63. mEtaBins, mEtaRange[0], mEtaRange[1]);
  64. hDenom->Sumw2();
  65. mDenominator.push_back(hDenom);
  66. } // for (int iKtBin=0; iKtBin<mNKtBins; iKtBin++ )
  67. }
  68. //_________________
  69. MpdFemtoDeltaEtaDeltaPhiStarMinKt::MpdFemtoDeltaEtaDeltaPhiStarMinKt(const MpdFemtoDeltaEtaDeltaPhiStarMinKt& copy) :
  70. MpdFemtoBaseCorrFctn(copy) {
  71. // Copy constructor
  72. setHistoParameters(copy.mEtaBins, copy.mEtaRange[0], copy.mEtaRange[1],
  73. copy.mPhiBins, copy.mPhiRange[0], copy.mPhiRange[1]);
  74. setKtRange(copy.mNKtBins, copy.mKtRange[0], copy.mKtRange[1]);
  75. setCharges(copy.mQ1, copy.mQ2);
  76. setBField(copy.mBField);
  77. mRadiusStep = copy.mRadiusStep;
  78. mRadii[0] = copy.mRadii[0];
  79. mRadii[1] = copy.mRadii[1];
  80. // Loop over kT bins and make histogram copies
  81. for (int iKtBin = 0; iKtBin < mNKtBins; iKtBin++) {
  82. mNumerator.push_back(copy.mNumerator.at(iKtBin));
  83. mNumerator.back()->Sumw2();
  84. mDenominator.push_back(copy.mDenominator.at(iKtBin));
  85. mDenominator.back()->Sumw2();
  86. } // for (int iKtBin=0; iKtBin<mNKtBins; iKtBin++)
  87. }
  88. //____________________________
  89. MpdFemtoDeltaEtaDeltaPhiStarMinKt::~MpdFemtoDeltaEtaDeltaPhiStarMinKt() {
  90. // Destructor
  91. if (!mNumerator.empty()) {
  92. mNumerator.clear();
  93. }
  94. if (!mDenominator.empty()) {
  95. mDenominator.clear();
  96. }
  97. }
  98. //_________________________
  99. MpdFemtoDeltaEtaDeltaPhiStarMinKt& MpdFemtoDeltaEtaDeltaPhiStarMinKt::operator=(const MpdFemtoDeltaEtaDeltaPhiStarMinKt& copy) {
  100. // Assignment operator
  101. if (this != &copy) {
  102. mPhiBins = copy.mPhiBins;
  103. mEtaBins = copy.mEtaBins;
  104. mNKtBins = copy.mNKtBins;
  105. mKtStep = copy.mKtStep;
  106. for (int i = 0; i < 2; i++) {
  107. mEtaRange[i] = copy.mEtaRange[i];
  108. mPhiRange[i] = copy.mPhiRange[i];
  109. mKtRange[i] = copy.mKtRange[i];
  110. }
  111. // Loop over kT bins and make histogram copies
  112. for (int iKtBin = 0; iKtBin < mNKtBins; iKtBin++) {
  113. mNumerator.push_back(copy.mNumerator.at(iKtBin));
  114. mNumerator.back()->Sumw2();
  115. mDenominator.push_back(copy.mDenominator.at(iKtBin));
  116. mDenominator.back()->Sumw2();
  117. } // for (int iKtBin=0; iKtBin<mNKtBins; iKtBin++)
  118. } // if (this != &copy)
  119. return *this;
  120. }
  121. //_________________
  122. void MpdFemtoDeltaEtaDeltaPhiStarMinKt::eventBegin(const MpdFemtoEvent* event) {
  123. setBField(event->bField());
  124. }
  125. //_________________
  126. void MpdFemtoDeltaEtaDeltaPhiStarMinKt::eventEnd(const MpdFemtoEvent* /* event */) {
  127. /* empty */
  128. }
  129. //_________________________
  130. void MpdFemtoDeltaEtaDeltaPhiStarMinKt::writeOutHistos() {
  131. // Write out all histograms to file
  132. for (unsigned int i = 0; i < mNumerator.size(); i++) {
  133. mNumerator.at(i)->Write();
  134. }
  135. for (unsigned int i = 0; i < mDenominator.size(); i++) {
  136. mDenominator.at(i)->Write();
  137. }
  138. }
  139. //______________________________
  140. TList* MpdFemtoDeltaEtaDeltaPhiStarMinKt::getOutputList() {
  141. // Prepare the list of objects to be written to the output
  142. TList *tOutputList = new TList();
  143. for (int i = 0; i < mNKtBins; i++) {
  144. tOutputList->Add(mNumerator.at(i));
  145. tOutputList->Add(mDenominator.at(i));
  146. } // for (int i=0; i<mNKtBins; i++)
  147. return tOutputList;
  148. }
  149. //_________________________
  150. void MpdFemtoDeltaEtaDeltaPhiStarMinKt::finish() {
  151. /* empty */
  152. }
  153. //____________________________
  154. MpdFemtoString MpdFemtoDeltaEtaDeltaPhiStarMinKt::report() {
  155. // Construct the report
  156. TString report("DeltaEta vs. DeltaPhiStarMin with kT binning report:\n");
  157. for (int iKt = 0; iKt < mNKtBins; iKt++) {
  158. report += TString::Format("Number of entries in %d-th numerator :\t%E\n",
  159. iKt, mNumerator.at(iKt)->GetEntries());
  160. report += TString::Format("Number of entries in %d-th denominator :\t%E\n",
  161. iKt, mDenominator.at(iKt)->GetEntries());
  162. } // for (int iKt=0; iKt<mNKtBins; iKt++)
  163. if (mPairCut) {
  164. report += "Report from the front-loaded PairCut specific to this CorrFctn\n";
  165. report += mPairCut->report();
  166. } else {
  167. report += "No front-loaded PairCut specific to this CorrFctn\n";
  168. }
  169. return MpdFemtoString((const char *) report);
  170. }
  171. //_________________
  172. void MpdFemtoDeltaEtaDeltaPhiStarMinKt::setRadiiParameters(const int& npoints,
  173. const double& radiusLo,
  174. const double& radiusHi) {
  175. if (radiusLo < radiusHi) {
  176. mRadii[0] = radiusLo;
  177. mRadii[1] = radiusHi;
  178. } else {
  179. mRadii[0] = 0.6;
  180. mRadii[1] = 2.;
  181. }
  182. if (npoints > 0) {
  183. mRadiusStep = (mRadii[1] - mRadii[0]) / npoints;
  184. } else {
  185. mRadiusStep = (mRadii[1] - mRadii[0]) / 2;
  186. }
  187. }
  188. //_________________
  189. void MpdFemtoDeltaEtaDeltaPhiStarMinKt::setCharges(const int& charge1, const int& charge2) {
  190. mQ1 = charge1;
  191. mQ2 = charge2;
  192. }
  193. //_________________
  194. void MpdFemtoDeltaEtaDeltaPhiStarMinKt::setBField(const double& bfield) {
  195. mBField = bfield;
  196. }
  197. //_________________
  198. void MpdFemtoDeltaEtaDeltaPhiStarMinKt::addRealPair(MpdFemtoPair* pair) {
  199. // Perform operations on real pairs
  200. // Check front-loaded pair cut
  201. if (mPairCut && !mPairCut->pass(pair)) {
  202. return;
  203. }
  204. // Find index for the current kT value
  205. int mIndexKt = (int) ((pair->kT() - mKtRange[0]) / mKtStep);
  206. // Check that index is in the requested kT range
  207. if ((mIndexKt >= 0) && (mIndexKt < mNKtBins) && pair->qInv() < 0.15) {
  208. double dPhiStarMin = -1000.;
  209. mMomentum1 = pair->track1()->momentum();
  210. mMomentum2 = pair->track2()->momentum();
  211. dPhiStarMin = MpdFemtoPair::calculateDPhiStarMin(mMomentum1, mQ1,
  212. mMomentum2, mQ2,
  213. mRadiusStep, mRadii[0], mRadii[1],
  214. mBField);
  215. // Check that histrogram exists
  216. if (mNumerator.at(mIndexKt)) {
  217. mNumerator.at(mIndexKt)->Fill(dPhiStarMin, pair->deltaEta(), 1.);
  218. } // if ( mNumerator.at(mIndexKt) )
  219. } // if ( (mIndexKt>=0) && (mIndexKt<mNKtBins) )
  220. }
  221. //_________________
  222. void MpdFemtoDeltaEtaDeltaPhiStarMinKt::addMixedPair(MpdFemtoPair* pair) {
  223. // Check front-loaded pair cut
  224. if (mPairCut && !mPairCut->pass(pair)) {
  225. return;
  226. }
  227. // Find index for the current kT value
  228. int mIndexKt = (int) ((pair->kT() - mKtRange[0]) / mKtStep);
  229. // Check that index is in the requested kT range
  230. if ((mIndexKt >= 0) && (mIndexKt < mNKtBins) && pair->qInv() < 0.15) {
  231. double dPhiStarMin = -1000.;
  232. mMomentum1 = pair->track1()->momentum();
  233. mMomentum2 = pair->track2()->momentum();
  234. dPhiStarMin = MpdFemtoPair::calculateDPhiStarMin(mMomentum1, mQ1,
  235. mMomentum2, mQ2,
  236. mRadiusStep, mRadii[0], mRadii[1],
  237. mBField);
  238. // Check that histrogram exists
  239. if (mDenominator.at(mIndexKt)) {
  240. mDenominator.at(mIndexKt)->Fill(dPhiStarMin, pair->deltaEta(), 1.);
  241. } // if ( mDenominator.at(mIndexKt) )
  242. } // if ( (mIndexKt>=0) && (mIndexKt<mNKtBins) )
  243. }
  244. //_________________
  245. void MpdFemtoDeltaEtaDeltaPhiStarMinKt::setKtRange(const int& nKtBins, const double& kTLow, const double& kTHi) {
  246. if (nKtBins >= 1) {
  247. mNKtBins = nKtBins;
  248. } else {
  249. std::cout << "[WARNING} void MpdFemtoDeltaEtaDeltaPhiStarMinKt::setKtRange - "
  250. << "Number of kT bins must be positive" << std::endl;
  251. mNKtBins = 1;
  252. }
  253. if (kTLow < kTHi) {
  254. mKtRange[0] = kTLow;
  255. mKtRange[1] = kTHi;
  256. } else {
  257. std::cout << "[WARNING} void MpdFemtoDeltaEtaDeltaPhiStarMinKt::setKtRange - "
  258. << "kTLow < kTHi. Resetting to the integral" << std::endl;
  259. mKtRange[0] = 0.;
  260. mKtRange[1] = 1.5;
  261. }
  262. mKtStep = (mKtRange[1] - mKtRange[0]) / mNKtBins;
  263. }
  264. //_________________
  265. void MpdFemtoDeltaEtaDeltaPhiStarMinKt::setHistoParameters(const int& nBinsEta, const double& etaLo, const double& etaHi,
  266. const int& nBinsPhi, const double& phiLo, const double& phiHi) {
  267. if (etaLo > etaHi) {
  268. std::cout << "[WARNING] void MpdFemtoDeltaEtaDeltaPhiStarMinKt::void setHistoParameters - "
  269. << " etaLo > etaHi. Resetting to defaults." << std::endl;
  270. mEtaRange[0] = -0.5;
  271. mEtaRange[1] = 0.5;
  272. } else {
  273. mEtaRange[0] = etaLo;
  274. mEtaRange[1] = etaHi;
  275. }
  276. if (nBinsEta < 1) {
  277. std::cout << "[WARNING] void MpdFemtoDeltaEtaDeltaPhiStarMinKt::void setHistoParameters - "
  278. << " nBinsEta<1. Resetting to defaults." << std::endl;
  279. mEtaBins = 200;
  280. } else {
  281. mEtaBins = nBinsEta;
  282. }
  283. if (phiLo > phiHi) {
  284. std::cout << "[WARNING] void MpdFemtoDeltaEtaDeltaPhiStarMinKt::void setHistoParameters - "
  285. << " phiLo > phiHi. Resetting to defaults." << std::endl;
  286. mPhiRange[0] = -0.5;
  287. mPhiRange[1] = 0.5;
  288. } else {
  289. mPhiRange[0] = phiLo;
  290. mPhiRange[1] = phiHi;
  291. }
  292. if (nBinsPhi < 1) {
  293. std::cout << "[WARNING] void MpdFemtoDeltaEtaDeltaPhiStarMinKt::void setHistoParameters - "
  294. << " nBinsPhi<1. Resetting to defaults." << std::endl;
  295. mPhiBins = 200;
  296. } else {
  297. mPhiBins = nBinsPhi;
  298. }
  299. }