MpdFemtoMcDstReader.cxx 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560
  1. //
  2. // Reader for the McDst format that takes McDstReader
  3. //
  4. // MpdFemtoMaker headers
  5. // Base
  6. #include "MpdFemtoBaseEventCut.h"
  7. #include "MpdFemtoBaseTrackCut.h"
  8. #include "MpdFemtoBaseV0Cut.h"
  9. #include "MpdFemtoBaseKinkCut.h"
  10. #include "MpdFemtoBaseXiCut.h"
  11. // Infrastructure
  12. #include "MpdFemtoModelHiddenInfo.h"
  13. // MpdFemtoMakerUser headers
  14. #include "MpdFemtoMcDstReader.h"
  15. // ROOT headers
  16. #include "TMath.h"
  17. #include "TMatrixDSymEigen.h"
  18. #include "TVectorT.h"
  19. ClassImp(MpdFemtoMcDstReader)
  20. //_________________
  21. MpdFemtoMcDstReader::MpdFemtoMcDstReader() :
  22. MpdFemtoBaseEventReader(),
  23. mMcDstReader(nullptr),
  24. mMagField(-5.),
  25. mHbtEvent(nullptr),
  26. mDoRotate(false),
  27. mEventPlaneResolution(1.),
  28. mRefMult(0),
  29. mRefMultPos(0),
  30. mRefMult2(0),
  31. mRefMult2Pos(0),
  32. mSphericity(-1.),
  33. mSphericity2(-1.),
  34. mEventsPassed(0),
  35. mMatrix(nullptr),
  36. mMatrix2(nullptr) {
  37. // Default constructor
  38. mReaderStatus = 0;
  39. mRandom = new TRandom3(0);
  40. mPhi[0] = -100. * TMath::Pi();
  41. mPhi[1] = 100. * TMath::Pi();
  42. mMatrix = new TMatrixTSym<double>(2);
  43. mMatrix2 = new TMatrixTSym<double>(2);
  44. }
  45. //_________________
  46. MpdFemtoMcDstReader::MpdFemtoMcDstReader(McDstReader *reader, int debug) :
  47. MpdFemtoBaseEventReader(),
  48. mMcDstReader(reader),
  49. mMagField(-5.),
  50. mHbtEvent(nullptr),
  51. mDoRotate(false),
  52. mEventPlaneResolution(1.),
  53. mRefMult(0),
  54. mRefMultPos(0),
  55. mRefMult2(0),
  56. mRefMult2Pos(0),
  57. mSphericity(-1.),
  58. mSphericity2(-1.),
  59. mEventsPassed(0),
  60. mMatrix(nullptr),
  61. mMatrix2(nullptr) {
  62. // Parametrized constructor
  63. mRandom = new TRandom3(0);
  64. // mDebug (defined in the base class)
  65. mDebug = debug;
  66. mPhi[0] = -100. * TMath::Pi();
  67. mPhi[1] = 100. * TMath::Pi();
  68. mMatrix = new TMatrixTSym<double>(2);
  69. mMatrix2 = new TMatrixTSym<double>(2);
  70. }
  71. //_________________
  72. MpdFemtoMcDstReader::MpdFemtoMcDstReader(const MpdFemtoMcDstReader& copy) :
  73. MpdFemtoBaseEventReader(),
  74. mMcDstReader(copy.mMcDstReader),
  75. mMagField(copy.mMagField),
  76. mHbtEvent(nullptr),
  77. mDoRotate(copy.mDoRotate),
  78. mEventPlaneResolution(copy.mEventPlaneResolution),
  79. mRefMult(0),
  80. mRefMultPos(0),
  81. mRefMult2(0),
  82. mRefMult2Pos(0),
  83. mSphericity(-1.),
  84. mSphericity2(-1.),
  85. mEventsPassed(0),
  86. mMatrix(nullptr),
  87. mMatrix2(nullptr) {
  88. // Copy constructor
  89. if (mRandom) {
  90. delete mRandom;
  91. mRandom = new TRandom3(0);
  92. }
  93. mPhi[0] = copy.mPhi[0];
  94. mPhi[1] = copy.mPhi[1];
  95. mMatrix = new TMatrixTSym<double>(2);
  96. mMatrix2 = new TMatrixTSym<double>(2);
  97. }
  98. //_________________
  99. MpdFemtoMcDstReader& MpdFemtoMcDstReader::operator=(const MpdFemtoMcDstReader& copy) {
  100. // Assignment operator
  101. if (this != &copy) {
  102. mMcDstReader = copy.mMcDstReader;
  103. mMagField = copy.mMagField;
  104. mHbtEvent = nullptr;
  105. mDoRotate = copy.mDoRotate;
  106. if (mRandom) {
  107. delete mRandom;
  108. mRandom = new TRandom3(0);
  109. }
  110. mPhi[0] = copy.mPhi[0];
  111. mPhi[1] = copy.mPhi[1];
  112. mEventPlaneResolution = copy.mEventPlaneResolution;
  113. mRefMult = 0;
  114. mRefMultPos = 0;
  115. mRefMult2 = 0;
  116. mRefMult2Pos = 0;
  117. mSphericity = -1.;
  118. mSphericity2 = -1.;
  119. mEventsPassed = 0;
  120. if (mMatrix) {
  121. delete mMatrix;
  122. mMatrix = new TMatrixTSym<double>(2);
  123. }
  124. if (mMatrix2) {
  125. delete mMatrix2;
  126. mMatrix = new TMatrixTSym<double>(2);
  127. }
  128. }
  129. return *this;
  130. }
  131. //_________________
  132. MpdFemtoMcDstReader::~MpdFemtoMcDstReader() {
  133. // Destructor
  134. if (mMcDstReader) {
  135. delete mMcDstReader;
  136. mMcDstReader = nullptr;
  137. }
  138. if (mHbtEvent) {
  139. delete mHbtEvent;
  140. mHbtEvent = nullptr;
  141. }
  142. if (mRandom) {
  143. delete mRandom;
  144. mRandom = nullptr;
  145. }
  146. if (mMatrix) {
  147. delete mMatrix;
  148. mMatrix = nullptr;
  149. }
  150. if (mMatrix2) {
  151. delete mMatrix2;
  152. mMatrix = nullptr;
  153. }
  154. }
  155. //_________________
  156. MpdFemtoEvent* MpdFemtoMcDstReader::returnHbtEvent() {
  157. // Read McDst and construct MpdFemtoEvent
  158. // Clean-up mHbtEvent from previous reading
  159. mHbtEvent = nullptr;
  160. // Check that McDstReader exists
  161. if (!mMcDstReader) {
  162. std::cout << "[ERROR] MpdFemtoEvent* MpdFemtoMcDstReader::returnHbtEvent() - no McDstReader is provided" << std::endl;
  163. mReaderStatus = 1;
  164. return mHbtEvent;
  165. }
  166. // Load event
  167. mMcDstReader->loadEntry(mEventsPassed);
  168. // Increment counter
  169. mEventsPassed++;
  170. // Check that McDst exists
  171. McDst *mcDst = mMcDstReader->mcDst();
  172. if (!mcDst) {
  173. std::cout << "[ERROR] MpdFemtoEvent* MpdFemtoMcDstReader::returnHbtEvent() - no McDst is in McDstReader" << std::endl;
  174. mReaderStatus = 1;
  175. return mHbtEvent;
  176. }
  177. // Check that McEvent exists
  178. if (!mcDst->event() || mcDst->numberOfParticles() <= 0) {
  179. if (!mcDst->event()) {
  180. std::cout << "[WARNING] MpdFemtoEvent* MpdFemtoMcDstReader::returnHbtEvent()"
  181. << " - McEvent does not exist" << std::endl;
  182. }
  183. if (mcDst->numberOfParticles() <= 0) {
  184. std::cout << "[WARNING] MpdFemtoEvent* MpdFemtoMcDstReader::returnHbtEvent()"
  185. << " - No particles in the event were found" << std::endl;
  186. }
  187. return mHbtEvent;
  188. }
  189. // Create empty HBT event and then fill it
  190. mHbtEvent = new MpdFemtoEvent();
  191. // Retrieve Monte Carlo event
  192. McEvent *mcEvent = (McEvent*) mcDst->event();
  193. // Next matrices will be used for the sphericity calculation
  194. mMatrix->Zero();
  195. mMatrix2->Zero();
  196. mHbtEvent->setEventNumber(mcEvent->eventNr());
  197. mHbtEvent->setRunNumber(0);
  198. mHbtEvent->setMagneticField(mMagField);
  199. mHbtEvent->setNumberOfPrimaryTracks(0);
  200. mHbtEvent->setCent16(-1);
  201. mHbtEvent->setImpactParameter(mcEvent->impact());
  202. double azimuthalAngle = 0.;
  203. double phi = mcEvent->phi();
  204. // Generate rotation angle by request
  205. if (mDoRotate) {
  206. phi = mRandom->Uniform(mPhi[0], mPhi[1]);
  207. }
  208. mHbtEvent->setEventPlaneAngle(phi);
  209. // Set ideal resolution
  210. mHbtEvent->setEventPlaneResolution(mEventPlaneResolution);
  211. // Initialize values
  212. mRefMult = 0;
  213. mRefMultPos = 0;
  214. mRefMult2 = 0;
  215. mRefMult2Pos = 0;
  216. mSphericity = -1.;
  217. mSphericity2 = -1.;
  218. float pTsum = 0.;
  219. float pTsum2 = 0.;
  220. float px = 0.;
  221. float py = 0.;
  222. float pz = 0.;
  223. float pt = 0.;
  224. // Track loop
  225. for (unsigned int iTrk = 0; iTrk < mcDst->numberOfParticles(); iTrk++) {
  226. // Retrieve i-th particle
  227. McParticle *particle = (McParticle*) mcDst->particle(iTrk);
  228. if (!particle) {
  229. std::cout << "[WARNING] MpdFemtoEvent* MpdFemtoMcDstReader::returnHbtEvent() - Particle does not exist"
  230. << std::endl;
  231. continue;
  232. } // if(!trk)
  233. // Create new MpdFemtoTrack and fill it with required information
  234. MpdFemtoTrack *hbtTrack = new MpdFemtoTrack();
  235. hbtTrack->setId(particle->index());
  236. // Set number of hits (45 correspond to the number
  237. // of pad row in STAR TPC)
  238. // IMPORTANT!!!! nHits also stores charge as: nHits * charge
  239. // but charge in the TDataBasePDF is 3 for positive and -3 for negative
  240. if (particle->charge() > 0) {
  241. hbtTrack->setNHits(45);
  242. } else if (particle->charge() == 0) {
  243. hbtTrack->setNHits(0);
  244. } else {
  245. hbtTrack->setNHits(-45);
  246. }
  247. hbtTrack->setNHitsPossible(45);
  248. hbtTrack->setNHitsDedx(45);
  249. hbtTrack->setChi2(1.);
  250. hbtTrack->setDedx(dEdxMean(particle->mass(), particle->ptot()));
  251. if (TMath::Abs(particle->pdg()) == 211) {
  252. // Pion
  253. hbtTrack->setNSigmaElectron(-65.);
  254. hbtTrack->setNSigmaPion(0.);
  255. hbtTrack->setNSigmaKaon(65.);
  256. hbtTrack->setNSigmaProton(65.);
  257. hbtTrack->setPidProbElectron(0.);
  258. hbtTrack->setPidProbPion(1.);
  259. hbtTrack->setPidProbKaon(0.);
  260. hbtTrack->setPidProbProton(0.);
  261. } else if (TMath::Abs(particle->pdg()) == 321) {
  262. // Kaon
  263. hbtTrack->setNSigmaElectron(-65.);
  264. hbtTrack->setNSigmaPion(-65.);
  265. hbtTrack->setNSigmaKaon(0.);
  266. hbtTrack->setNSigmaProton(65.);
  267. hbtTrack->setPidProbElectron(0.);
  268. hbtTrack->setPidProbPion(0.);
  269. hbtTrack->setPidProbKaon(1.);
  270. hbtTrack->setPidProbProton(0.);
  271. } else if (TMath::Abs(particle->pdg()) == 2212) {
  272. // Proton
  273. hbtTrack->setNSigmaElectron(-65.);
  274. hbtTrack->setNSigmaPion(-65.);
  275. hbtTrack->setNSigmaKaon(-65.);
  276. hbtTrack->setNSigmaProton(0.);
  277. hbtTrack->setPidProbElectron(0.);
  278. hbtTrack->setPidProbPion(0.);
  279. hbtTrack->setPidProbKaon(0.);
  280. hbtTrack->setPidProbProton(1.);
  281. } else if (TMath::Abs(particle->pdg()) == 11) {
  282. // Electron
  283. hbtTrack->setNSigmaElectron(0.);
  284. hbtTrack->setNSigmaPion(65.);
  285. hbtTrack->setNSigmaKaon(65.);
  286. hbtTrack->setNSigmaProton(65.);
  287. hbtTrack->setPidProbElectron(1.);
  288. hbtTrack->setPidProbPion(0.);
  289. hbtTrack->setPidProbKaon(0.);
  290. hbtTrack->setPidProbProton(0.);
  291. } else {
  292. // Other particles
  293. hbtTrack->setNSigmaElectron(-65.);
  294. hbtTrack->setNSigmaPion(-65.);
  295. hbtTrack->setNSigmaKaon(-65.);
  296. hbtTrack->setNSigmaProton(-65.);
  297. hbtTrack->setPidProbElectron(0.);
  298. hbtTrack->setPidProbPion(0.);
  299. hbtTrack->setPidProbKaon(0.);
  300. hbtTrack->setPidProbProton(0.);
  301. }
  302. // McParticle stores emission point (last or decay point) in fm,
  303. // but DCA is stored in cm
  304. hbtTrack->setDca(particle->x() * 1e-13,
  305. particle->y() * 1e-13,
  306. particle->x() * 1e-13);
  307. px = particle->px();
  308. py = particle->py();
  309. pz = particle->pz();
  310. pt = particle->pt();
  311. if (mDoRotate) {
  312. hbtTrack->setPz(pz);
  313. hbtTrack->setGlobalPz(pz);
  314. azimuthalAngle = particle->phi();
  315. azimuthalAngle += phi;
  316. hbtTrack->setPx(pt * TMath::Cos(azimuthalAngle));
  317. hbtTrack->setPy(pt * TMath::Sin(azimuthalAngle));
  318. hbtTrack->setGlobalPx(pt * TMath::Cos(azimuthalAngle));
  319. hbtTrack->setGlobalPy(pt * TMath::Sin(azimuthalAngle));
  320. } else {
  321. hbtTrack->setP(px, py, pz);
  322. hbtTrack->setGlobalP(px, pt, pz);
  323. }
  324. // McDst assumes that event happen at (0.,0.,0.)
  325. hbtTrack->setPrimaryVertex(0., 0., 0.);
  326. hbtTrack->setMagneticField(mMagField);
  327. // Set topology map (probably can insert probabilty function from data)
  328. hbtTrack->setTopologyMap(0);
  329. hbtTrack->setBeta(particle->ptot() / particle->energy());
  330. // Create model hidden info and fill it
  331. MpdFemtoModelHiddenInfo *hiddenInfo = new MpdFemtoModelHiddenInfo();
  332. // px and py values are not a mistake and take into
  333. // account event plane angle rotation
  334. hiddenInfo->setTrueMomentum(hbtTrack->p().X(),
  335. hbtTrack->p().Y(),
  336. pz);
  337. // Check if rotate event plane angle (IS IT NEEDED??)
  338. if (mDoRotate) {
  339. double radialPos = TMath::Sqrt(particle->x() * particle->x() +
  340. particle->y() * particle->y());
  341. azimuthalAngle = TMath::ATan2(particle->y(), particle->x());
  342. azimuthalAngle += phi;
  343. hiddenInfo->setEmissionPoint(radialPos * TMath::Cos(azimuthalAngle),
  344. radialPos * TMath::Sin(azimuthalAngle),
  345. particle->z(),
  346. particle->t());
  347. } else {
  348. hiddenInfo->setEmissionPoint(particle->x(),
  349. particle->y(),
  350. particle->z(),
  351. particle->t());
  352. }
  353. hiddenInfo->setPdgPid(particle->pdg());
  354. hiddenInfo->setOrigin(0);
  355. hbtTrack->setHiddenInfo(hiddenInfo);
  356. // Check if a front-loaded cut exists. The mTrackCut
  357. // is inherited from MpdFemtoBaseReader
  358. if (mTrackCut && !mTrackCut->pass(hbtTrack)) {
  359. delete hbtTrack;
  360. continue;
  361. } // if ( mTrackCut )
  362. // Add track to track collection
  363. mHbtEvent->trackCollection()->push_back(hbtTrack);
  364. //
  365. // Calculate event properties: refMult, sphericity
  366. //
  367. // In the experiment, one works only with charged particles
  368. if (particle->charge() == 0) continue;
  369. // Remove nucleons from incoming nuclei (should work for UrQMD)
  370. //if (particle->parent() == 0) continue;
  371. // Particle must have pT>0.15 GeV/c and |eta|<1 (STAR acceptance)
  372. if (TMath::Abs(particle->eta()) > 1. ||
  373. pt < 0.15) continue;
  374. mRefMult2++;
  375. if (particle->charge() > 0) {
  376. mRefMult2Pos++;
  377. }
  378. // Event properties in |eta|<1.
  379. px = particle->px();
  380. py = particle->py();
  381. pt = particle->pt();
  382. (*mMatrix2)(0, 0) += px * px / pt;
  383. (*mMatrix2)(1, 1) += py * py / pt;
  384. (*mMatrix2)(0, 1) += px * py / pt;
  385. (*mMatrix2)(1, 0) += px * py / pt;
  386. pTsum2 += pt;
  387. // Event properties in |eta|<0.5
  388. if (TMath::Abs(particle->eta()) <= 0.5) {
  389. mRefMult++;
  390. if (particle->charge() > 0) {
  391. mRefMultPos++;
  392. }
  393. (*mMatrix)(0, 0) += px * px / pt;
  394. (*mMatrix)(1, 1) += py * py / pt;
  395. (*mMatrix)(0, 1) += px * py / pt;
  396. (*mMatrix)(1, 0) += px * py / pt;
  397. pTsum += pt;
  398. } // if ( TMath::Abs( particle->eta() ) <= 0.5 )
  399. } // for(int iTrk=0; iTrk<mcDst->numberOfParticles(); iTrk++)
  400. if (pTsum != 0) {
  401. *mMatrix *= 1. / pTsum;
  402. TMatrixDSymEigen eigenEstimator(*mMatrix);
  403. TVectorD eigen = eigenEstimator.GetEigenValues();
  404. mSphericity = 2. * eigen.Min() / (eigen[0] + eigen[1]);
  405. }
  406. if (pTsum2 != 0) {
  407. *mMatrix2 *= 1. / pTsum2;
  408. TMatrixDSymEigen eigenEstimator2(*mMatrix2);
  409. TVectorD eigen2 = eigenEstimator2.GetEigenValues();
  410. mSphericity2 = 2. * eigen2.Min() / (eigen2[0] + eigen2[1]);
  411. }
  412. mHbtEvent->setPrimaryVertex(0., 0., 0.);
  413. mHbtEvent->setRefMult(mRefMult);
  414. mHbtEvent->setRefMultPos(mRefMultPos);
  415. mHbtEvent->setRefMult2(mRefMult2);
  416. mHbtEvent->setRefMult2Pos(mRefMult2Pos);
  417. mHbtEvent->setGRefMult(mRefMult);
  418. mHbtEvent->setGRefMultPos(mRefMultPos);
  419. mHbtEvent->setBTofTrayMult(mRefMult);
  420. mHbtEvent->setNumberOfBTofMatched(mRefMult);
  421. mHbtEvent->setNumberOfBEMCMatched(mRefMult);
  422. mHbtEvent->setSphericity(mSphericity);
  423. mHbtEvent->setSphericity2(mSphericity2);
  424. // Check front-loaded cuts here, because now all event information
  425. // should be filled making event cut on it possible. The mEventCut
  426. // is inherited from MpdFemtoBaseReader
  427. if (mEventCut && !mEventCut->pass(mHbtEvent)) {
  428. delete mHbtEvent;
  429. mHbtEvent = nullptr;
  430. // return mHbtEvent; // Will be deleted 6 lines below
  431. } // if (mEventCut)
  432. return mHbtEvent;
  433. }
  434. //_________________
  435. double MpdFemtoMcDstReader::dEdxMean(Double_t mass, Double_t momentum) {
  436. Double_t dedxMean = 0.;
  437. Double_t tpcDedxGain = 0.174325e-06;
  438. Double_t tpcDedxOffset = -2.71889;
  439. Double_t tpcDedxRise = 776.626;
  440. Double_t gamma = TMath::Sqrt(momentum * momentum / (mass * mass) + 1.);
  441. Double_t beta = TMath::Sqrt(1. - 1. / (gamma * gamma));
  442. Double_t rise = tpcDedxRise * beta * beta * gamma*gamma;
  443. if (beta > 0) {
  444. dedxMean = tpcDedxGain / (beta * beta) * (0.5 * TMath::Log(rise) -
  445. beta * beta - tpcDedxOffset);
  446. } else {
  447. dedxMean = 1000.;
  448. }
  449. return dedxMean;
  450. }
  451. //_________________
  452. MpdFemtoString MpdFemtoMcDstReader::report() {
  453. // Make a report
  454. MpdFemtoString repstr = "\nMpdFemtoString MpdFemtoMcDstReader::report() - reporting\n";
  455. repstr += "---> Event cut in the reader: ";
  456. if (mEventCut) {
  457. repstr += mEventCut->report();
  458. } else {
  459. repstr += "\tNONE";
  460. }
  461. repstr += "\n---> Track cut in the reader: ";
  462. if (mTrackCut) {
  463. repstr += mTrackCut->report();
  464. } else {
  465. repstr += "\tNONE";
  466. }
  467. repstr += "\n---> V0 cut in the reader: ";
  468. if (mV0Cut) {
  469. repstr += mV0Cut->report();
  470. } else {
  471. repstr += "\tNONE";
  472. }
  473. repstr += "\n---> V0 cut in the reader: ";
  474. if (mV0Cut) {
  475. repstr += mV0Cut->report();
  476. } else {
  477. repstr += "\tNONE";
  478. }
  479. repstr += "\n---> Kink cut in the reader: ";
  480. if (mKinkCut) {
  481. repstr += mKinkCut->report();
  482. } else {
  483. repstr += "\tNONE";
  484. }
  485. repstr += "\n---> Xi cut in the reader: ";
  486. if (mXiCut) {
  487. repstr += mXiCut->report();
  488. } else {
  489. repstr += "\tNONE";
  490. }
  491. repstr += "\n";
  492. return repstr;
  493. }