MpdFemtoBasicPairCut.cxx 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328
  1. //
  2. // A pair cut that keeps most of the pair observables
  3. //
  4. // MpdFemtoMaker headers
  5. #include "MpdFemtoPair.h"
  6. // MpdFemtoMakerUser headers
  7. #include "MpdFemtoBasicPairCut.h"
  8. // ROOT headers
  9. #include "TMath.h"
  10. #include "TString.h"
  11. // C++ headers
  12. #include <limits>
  13. ClassImp(MpdFemtoBasicPairCut)
  14. //_________________
  15. MpdFemtoBasicPairCut::MpdFemtoBasicPairCut() : MpdFemtoBasePairCut() {
  16. // Default constructor
  17. mPrimaryVertex.SetXYZ(0., 0., 0.);
  18. mQuality[0] = -1e6;
  19. mQuality[1] = 1e6;
  20. mKt[0] = -1e6;
  21. mKt[1] = 1e6;
  22. mPt[0] = -1e6;
  23. mPt[1] = 1e6;
  24. mOpeningAngle[0] = -1e6;
  25. mOpeningAngle[1] = 1e6;
  26. mRapidity[0] = -1e6;
  27. mRapidity[1] = 1e6;
  28. mEta[0] = -1e6;
  29. mEta[1] = 1e6;
  30. mQinv[0] = -1e6;
  31. mQinv[1] = 1e6;
  32. mMinv[0] = -1e6;
  33. mMinv[1] = 1e6;
  34. mEntranceSeparation[0] = -1e6;
  35. mEntranceSeparation[1] = 1e6;
  36. mAngleToPrimaryVertex[0] = -1e6;
  37. mAngleToPrimaryVertex[1] = 1e6;
  38. mFracOfMergedRow[0] = -1e6;
  39. mFracOfMergedRow[1] = 1e6;
  40. mClosestRowAtDCA[0] = -1e6;
  41. mClosestRowAtDCA[1] = 1e6;
  42. mWeightedAvSep[0] = -1e6;
  43. mWeightedAvSep[1] = 1e6;
  44. mAverageSeparation[0] = -1e6;
  45. mAverageSeparation[1] = 1e6;
  46. mRValueLo = -1e6;
  47. mDPhiStarMin[0] = 0.;
  48. mDPhiStarMin[1] = 0.;
  49. mDEta[0] = 0.;
  50. mDEta[1] = 0.;
  51. mNPairsPassed = 0;
  52. mNPairsFailed = 0;
  53. mVerbose = false;
  54. }
  55. //_________________
  56. MpdFemtoBasicPairCut::MpdFemtoBasicPairCut(const MpdFemtoBasicPairCut& c) : MpdFemtoBasePairCut(c) {
  57. // Copy constructor
  58. mNPairsPassed = 0;
  59. mNPairsFailed = 0;
  60. mPrimaryVertex = c.mPrimaryVertex;
  61. mQuality[0] = c.mQuality[0];
  62. mQuality[1] = c.mQuality[1];
  63. mKt[0] = c.mKt[0];
  64. mKt[1] = c.mKt[1];
  65. mPt[0] = c.mPt[0];
  66. mPt[1] = c.mPt[1];
  67. mOpeningAngle[0] = c.mOpeningAngle[0];
  68. mOpeningAngle[1] = c.mOpeningAngle[1];
  69. mQinv[0] = c.mQinv[0];
  70. mQinv[1] = c.mQinv[1];
  71. mMinv[0] = c.mMinv[0];
  72. mMinv[1] = c.mMinv[1];
  73. mRapidity[0] = c.mRapidity[0];
  74. mRapidity[1] = c.mRapidity[1];
  75. mEta[0] = c.mEta[0];
  76. mEta[1] = c.mEta[1];
  77. mEntranceSeparation[0] = c.mEntranceSeparation[0];
  78. mEntranceSeparation[1] = c.mEntranceSeparation[1];
  79. mAngleToPrimaryVertex[0] = c.mAngleToPrimaryVertex[0];
  80. mAngleToPrimaryVertex[1] = c.mAngleToPrimaryVertex[1];
  81. mFracOfMergedRow[0] = c.mFracOfMergedRow[0];
  82. mFracOfMergedRow[1] = c.mFracOfMergedRow[1];
  83. mClosestRowAtDCA[0] = c.mClosestRowAtDCA[0];
  84. mClosestRowAtDCA[1] = c.mClosestRowAtDCA[1];
  85. mWeightedAvSep[0] = c.mWeightedAvSep[0];
  86. mWeightedAvSep[1] = c.mWeightedAvSep[1];
  87. mAverageSeparation[0] = c.mAverageSeparation[0];
  88. mAverageSeparation[1] = c.mAverageSeparation[1];
  89. mRValueLo = c.mRValueLo;
  90. mDPhiStarMin[0] = c.mDPhiStarMin[0];
  91. mDPhiStarMin[1] = c.mDPhiStarMin[1];
  92. mVerbose = c.mVerbose;
  93. }
  94. //_________________
  95. MpdFemtoBasicPairCut& MpdFemtoBasicPairCut::operator=(const MpdFemtoBasicPairCut& c) {
  96. // Assignment operator
  97. if (this != &c) {
  98. MpdFemtoBasePairCut::operator=(c);
  99. mNPairsPassed = 0;
  100. mNPairsFailed = 0;
  101. mPrimaryVertex = c.mPrimaryVertex;
  102. mQuality[0] = c.mQuality[0];
  103. mQuality[1] = c.mQuality[1];
  104. mKt[0] = c.mKt[0];
  105. mKt[1] = c.mKt[1];
  106. mPt[0] = c.mPt[0];
  107. mPt[1] = c.mPt[1];
  108. mOpeningAngle[0] = c.mOpeningAngle[0];
  109. mOpeningAngle[1] = c.mOpeningAngle[1];
  110. mQinv[0] = c.mQinv[0];
  111. mQinv[1] = c.mQinv[1];
  112. mMinv[0] = c.mMinv[0];
  113. mMinv[1] = c.mMinv[1];
  114. mRapidity[0] = c.mRapidity[0];
  115. mRapidity[1] = c.mRapidity[1];
  116. mEta[0] = c.mEta[0];
  117. mEta[1] = c.mEta[1];
  118. mEntranceSeparation[0] = c.mEntranceSeparation[0];
  119. mEntranceSeparation[1] = c.mEntranceSeparation[1];
  120. mAngleToPrimaryVertex[0] = c.mAngleToPrimaryVertex[0];
  121. mAngleToPrimaryVertex[1] = c.mAngleToPrimaryVertex[1];
  122. mFracOfMergedRow[0] = c.mFracOfMergedRow[0];
  123. mFracOfMergedRow[1] = c.mFracOfMergedRow[1];
  124. mClosestRowAtDCA[0] = c.mClosestRowAtDCA[0];
  125. mClosestRowAtDCA[1] = c.mClosestRowAtDCA[1];
  126. mWeightedAvSep[0] = c.mWeightedAvSep[0];
  127. mWeightedAvSep[1] = c.mWeightedAvSep[1];
  128. mAverageSeparation[0] = c.mAverageSeparation[0];
  129. mAverageSeparation[1] = c.mAverageSeparation[1];
  130. mRValueLo = c.mRValueLo;
  131. mDPhiStarMin[0] = c.mDPhiStarMin[0];
  132. mDPhiStarMin[1] = c.mDPhiStarMin[1];
  133. mVerbose = c.mVerbose;
  134. } // if (this != &c)
  135. return *this;
  136. }
  137. //_________________
  138. MpdFemtoBasicPairCut::~MpdFemtoBasicPairCut() {
  139. // Destructor
  140. /* no-op */
  141. }
  142. //__________________
  143. bool MpdFemtoBasicPairCut::pass(const MpdFemtoPair* pair) {
  144. // Check if pair passes pair-cut
  145. // For charged TRACKS only
  146. float dPhiStarMin = MpdFemtoPair::calculateDPhiStarMin(
  147. pair->track1()->momentum(),
  148. pair->track1()->track()->charge(),
  149. pair->track2()->momentum(),
  150. pair->track2()->track()->charge(),
  151. 0.1, 0.6, 2., /* in meters */
  152. pair->track1()->track()->bField() * 0.1 /* in Tesla */
  153. );
  154. // Selecting dEta_dPhiStarMin square ...
  155. Double_t dEta = pair->deltaEta();
  156. Bool_t isInside = (dEta >= mDEta[0] && dEta <= mDEta[1] && dPhiStarMin >= mDPhiStarMin[0] && dPhiStarMin <= mDPhiStarMin[1]);
  157. bool mGoodPair = (
  158. (pair->quality() >= mQuality[0]) &&
  159. (pair->quality() <= mQuality[1]) &&
  160. (pair->kT() >= mKt[0]) &&
  161. (pair->kT() <= mKt[1]) &&
  162. (pair->pT() >= mPt[0]) &&
  163. (pair->pT() <= mPt[1]) &&
  164. (pair->openingAngle() >= mOpeningAngle[0]) &&
  165. (pair->openingAngle() <= mOpeningAngle[1]) &&
  166. (pair->rapidity() >= mRapidity[0]) &&
  167. (pair->rapidity() <= mRapidity[1]) &&
  168. (pair->eta() >= mEta[0]) &&
  169. (pair->eta() <= mEta[1]) &&
  170. (pair->qInv() >= mQinv[0]) &&
  171. (pair->qInv() <= mQinv[1]) &&
  172. (pair->mInv() >= mMinv[0]) &&
  173. (pair->mInv() <= mMinv[1]) &&
  174. (pair->emissionAngle() >= mAngleToPrimaryVertex[0]) &&
  175. (pair->emissionAngle() <= mAngleToPrimaryVertex[1]) &&
  176. (pair->nominalTpcEntranceSeparation() >= mEntranceSeparation[0]) &&
  177. (pair->nominalTpcEntranceSeparation() <= mEntranceSeparation[1]) &&
  178. (pair->fractionOfMergedRow() >= mFracOfMergedRow[0]) &&
  179. (pair->fractionOfMergedRow() <= mFracOfMergedRow[1]) &&
  180. (pair->closestRowAtDCA() >= mClosestRowAtDCA[0]) &&
  181. (pair->closestRowAtDCA() <= mClosestRowAtDCA[1]) &&
  182. (pair->weightedAvSep() >= mWeightedAvSep[0]) &&
  183. (pair->weightedAvSep() <= mWeightedAvSep[1]) &&
  184. (pair->nominalTpcAverageSeparation() >= mAverageSeparation[0]) &&
  185. (pair->nominalTpcAverageSeparation() <= mAverageSeparation[1]) &&
  186. (pair->rValue() >= mRValueLo) &&
  187. !isInside
  188. );
  189. // Print pair parameters if requested
  190. if (mVerbose) {
  191. std::cout << "\n --MpdFemtoBasicPairCut-- \n"
  192. << Form("Pair cut %s \n", (mGoodPair) ? "PASSED" : "FAILED")
  193. << Form("Quality : %4.2f <= %4.2f <= %4.2f \t pass: %s \n",
  194. mQuality[0], pair->quality(), mQuality[0],
  195. (pair->quality() >= mQuality[0] && pair->quality() <= mQuality[1]) ? "true" : "false")
  196. << Form("kT : %4.2f <= %4.2f <= %4.2f \t pass: %s \n",
  197. mKt[0], pair->kT(), mKt[1],
  198. (pair->kT() >= mKt[0] && pair->kT() <= mKt[1]) ? "true" : "false")
  199. << Form("pT : %4.2f <= %4.2f <= %4.2f \t pass: %s \n",
  200. mPt[0], pair->pT(), mPt[1],
  201. (pair->pT() >= mPt[0] && pair->pT() <= mPt[1]) ? "true" : "false")
  202. << Form("Opening angle : %4.2f <= %4.2f <= %4.2f \t pass: %s \n",
  203. mOpeningAngle[0], pair->openingAngle(), mOpeningAngle[1],
  204. (pair->openingAngle() >= mOpeningAngle[0] && pair->openingAngle() <= mOpeningAngle[1]) ?
  205. "true" : "false")
  206. << Form("Rapidity : %4.3f <= %4.3f <= %4.3f \t pass: %s \n",
  207. mRapidity[0], pair->rapidity(), mRapidity[1],
  208. (pair->rapidity() >= mRapidity[0] && pair->rapidity() <= mRapidity[1]) ? "true" : "false")
  209. << Form("Eta : %4.3f <= %4.3f <= %4.3f \t pass: %s \n",
  210. mEta[0], pair->eta(), mEta[1],
  211. (pair->eta() >= mEta[0] && pair->eta() <= mEta[1]) ? "true" : "false")
  212. << Form("qInv : %4.3f <= %4.3f <= %4.3f \t pass: %s \n",
  213. mQinv[0], pair->qInv(), mQinv[1],
  214. (pair->qInv() >= mQinv[0] && pair->qInv() <= mQinv[1]) ? "true" : "false")
  215. << Form("mInv : %4.3f <= %4.3f <= %4.3f \t pass: %s \n",
  216. mMinv[0], pair->mInv(), mMinv[1],
  217. (pair->mInv() >= mMinv[0] && pair->mInv() <= mMinv[1]) ? "true" : "false")
  218. << Form("Emission angle : %4.2f <= %4.2f <= %4.2f \t pass: %s \n",
  219. mAngleToPrimaryVertex[0], pair->emissionAngle(), mAngleToPrimaryVertex[1],
  220. (pair->emissionAngle() >= mAngleToPrimaryVertex[0] &&
  221. pair->emissionAngle() <= mAngleToPrimaryVertex[1]) ? "true" : "false")
  222. << Form("TpcEntranceSep : %5.2f <= %5.2f <= %5.2f \t pass: %s \n",
  223. mEntranceSeparation[0], pair->nominalTpcEntranceSeparation(), mEntranceSeparation[1],
  224. (pair->nominalTpcEntranceSeparation() >= mEntranceSeparation[0] &&
  225. pair->nominalTpcEntranceSeparation() <= mEntranceSeparation[1]) ? "true" : "false")
  226. << Form("FracOfMergedRow: %4.3f <= %4.3f <= %4.3f \t pass: %s \n",
  227. mFracOfMergedRow[0], pair->fractionOfMergedRow(), mFracOfMergedRow[1],
  228. (pair->fractionOfMergedRow() >= mFracOfMergedRow[0] &&
  229. pair->fractionOfMergedRow() <= mFracOfMergedRow[1]) ? "true" : "false")
  230. << Form("ClosesRowAtDCA : %4.3f <= %4.3f <= %4.3f \t pass: %s \n",
  231. mClosestRowAtDCA[0], pair->closestRowAtDCA(), mClosestRowAtDCA[1],
  232. (pair->closestRowAtDCA() >= mClosestRowAtDCA[0] &&
  233. pair->closestRowAtDCA() <= mClosestRowAtDCA[1]) ? "true" : "false")
  234. << Form("WeightedAverSep: %5.2f <= %5.2f <= %5.2f \t pass: %s \n",
  235. mWeightedAvSep[0], pair->weightedAvSep(), mWeightedAvSep[1],
  236. (pair->weightedAvSep() >= mWeightedAvSep[0] &&
  237. pair->weightedAvSep() <= mWeightedAvSep[1]) ? "true" : "false")
  238. << Form("AverageSep : %5.2f <= %5.2f <= %5.2f \t pass: %s \n",
  239. mAverageSeparation[0], pair->nominalTpcAverageSeparation(), mAverageSeparation[1],
  240. (pair->nominalTpcAverageSeparation() >= mAverageSeparation[0] &&
  241. pair->nominalTpcAverageSeparation() <= mAverageSeparation[1]) ? "true" : "false")
  242. << Form("R value : %4.2f <= %4.2f \t pass: %s \n",
  243. mRValueLo, pair->rValue(), (pair->rValue() >= mRValueLo) ? "true" : "false")
  244. << Form("dPhi* min : %4.3f <= %4.3f <= %4.3f \t pass: %s \n",
  245. mDPhiStarMin[0], dPhiStarMin, mDPhiStarMin[1],
  246. (dPhiStarMin >= mDPhiStarMin[0] &&
  247. dPhiStarMin <= mDPhiStarMin[1]) ? "true" : "false")
  248. << std::endl << std::endl;
  249. } // if ( mVerbose )
  250. mGoodPair ? mNPairsPassed++ : mNPairsFailed++;
  251. return mGoodPair;
  252. }
  253. //__________________
  254. MpdFemtoString MpdFemtoBasicPairCut::report() {
  255. // Construct the report
  256. TString report;
  257. // Interesting option to print values
  258. // #define PRINTVAR(var) Form("%.4e <= " #var " <= %.4e\n", var[0], var[1])
  259. // report += PRINTVAR(mQuality);
  260. // #undef PRINTVAR
  261. report = "\n-- MpdFemtoBasicPairCut --\n";
  262. report += Form("%5.2f <= quality <= %5.2f\n", mQuality[0], mQuality[1]);
  263. report += Form("%5.2f <= kT <= %5.2f\n", mKt[0], mKt[1]);
  264. report += Form("%5.2f <= pT <= %5.2f\n", mPt[0], mPt[1]);
  265. report += Form("%5.2f <= opening angle <= %5.2f\n", mOpeningAngle[0], mOpeningAngle[1]);
  266. report += Form("%5.2f <= rapidity <= %5.2f\n", mRapidity[0], mRapidity[1]);
  267. report += Form("%5.2f <= eta <= %5.2f\n", mEta[0], mEta[1]);
  268. report += Form("%5.2f <= qInv <= %5.2f\n", mQinv[0], mQinv[1]);
  269. report += Form("%5.2f <= mInv <= %5.2f\n", mMinv[0], mMinv[1]);
  270. report += Form("%5.2f <= eta <= %5.2f\n", mEta[0], mEta[1]);
  271. report += Form("%5.2f <= entrance separation <= %5.2f\n", mEntranceSeparation[0], mEntranceSeparation[1]);
  272. report += Form("%5.2f <= angle to primary vertex <= %5.2f\n", mAngleToPrimaryVertex[0], mAngleToPrimaryVertex[1]);
  273. report += Form("%5.2f <= FMR <= %5.2f\n", mFracOfMergedRow[0], mFracOfMergedRow[1]);
  274. report += Form("%5.2f <= closest row at DCA <= %5.2f\n", mClosestRowAtDCA[0], mClosestRowAtDCA[1]);
  275. report += Form("%5.2f <= weighted average separation <= %5.2f\n", mWeightedAvSep[0], mWeightedAvSep[1]);
  276. report += Form("%5.2f <= average separation <= %5.2f\n", mAverageSeparation[0], mAverageSeparation[1]);
  277. report += Form("R >= %5.2f\n", mRValueLo);
  278. report += Form("%5.2f <= dphi* min <= %5.2f\n", mDPhiStarMin[0], mDPhiStarMin[1]);
  279. report += Form("NPairsPassed = %li\nNPairsFailed = %li\n", mNPairsPassed, mNPairsFailed);
  280. return MpdFemtoString((const char *) report);
  281. }
  282. //_________________
  283. TList *MpdFemtoBasicPairCut::listSettings() {
  284. // Return a list of settings in a writable form
  285. TList *settings_list = new TList();
  286. #define NEWOBJ(var) new TObjString( TString::Format( "MpdFemtoBasicPairCut." #var ".min=%f", var[0] ) ), \
  287. new TObjString( TString::Format( "MpdFemtoBasicPairCut." #var ".max=%f", var[1] ) )
  288. settings_list->AddVector(NEWOBJ(mQuality),
  289. NEWOBJ(mKt),
  290. NEWOBJ(mPt),
  291. NEWOBJ(mOpeningAngle),
  292. NEWOBJ(mRapidity),
  293. NEWOBJ(mEta),
  294. NEWOBJ(mQinv),
  295. NEWOBJ(mMinv),
  296. NEWOBJ(mEntranceSeparation),
  297. NEWOBJ(mAverageSeparation),
  298. new TObjString(TString::Format("MpdFemtoBasicPairCut.mRValueLo=%f", mRValueLo)),
  299. NULL);
  300. return settings_list;
  301. }