MpdStack.cxx 22 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691
  1. // -------------------------------------------------------------------------
  2. // ----- FairStack source file -----
  3. // ----- Created 10/08/04 by D. Bertini / V. Friese -----
  4. // ----- adopted for NICA/MPD 29/03/10 (litvin) -----
  5. // ----- adopted for NICA/MPD 20/12/19 (ABychkov) -----
  6. // ----- added external decayer 30/04/20 (AZinchenko) -----
  7. // -------------------------------------------------------------------------
  8. #include "MpdStack.h"
  9. #include "MpdDecayer.h"
  10. #include "MpdDecayerPyt8.h"
  11. #include "MpdMCTrack.h"
  12. #include "FairDetector.h"
  13. #include "FairMCPoint.h"
  14. #include "FairRootManager.h"
  15. #include "TError.h"
  16. #include "TLorentzVector.h"
  17. #include "TParticle.h"
  18. #include "TRefArray.h"
  19. #include "TParticlePDG.h"
  20. #include "TVirtualMC.h"
  21. #include <list>
  22. #include <iostream>
  23. using std::cout;
  24. using std::endl;
  25. using std::pair;
  26. using std::vector;
  27. // ----- Default constructor -------------------------------------------
  28. MpdStack::MpdStack(Int_t size) // {
  29. // fParticles = new TClonesArray("TParticle", size);
  30. // fTracks = new TClonesArray("MpdMCTrack", size);
  31. // fCurrentTrack = -1;
  32. // fNPrimaries = fNParticles = fNTracks = 0;
  33. // fIndex = 0;
  34. // fStoreSecondaries = kTRUE;
  35. // fMinPoints = 1;
  36. // fEnergyCut = 0.;
  37. // fStoreMothers = kTRUE;
  38. //}
  39. : FairGenericStack(),
  40. fStack(),
  41. fParticles(new TClonesArray("TParticle", size)),
  42. fTracks(new TClonesArray("MpdMCTrack", size)),
  43. fStoreMap(),
  44. fStoreIter(),
  45. fIndexMap(),
  46. fIndexIter(),
  47. fPointsMap(),
  48. fCurrentTrack(-1),
  49. fNPrimaries(0),
  50. fNParticles(0),
  51. fNTracks(0),
  52. fIndex(0),
  53. fDecayFlag(0),
  54. fStoreSecondaries(kTRUE),
  55. fMinPoints(1),
  56. //fEnergyCut(0.),
  57. fEnergyCut(-0.001), //AZ
  58. fStoreMothers(kTRUE),
  59. fMinMotherMass(0.4),
  60. fMaxMotherMass(6.1),
  61. fRadiusCut(172), //175 NG ecal outer XYradius, some mm smaller
  62. fVzCut(311), //296 NG barrel half length, TODO does not include EndCaps
  63. fNoZDC(kTRUE) //NG
  64. {
  65. fDecays = new TClonesArray("TParticle");
  66. }
  67. // -------------------------------------------------------------------------
  68. // ----- Destructor ----------------------------------------------------
  69. MpdStack::~MpdStack() {
  70. if (fParticles) {
  71. fParticles->Delete();
  72. delete fParticles;
  73. }
  74. if (fTracks) {
  75. fTracks->Delete();
  76. delete fTracks;
  77. }
  78. if (fDecays) {
  79. fDecays->Delete();
  80. delete fDecays;
  81. }
  82. }
  83. // -------------------------------------------------------------------------
  84. // ----- Virtual public method PushTrack -------------------------------
  85. void MpdStack::PushTrack(Int_t toBeDone, Int_t parentId, Int_t pdgCode,
  86. Double_t px, Double_t py, Double_t pz,
  87. Double_t e, Double_t vx, Double_t vy, Double_t vz,
  88. Double_t time, Double_t polx, Double_t poly,
  89. Double_t polz, TMCProcess proc, Int_t& ntr,
  90. Double_t weight, Int_t is) {
  91. // --> Get TParticle array
  92. TClonesArray& partArray = *fParticles;
  93. // --> Create new TParticle and add it to the TParticle array
  94. Int_t trackId = fNParticles;
  95. Int_t nPoints = 0;
  96. Int_t daughter1Id = -1;
  97. Int_t daughter2Id = -1;
  98. TParticle* particle =
  99. new(partArray[fNParticles++]) TParticle(pdgCode, trackId, parentId,
  100. nPoints, daughter1Id,
  101. daughter2Id, px, py, pz, e,
  102. vx, vy, vz, time);
  103. particle->SetPolarisation(polx, poly, polz);
  104. particle->SetWeight(weight);
  105. particle->SetUniqueID(proc);
  106. //AZ - global polarization transfer to decay Lambda
  107. if (TMath::Abs(pdgCode) == 3122 && proc == kPDecay) {
  108. TParticle *moth = (TParticle*) partArray.UncheckedAt(parentId);
  109. TVector3 polar;
  110. moth->GetPolarisation(polar);
  111. Double_t value = moth->GetWeight();
  112. if (TMath::Abs(moth->GetPdgCode()) == 3312) value *= 0.927; // Xi+-
  113. else if (TMath::Abs(moth->GetPdgCode()) == 3322) value *= 0.900; // Xi0
  114. else if (TMath::Abs(moth->GetPdgCode()) == 3212) {
  115. // Sigma0
  116. value *= (1./3);
  117. polar *= -1.0;
  118. }
  119. particle->SetPolarisation(polar);
  120. particle->SetWeight(value);
  121. }
  122. // --> Increment counter
  123. if (parentId < 0) fNPrimaries++;
  124. // --> Set argument variable
  125. ntr = trackId;
  126. //AZ
  127. if (fDecayFlag && (!particle->GetPDG()->Stable() || (particle->GetPdgCode() == 221) || (particle->GetPdgCode() == 111)) && particle->GetPDG()->DecayList()) {
  128. // Decay unstable particle
  129. //MpdDecayer::Instance()->Decay(particle);
  130. //Int_t npart = MpdDecayer::Instance()->ImportParticles(fDecays);
  131. MpdDecayerPyt8::Instance()->Decay(particle);
  132. Int_t npart = MpdDecayerPyt8::Instance()->ImportParticles(fDecays);
  133. // Copy decay products to vector to avoid TClonesArray overwriting
  134. // for recursive decays
  135. vector<TParticle> vDecays;
  136. for (Int_t j = 0; j < npart; ++j) {
  137. TParticle tpart = *((TParticle*) fDecays->UncheckedAt(j));
  138. tpart.Print();
  139. vDecays.push_back(tpart);
  140. }
  141. for (Int_t j = 0; j < npart; ++j) { // skip mother particle
  142. if (j == 0) {
  143. particle->SetStatusCode(11); // decayed particle
  144. continue; // skip mother particle
  145. }
  146. //TParticle *part = (TParticle*) fDecays->UncheckedAt(j);
  147. TParticle *part = &vDecays[j];
  148. Int_t toDo = 1;
  149. if (TString(gMC->GetName()) == "TGeant4") toDo = 0;
  150. //PushTrack(toBeDone, trackId, part->GetPdgCode(),
  151. PushTrack(toDo, trackId, part->GetPdgCode(),
  152. part->Px(), part->Py(), part->Pz(),
  153. part->Energy(), part->Vx(), part->Vy(), part->Vz(),
  154. time, polx, poly, polz, proc, ntr,
  155. weight, is);
  156. }
  157. if (parentId < 0) fNPrimaries += (npart-1); // treat decay products as primaries (dirty trick)
  158. } else {
  159. // --> Push particle on the stack if toBeDone is set
  160. if (toBeDone == 1) fStack.push(particle);
  161. }
  162. // --> Push particle on the stack if toBeDone is set
  163. //AZ if (toBeDone == 1) fStack.push(particle);
  164. }
  165. void MpdStack::PushTrack(Int_t toBeDone, Int_t parentId, Int_t pdgCode,
  166. Double_t px, Double_t py, Double_t pz,
  167. Double_t e, Double_t vx, Double_t vy, Double_t vz,
  168. Double_t time, Double_t polx, Double_t poly,
  169. Double_t polz, TMCProcess proc, Int_t& ntr,
  170. Double_t weight, Int_t is, Int_t secondMotherID) {
  171. // --> Get TParticle array
  172. TClonesArray& partArray = *fParticles;
  173. // --> Create new TParticle and add it to the TParticle array
  174. Int_t trackId = fNParticles;
  175. Int_t nPoints = 0;
  176. Int_t daughter1Id = -1;
  177. Int_t daughter2Id = -1;
  178. TParticle* particle =
  179. // new(partArray[fNParticles++]) TParticle(pdgCode, trackId, parentId,
  180. new(partArray[fNParticles++]) TParticle(pdgCode, trackId, secondMotherID,
  181. nPoints, daughter1Id,
  182. daughter2Id, px, py, pz, e,
  183. vx, vy, vz, time);
  184. particle->SetPolarisation(polx, poly, polz);
  185. particle->SetWeight(weight);
  186. particle->SetUniqueID(proc);
  187. //AZ - global polarization transfer to decay Lambda
  188. if (TMath::Abs(pdgCode) == 3122 && proc == kPDecay) {
  189. TParticle *moth = (TParticle*) partArray.UncheckedAt(parentId);
  190. TVector3 polar;
  191. moth->GetPolarisation(polar);
  192. Double_t value = moth->GetWeight();
  193. if (TMath::Abs(moth->GetPdgCode()) == 3312) value *= 0.927; // Xi+-
  194. else if (TMath::Abs(moth->GetPdgCode()) == 3322) value *= 0.900; // Xi0
  195. else if (TMath::Abs(moth->GetPdgCode()) == 3212) {
  196. // Sigma0
  197. value *= (1./3);
  198. polar *= -1.0;
  199. }
  200. particle->SetPolarisation(polar);
  201. particle->SetWeight(value);
  202. }
  203. // --> Increment counter
  204. if (parentId < 0) fNPrimaries++;
  205. // --> Set argument variable
  206. ntr = trackId;
  207. //AZ
  208. if (fDecayFlag && (!particle->GetPDG()->Stable() || (particle->GetPdgCode() == 221) || (particle->GetPdgCode() == 111)) && particle->GetPDG()->DecayList()) {
  209. // Decay unstable particle
  210. //MpdDecayer::Instance()->Decay(particle);
  211. //Int_t npart = MpdDecayer::Instance()->ImportParticles(fDecays);
  212. MpdDecayerPyt8::Instance()->Decay(particle);
  213. Int_t npart = MpdDecayerPyt8::Instance()->ImportParticles(fDecays);
  214. // Copy decay products to vector to avoid TClonesArray overwriting
  215. // for recursive decays
  216. vector<TParticle> vDecays;
  217. for (Int_t j = 0; j < npart; ++j) {
  218. TParticle tpart = *((TParticle*) fDecays->UncheckedAt(j));
  219. tpart.Print();
  220. vDecays.push_back(tpart);
  221. }
  222. for (Int_t j = 0; j < npart; ++j) {
  223. if (j == 0) {
  224. particle->SetStatusCode(11); // decayed particle
  225. continue; // skip mother particle
  226. }
  227. //TParticle *part = (TParticle*) fDecays->UncheckedAt(j);
  228. TParticle *part = &vDecays[j];
  229. Int_t toDo = 1;
  230. if (TString(gMC->GetName()) == "TGeant4") toDo = 0;
  231. //PushTrack(toBeDone, trackId, part->GetPdgCode(),
  232. PushTrack(toDo, trackId, part->GetPdgCode(),
  233. part->Px(), part->Py(), part->Pz(),
  234. part->Energy(), part->Vx(), part->Vy(), part->Vz(),
  235. time, polx, poly, polz, proc, ntr,
  236. weight, is);
  237. }
  238. if (parentId < 0) fNPrimaries += (npart-1); // treat decay products as primaries (dirty trick)
  239. } else {
  240. // --> Push particle on the stack if toBeDone is set
  241. if (toBeDone == 1) fStack.push(particle);
  242. }
  243. // --> Push particle on the stack if toBeDone is set
  244. //AZ if (toBeDone == 1) fStack.push(particle);
  245. }
  246. // -------------------------------------------------------------------------
  247. // ----- Virtual method PopNextTrack -----------------------------------
  248. TParticle* MpdStack::PopNextTrack(Int_t& iTrack) {
  249. // If end of stack: Return empty pointer
  250. if (fStack.empty()) {
  251. iTrack = -1;
  252. return NULL;
  253. }
  254. // If not, get next particle from stack
  255. TParticle* thisParticle = fStack.top();
  256. fStack.pop();
  257. if ( !thisParticle) {
  258. iTrack = 0;
  259. return NULL;
  260. }
  261. fCurrentTrack = thisParticle->GetStatusCode();
  262. iTrack = fCurrentTrack;
  263. return thisParticle;
  264. }
  265. // -------------------------------------------------------------------------
  266. // ----- Virtual method PopPrimaryForTracking --------------------------
  267. TParticle* MpdStack::PopPrimaryForTracking(Int_t iPrim) {
  268. // Get the iPrimth particle from the fStack TClonesArray. This
  269. // should be a primary (if the index is correct).
  270. // Test for index
  271. //AZ-15.05.20 if (iPrim < 0 || iPrim >= fNPrimaries) {
  272. if (iPrim < 0 || iPrim >= fParticles->GetEntriesFast()) {
  273. cout << fParticles->GetEntriesFast() << endl;
  274. cout << "-E- MpdStack: Primary index out of range! " << iPrim << endl;
  275. Fatal("MpdStack::PopPrimaryForTracking", "Index out of range");
  276. }
  277. // Return the iPrim-th TParticle from the fParticle array. This should be
  278. // a primary.
  279. TParticle* part = (TParticle*)fParticles->At(iPrim);
  280. /*AZ-30.04.20
  281. if ( ! (part->GetMother(0) < 0) ) {
  282. cout << "-E- MpdStack:: Not a primary track! " << iPrim << endl;
  283. Fatal("MpdStack::PopPrimaryForTracking", "Not a primary track");
  284. }
  285. */
  286. if (part->GetStatusCode() == 11) return (TParticle*)NULL; //AZ-30.04.20
  287. return part;
  288. }
  289. // -------------------------------------------------------------------------
  290. // ----- Virtual public method GetCurrentTrack -------------------------
  291. TParticle* MpdStack::GetCurrentTrack() const {
  292. TParticle* currentPart = GetParticle(fCurrentTrack);
  293. if ( ! currentPart) {
  294. cout << "-W- MpdStack: Current track not found in stack!" << endl;
  295. Warning("MpdStack::GetCurrentTrack", "Track not found in stack");
  296. }
  297. return currentPart;
  298. }
  299. // -------------------------------------------------------------------------
  300. // ----- Public method AddParticle -------------------------------------
  301. void MpdStack::AddParticle(TParticle* oldPart) {
  302. TClonesArray& array = *fParticles;
  303. TParticle* newPart = new(array[fIndex]) TParticle(*oldPart);
  304. newPart->SetWeight(oldPart->GetWeight());
  305. newPart->SetUniqueID(oldPart->GetUniqueID());
  306. fIndex++;
  307. }
  308. // -------------------------------------------------------------------------
  309. // ----- Public method FillTrackArray ----------------------------------
  310. void MpdStack::FillTrackArray() {
  311. cout << "-I- MpdStack: Filling MCTrack array..." << endl;
  312. // --> Reset index map and number of output tracks
  313. fIndexMap.clear();
  314. fNTracks = 0;
  315. // --> Check tracks for selection criteria
  316. SelectTracks();
  317. // --> Loop over fParticles array and copy selected tracks
  318. for (Int_t iPart=0; iPart<fNParticles; iPart++) {
  319. fStoreIter = fStoreMap.find(iPart);
  320. if (fStoreIter == fStoreMap.end() ) {
  321. cout << "-E- MpdStack: Particle " << iPart
  322. << " not found in storage map!" << endl;
  323. Fatal("MpdStack::FillTrackArray",
  324. "Particle not found in storage map.");
  325. }
  326. Bool_t store = (*fStoreIter).second;
  327. if (store) {
  328. MpdMCTrack* track =
  329. new( (*fTracks)[fNTracks]) MpdMCTrack(GetParticle(iPart));
  330. fIndexMap[iPart] = fNTracks;
  331. // --> Set the number of points in the detectors for this track
  332. for (Int_t iDet=kSTS; iDet<=kMCORD; iDet++) {
  333. pair<Int_t, Int_t> a(iPart, iDet);
  334. track->SetNPoints(iDet, fPointsMap[a]);
  335. }
  336. fNTracks++;
  337. }
  338. else fIndexMap[iPart] = -2;
  339. }
  340. // --> Map index for primary mothers
  341. fIndexMap[-1] = -1;
  342. // --> Screen output
  343. Print(0);
  344. }
  345. // -------------------------------------------------------------------------
  346. // ----- Public method UpdateTrackIndex --------------------------------
  347. void MpdStack::UpdateTrackIndex(TRefArray* detList) {
  348. cout << "-I- MpdStack: Updating track indizes...";
  349. Int_t nColl = 0;
  350. // First update mother ID in MCTracks
  351. for (Int_t i=0; i<fNTracks; i++) {
  352. MpdMCTrack* track = (MpdMCTrack*)fTracks->At(i);
  353. Int_t iMotherOld = track->GetMotherId();
  354. if (iMotherOld < 0) continue; //AZ
  355. fIndexIter = fIndexMap.find(iMotherOld);
  356. if (fIndexIter == fIndexMap.end()) {
  357. cout << "-E- MpdStack: Particle index " << iMotherOld
  358. << " not found in dex map! " << endl;
  359. Fatal("MpdStack::UpdateTrackIndex",
  360. "Particle index not found in map");
  361. }
  362. track->SetMotherId( (*fIndexIter).second );
  363. }
  364. // Now iterate through all active detectors
  365. TIterator* detIter = detList->MakeIterator();
  366. detIter->Reset();
  367. FairDetector* det = NULL;
  368. while( (det = (FairDetector*)detIter->Next() ) ) {
  369. // --> Get hit collections from detector
  370. Int_t iColl = 0;
  371. TClonesArray* hitArray;
  372. while ( (hitArray = det->GetCollection(iColl++)) ) {
  373. nColl++;
  374. Int_t nPoints = hitArray->GetEntriesFast();
  375. // --> Update track index for all MCPoints in the collection
  376. for (Int_t iPoint=0; iPoint<nPoints; iPoint++) {
  377. FairMCPoint* point = (FairMCPoint*)hitArray->At(iPoint);
  378. Int_t iTrack = point->GetTrackID();
  379. fIndexIter = fIndexMap.find(iTrack);
  380. if (fIndexIter == fIndexMap.end()) {
  381. cout << "-E- MpdStack: Particle index " << iTrack
  382. << " not found in index map! " << endl;
  383. Fatal("MpdStack::UpdateTrackIndex",
  384. "Particle index not found in map");
  385. }
  386. point->SetTrackID((*fIndexIter).second);
  387. //point->SetLink(FairLink("MCTrack", (*fIndexIter).second));
  388. }
  389. } // Collections of this detector
  390. } // List of active detectors
  391. cout << "...stack and " << nColl << " collections updated." << endl;
  392. }
  393. // -------------------------------------------------------------------------
  394. // ----- Public method Reset -------------------------------------------
  395. void MpdStack::Reset() {
  396. fIndex = 0;
  397. fCurrentTrack = -1;
  398. fNPrimaries = fNParticles = fNTracks = 0;
  399. while (! fStack.empty() ) fStack.pop();
  400. fParticles->Clear();
  401. fTracks->Clear();
  402. fPointsMap.clear();
  403. }
  404. // -------------------------------------------------------------------------
  405. // ----- Public method Register ----------------------------------------
  406. void MpdStack::Register() {
  407. FairRootManager::Instance()->Register("MCTrack", "Stack", fTracks,kTRUE);
  408. }
  409. // -------------------------------------------------------------------------
  410. // ----- Public method Print --------------------------------------------
  411. void MpdStack::Print(Int_t iVerbose) const {
  412. cout << "-I- MpdStack: Number of primaries = "
  413. << fNPrimaries << endl;
  414. cout << " Total number of particles = "
  415. << fNParticles << endl;
  416. cout << " Number of tracks in output = "
  417. << fNTracks << endl;
  418. if (iVerbose) {
  419. for (Int_t iTrack=0; iTrack<fNTracks; iTrack++)
  420. ((MpdMCTrack*) fTracks->At(iTrack))->Print(iTrack);
  421. }
  422. }
  423. // -------------------------------------------------------------------------
  424. // ----- Public method AddPoint (for current track) --------------------
  425. void MpdStack::AddPoint(DetectorIdMPD detId) {
  426. Int_t iDet = detId;
  427. pair<Int_t, Int_t> a(fCurrentTrack, iDet);
  428. if ( fPointsMap.find(a) == fPointsMap.end() ) fPointsMap[a] = 1;
  429. else fPointsMap[a]++;
  430. }
  431. // -------------------------------------------------------------------------
  432. // ----- Public method AddPoint (for arbitrary track) -------------------
  433. void MpdStack::AddPoint(DetectorIdMPD detId, Int_t iTrack) {
  434. if ( iTrack < 0 ) return;
  435. Int_t iDet = detId;
  436. pair<Int_t, Int_t> a(iTrack, iDet);
  437. if ( fPointsMap.find(a) == fPointsMap.end() ) fPointsMap[a] = 1;
  438. else fPointsMap[a]++;
  439. }
  440. // -------------------------------------------------------------------------
  441. // ----- Virtual method GetCurrentParentTrackNumber --------------------
  442. Int_t MpdStack::GetCurrentParentTrackNumber() const {
  443. TParticle* currentPart = GetCurrentTrack();
  444. if ( currentPart ) return currentPart->GetFirstMother();
  445. else return -1;
  446. }
  447. // -------------------------------------------------------------------------
  448. // ----- Public method GetParticle -------------------------------------
  449. TParticle* MpdStack::GetParticle(Int_t trackID) const {
  450. if (trackID < 0 || trackID >= fNParticles) {
  451. cout << "-E- MpdStack: Particle index " << trackID
  452. << " out of range." << endl;
  453. Fatal("MpdStack::GetParticle", "Index out of range");
  454. }
  455. return (TParticle*)fParticles->At(trackID);
  456. }
  457. // -------------------------------------------------------------------------
  458. // ----- Private method SelectTracks -----------------------------------
  459. void MpdStack::SelectTracks() {
  460. //#define EDEBUG
  461. #ifdef EDEBUG
  462. static Int_t lEDEBUGcounter=0;
  463. if (lEDEBUGcounter<=20)
  464. std::cout << "EDEBUG-- MpdStack::SelectTracks(): entered " << fNParticles << endl;
  465. lEDEBUGcounter++;
  466. #endif
  467. // --> Clear storage map
  468. fStoreMap.clear();
  469. // --> Check particles in the fParticle array
  470. for (Int_t i=0; i<fNParticles; i++) {
  471. TParticle* thisPart = GetParticle(i);
  472. Bool_t store = kTRUE;
  473. #ifdef EDEBUG
  474. if (lEDEBUGcounter<=20)
  475. thisPart->Print();
  476. lEDEBUGcounter++;
  477. #endif
  478. #undef EDEBUG
  479. // --> Get track parameters
  480. Int_t iMother = thisPart->GetMother(0);
  481. TLorentzVector p;
  482. thisPart->Momentum(p);
  483. Double_t energy = p.E();
  484. Double_t mass = thisPart->GetMass();
  485. Double_t eKin = energy - mass;
  486. // --> Calculate number of points
  487. Int_t nPoints = 0, nPointsZDC = 0, nPointsNoZDC = 0;
  488. for (Int_t iDet=kSTS; iDet<=kMCORD; iDet++) {
  489. pair<Int_t, Int_t> a(i, iDet);
  490. if ( fPointsMap.find(a) != fPointsMap.end() ) {
  491. nPoints += fPointsMap[a];
  492. if (iDet == kZDC) nPointsZDC += fPointsMap[a]; //NG save points in ZDC
  493. else nPointsNoZDC += fPointsMap[a]; // points not in ZDC
  494. }
  495. }
  496. // --> Check for cuts (store primaries in any case)
  497. if (iMother < 0) store = kTRUE;
  498. else {
  499. if (!fStoreSecondaries) store = kFALSE;
  500. if (nPoints < fMinPoints) store = kFALSE;
  501. if (eKin < fEnergyCut) store = kFALSE;
  502. if (fNoZDC && nPointsZDC && !nPointsNoZDC) store = kFALSE; //NG cut many secondary zdc tracks
  503. // !!!!!AZ - store products of potentially interesting decays
  504. if (nPoints == 0) {
  505. store = kFALSE;
  506. if (fStoreMothers) {
  507. TParticle* mother = GetParticle(iMother);
  508. if (mass < 0.1 && thisPart->P() < 0.01) store = kFALSE; // P < cut on e+- and gamma
  509. else if (mother->GetPDG()->Mass() > fMinMotherMass && mother->GetPDG()->Mass() < fMaxMotherMass && //0.4 and 6.1
  510. TMath::Abs(mother->GetPDG()->Mass()-0.940) > 0.004 && mother->R() < fRadiusCut && abs(mother->Vz()) < fVzCut)
  511. store = kTRUE; // except nucleons and heavy fragments
  512. }
  513. } //no points
  514. }
  515. // --> Set storage flag
  516. fStoreMap[i] = store;
  517. }
  518. //AZ If flag is set, flag recursively mothers of selected tracks and their sisters
  519. if (fStoreMothers) {
  520. // Fill mother map
  521. std::multimap<Int_t,Int_t> moths;
  522. std::multimap<Int_t,Int_t>::iterator it;
  523. pair<std::multimap<Int_t,Int_t>::iterator,std::multimap<Int_t,Int_t>::iterator> ret;
  524. std::map<Int_t,Bool_t> copyMap = fStoreMap;
  525. for (Int_t i=0; i<fNParticles; i++) {
  526. Int_t iMother = GetParticle(i)->GetMother(0);
  527. if (iMother >= 0) moths.insert(pair<Int_t,Int_t>(iMother,i));
  528. }
  529. for (Int_t i=0; i<fNParticles; i++) {
  530. if (copyMap[i]) { //if particle is to be stored
  531. Int_t iMother = GetParticle(i)->GetMother(0);
  532. while (iMother >= 0) { //and it's mother is not primary
  533. TParticle* mother = GetParticle(iMother);
  534. if (mother->GetPDG()->Mass() < 0.1 && mother->P() < 0.01) break;
  535. //if (mother->GetPDG()->Mass() < 0.1 && mother->P() < 0.001) break;
  536. if (abs(mother->Vz()) > fVzCut || mother->R() > fRadiusCut) break; //NG and doesn't originate within Rcut/Zcut
  537. fStoreMap[iMother] = kTRUE;//store the mother
  538. ret = moths.equal_range(iMother);
  539. for (it = ret.first; it != ret.second; ++it) {
  540. TParticle* daught = GetParticle(it->second);
  541. if (abs(daught->Vz()) > fVzCut || daught->R() > fRadiusCut) continue; //NG and doesn't originate within Rcut/Zcut
  542. if (daught->GetPDG()->Mass() < 0.1 && daught->P() < 0.01) continue;
  543. //if (daught->GetPDG()->Mass() < 0.1 && daught->P() < 0.001) continue;
  544. fStoreMap[it->second] = kTRUE; // sister
  545. }
  546. iMother = mother->GetMother(0); //mother of the mother, loop again
  547. } // while (iMother >= 0)
  548. } // if (copyMap[i])
  549. } // for (Int_t i=0; i<fNParticles;
  550. // AZ - If mother is not stored, find the nearest stored ancestor
  551. for (Int_t i = 0; i < fNParticles; ++i) {
  552. if (fStoreMap[i]) {
  553. Int_t iMother = GetParticle(i)->GetMother(0);
  554. while (iMother >= 0 && fStoreMap[iMother] == kFALSE) {
  555. TParticle* mother = GetParticle(iMother);
  556. iMother = mother->GetMother(0); //mother of the mother, loop again
  557. }
  558. if (iMother >= 0) GetParticle(i)->SetFirstMother(iMother);
  559. }
  560. }
  561. } // if (fStoreMothers)
  562. //*/
  563. }
  564. // -------------------------------------------------------------------------
  565. ClassImp(MpdStack)