MpdDecayer.cxx 23 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726
  1. // -------------------------------------------------------------------------
  2. // ----- MpdDecayer source file -----
  3. // ----- Created 22/04/2020 -----
  4. // ----- R. Akhat, A. Zinchenko -----
  5. // ----- External decayer for MPD -----
  6. // -------------------------------------------------------------------------
  7. #include "MpdDecayer.h"
  8. #include "TClonesArray.h"
  9. #include "TDatabasePDG.h"
  10. #include <TDecayChannel.h>
  11. #include "TF1.h"
  12. //#include "TGeant3.h"
  13. #include "TH1.h"
  14. #include <THashList.h>
  15. #include "TLorentzVector.h"
  16. #include "TParticle.h"
  17. #include <TParticlePDG.h>
  18. #include "TPDGCode.h"
  19. #include "TPythia6.h"
  20. #include "TVirtualMC.h"
  21. //#include "G4ParticleDefinition.hh"
  22. //#include "G4ParticleTable.hh"
  23. //#include "G4DecayTable.hh"
  24. #include <fstream>
  25. #include <iostream>
  26. #include <iomanip>
  27. using std::cout;
  28. using std::endl;
  29. using std::set;
  30. ClassImp(MpdDecayer);
  31. MpdDecayer* MpdDecayer::fgInstance = 0;
  32. ////////////////////////////////////////////////////////////////////////////////
  33. /// Get the singleton object.
  34. MpdDecayer* MpdDecayer::Instance()
  35. {
  36. if (!fgInstance) fgInstance = new MpdDecayer;
  37. return fgInstance;
  38. }
  39. ////////////////////////////////////////////////////////////////////////////////
  40. /// Constructor
  41. MpdDecayer::MpdDecayer()
  42. // : fDecay(kMaxDecay),
  43. : fDecay(501),
  44. fBraPart(501)
  45. {
  46. fBraPart.Reset(1);
  47. Init();
  48. }
  49. ////////////////////////////////////////////////////////////////////////////////
  50. /// Initialize the decayer
  51. void MpdDecayer::Init()
  52. {
  53. static Bool_t init = kFALSE;
  54. if (init) return;
  55. init = kTRUE;
  56. //ForceDecay();
  57. fParticles = new TClonesArray("TParticle");
  58. // Get branching of Lmabda to p + \pi-
  59. TPythia6 *pythia = TPythia6::Instance();
  60. Int_t ccLamb = pythia->Pycomp(3122); // compressed code for Lambda
  61. Int_t idcb = pythia->GetMDCY(ccLamb,2);
  62. Int_t idce = idcb + pythia->GetMDCY(ccLamb,3);;
  63. //for (Int_t idc = idcb; idc < idce; ++idc) cout << pythia->GetBRAT(idc) << " ";
  64. //cout << endl;
  65. fBranch = pythia->GetBRAT(idcb);
  66. // Define additional particles
  67. DefineParticles();
  68. // Since the pi0 lifetime is less than 1.E-15 sec, if pi0 is produced as a decay
  69. // product in pythia6, e.g. KL0 -> pi0 pi+ pi-, the pi0 will be immediately
  70. // decayed by pythia6 to 2 gammas, and the KL0 decay product list passed
  71. // back to the transport mechanism will be "gamma gamma pi+ pi-", i.e.
  72. // the pi0 will not be visible in the list of secondaries passed back to
  73. // the transport mechanism and will not be pushed to the stack for possible
  74. // storage to the stdhep output array.
  75. // To avoid this, the pi0 is set to stable in pythia6, and its decay
  76. // will be handled by Geant.
  77. pythia->SetMDCY(pythia->Pycomp(111),1,0);
  78. if (TString(gMC->GetName()).Contains("TGeant4")) return;
  79. /*
  80. TGeant3 *g3 = (TGeant3*) gMC;
  81. // Get decay modes
  82. Int_t ipart = g3->IdFromPDG(3122);
  83. Int_t jpart = g3->Gclink()->jpart;
  84. Int_t jpa = g3->Lq()[jpart-ipart];
  85. Int_t jpa1 = g3->Lq()[jpa-1];
  86. Int_t jpa2 = g3->Lq()[jpa-2];
  87. for (Int_t i = 1; i <= 6; ++i) {
  88. fBratio[i-1] = g3->Q()[jpa1+i];
  89. fMode[i-1] = g3->Iq()[jpa2+i];
  90. std::cout << i << " " << fBratio[i-1] << " " << fMode[i-1] << std::endl;
  91. }
  92. */
  93. }
  94. ////////////////////////////////////////////////////////////////////////////////
  95. /// Decay a particle of type IDPART (PDG code) and momentum P.
  96. void MpdDecayer::Decay(Int_t idpart, TLorentzVector* p)
  97. {
  98. // Decay function
  99. //{ cout << " !!! Out !!! " << endl; exit(0); }
  100. // Reset internal particle container
  101. fParticles->Delete();
  102. fSourceFlag = kPythia;
  103. fMother.SetXYZT(0,0,0,0);
  104. if (!p) return;
  105. TPythia6 *pythia = TPythia6::Instance();
  106. TParticle *part = gMC->GetStack()->GetCurrentTrack();
  107. //part->Print();
  108. if (fMothersPdg.find(idpart) == fMothersPdg.end()) {
  109. // Particle is not defined
  110. new ((*fParticles)[0]) TParticle(*part); // store mother particle
  111. fSourceFlag = kCustom;
  112. return;
  113. }
  114. TVector3 polar;
  115. part->GetPolarisation(polar);
  116. //polar.SetX(1.0); // just for test
  117. //polar.Print();
  118. //p->Print();
  119. //exit(0);
  120. if (TMath::Abs(polar.X()) < 0.0001 && TMath::Abs(polar.Z()) < 0.0001) {
  121. //if (polar.X() < -0.0001 && polar.Z() < 0.0001) {
  122. // No polarisation - check this !!!
  123. pythia->Py1ent(0, idpart, p->Energy(), p->Theta(), p->Phi());
  124. pythia->GetPrimaries();
  125. /*
  126. Int_t kcLamb = pythia->Pycomp(3122); // compressed code for Lambda
  127. cout << " ----------------- " << pythia->GetMDCY(kcLamb,1) << " " << pythia->GetMDCY(kcLamb,2) << " "
  128. << pythia->GetMDCY(kcLamb,3) << endl;
  129. Int_t idcb = pythia->GetMDCY(kcLamb,2);
  130. Int_t idce = idcb + pythia->GetMDCY(kcLamb,3);;
  131. for (Int_t idc = idcb; idc < idce; ++idc) cout << pythia->GetBRAT(idc) << " ";
  132. cout << endl;
  133. */
  134. } else {
  135. // Polarized lambda - simulate anysotropic decay only to p + \pi-
  136. if (gRandom->Rndm() < fBranch) {
  137. // Anysotropic decay
  138. //cout << " ----------------- " << fBranch << endl;
  139. new ((*fParticles)[0]) TParticle(*part); // store mother particle
  140. Gdecay (idpart, p);
  141. } else {
  142. // Force the other channels
  143. Int_t idcb = pythia->GetMDCY(pythia->Pycomp(3122),2);
  144. pythia->SetMDME(idcb,1,0); // disable decay to p + \pi-
  145. pythia->Py1ent(0, idpart, p->Energy(), p->Theta(), p->Phi());
  146. pythia->GetPrimaries();
  147. pythia->SetMDME(idcb,1,1); // enable decay to p + \pi-
  148. }
  149. }
  150. /*
  151. // Restore decay modes - to do the decay in Geant3
  152. TGeant3 *g3 = (TGeant3*) gMC;
  153. g3->Gsdk(g3->IdFromPDG(3122), fBratio, fMode);
  154. //g3decay_ ();
  155. g3decayaz_();
  156. g3->SetUserDecay(3122);
  157. */
  158. }
  159. ////////////////////////////////////////////////////////////////////////////////
  160. /// Get the decay products into the passed PARTICLES TClonesArray of
  161. /// TParticles
  162. Int_t MpdDecayer::ImportParticles(TClonesArray *particles)
  163. {
  164. Int_t npart = 0;
  165. if (fSourceFlag == kPythia) {
  166. fParticles->Delete();
  167. //npart = TPythia6::Instance()->ImportParticles(particles,"All");
  168. npart = TPythia6::Instance()->ImportParticles(fParticles,"All");
  169. TLorentzVector lvDaugh;
  170. for (Int_t j = 0; j < npart; ++j) {
  171. TParticle *part = (TParticle*) fParticles->UncheckedAt(j);
  172. part->ProductionVertex(lvDaugh);
  173. lvDaugh += fMother;
  174. part->SetProductionVertex(lvDaugh);
  175. new ((*particles)[j]) TParticle(*part);
  176. }
  177. } else {
  178. npart = fParticles->GetEntriesFast();
  179. for (Int_t j = 0; j < npart; ++j) new ((*particles)[j]) TParticle(*((TParticle*)fParticles->UncheckedAt(j)));
  180. }
  181. /*
  182. cout << " Number of particles: " << npart << endl;
  183. for (Int_t j = 0; j < npart; ++j) {
  184. TParticle *part = (TParticle*) particles->UncheckedAt(j);
  185. part->Print();
  186. cout << " " << part->GetFirstMother() << " " << part->GetSecondMother() << " "
  187. << part->GetFirstDaughter() << " " << part->GetLastDaughter() << " " << part->GetNDaughters()
  188. << " " << part->GetStatusCode() << endl;
  189. }
  190. */
  191. return npart;
  192. }
  193. ////////////////////////////////////////////////////////////////////////////////
  194. /// Force a particular decay type
  195. void MpdDecayer::SetForceDecay(Int_t type)
  196. {
  197. //if (type > kMaxDecay) {
  198. if (type > 501) {
  199. Warning("SetForceDecay", "Invalid decay mode: %d", type);
  200. return;
  201. }
  202. //fDecay = EDecayType(type);
  203. }
  204. ////////////////////////////////////////////////////////////////////////////////
  205. /// Force a particle decay mode
  206. void MpdDecayer::ForceDecay()
  207. {
  208. }
  209. ////////////////////////////////////////////////////////////////////////////////
  210. /// Get the partial branching ratio for a particle of type IPART (a
  211. /// PDG code).
  212. Float_t MpdDecayer::GetPartialBranchingRatio(Int_t ipart)
  213. {
  214. /*
  215. Int_t kc = TPythia6::Instance()->Pycomp(TMath::Abs(ipart));
  216. // return TPythia6::Instance()->GetBRAT(kc);
  217. return fBraPart[kc];
  218. */
  219. return -1;
  220. }
  221. ////////////////////////////////////////////////////////////////////////////////
  222. /// Get the life-time of a particle of type KF (a PDG code).
  223. Float_t MpdDecayer::GetLifetime(Int_t kf)
  224. {
  225. /*
  226. Int_t kc=TPythia6::Instance()->Pycomp(TMath::Abs(kf));
  227. return TPythia6::Instance()->GetPMAS(kc,4) * 3.3333e-12;
  228. */
  229. return -1;
  230. }
  231. ////////////////////////////////////////////////////////////////////////////////
  232. /// Define additional particles
  233. void MpdDecayer::DefineParticles()
  234. {
  235. SetDecayTableFile("pythia6.dat");
  236. std::ofstream fout(fDecayTableFile);
  237. // Lambda(1520)
  238. //PrintPDG(TDatabasePDG::Instance()->GetParticle(3124),fout);
  239. MakeDecayList(3124, fout);
  240. fout.close();
  241. ReadDecayTable();
  242. TPythia6 *pythia = TPythia6::Instance();
  243. Int_t ccLamb = pythia->Pycomp(2214); // compressed code for phi
  244. Int_t idcb = pythia->GetMDCY(ccLamb,2);
  245. Int_t idce = idcb + pythia->GetMDCY(ccLamb,3);;
  246. for (Int_t idc = idcb; idc < idce; ++idc) cout << pythia->GetBRAT(idc) << " ";
  247. cout << endl;
  248. for (Int_t idc = idcb; idc < idce; ++idc) cout << pythia->GetMDME(idc,2) << " ";
  249. cout << endl;
  250. }
  251. ////////////////////////////////////////////////////////////////////////////////
  252. /// Read in particle data from an ASCII file. The file name must
  253. /// previously have been set using the member function
  254. /// SetDecayTableFile.
  255. void MpdDecayer::ReadDecayTable()
  256. {
  257. if (fDecayTableFile.IsNull()) {
  258. Warning("ReadDecayTable", "No file set");
  259. return;
  260. }
  261. Int_t lun = 15;
  262. TPythia6::Instance()->OpenFortranFile(lun,
  263. const_cast<char*>(fDecayTableFile.Data()));
  264. TPythia6::Instance()->Pyupda(3,lun);
  265. TPythia6::Instance()->CloseFortranFile(lun);
  266. }
  267. // ===================================================================
  268. // BEGIN COMMENT
  269. //
  270. // It would be better if the particle and decay information could be
  271. // read from the current TDatabasePDG instance.
  272. //
  273. // However, it seems to me that some information is missing. In
  274. // particular
  275. //
  276. // - The broadning cut-off,
  277. // - Resonance width
  278. // - Color charge
  279. // - MWID (?)
  280. //
  281. // Further more, it's not clear to me at least, what all the
  282. // parameters Pythia needs are.
  283. //
  284. // Code like the below could be used to make a temporary file that
  285. // Pythia could then read in. Ofcourse, one could also manipulate
  286. // the data structures directly, but that's propably more dangerous.
  287. //
  288. #if 1
  289. void MpdDecayer::PrintPDG(TParticlePDG* pdg,std::ofstream &fout)
  290. {
  291. TParticlePDG* anti = pdg->AntiParticle();
  292. const char* antiName = (anti ? anti->GetName() : "");
  293. Int_t color = 0;
  294. switch (TMath::Abs(pdg->PdgCode())) {
  295. case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8: // Quarks
  296. color = 1; break;
  297. case 21: // Gluon
  298. color = 2; break;
  299. case 1103:
  300. case 2101: case 2103: case 2203:
  301. case 3101: case 3103: case 3201: case 3203: case 3303:
  302. case 4101: case 4103: case 4201: case 4203: case 4301: case 4303: case 4403:
  303. case 5101: case 5103: case 5201: case 5203: case 5301: case 5303: case 5401:
  304. case 5403: case 5503:
  305. // Quark combinations
  306. color = -1; break;
  307. case 1000001: case 1000002: case 1000003: case 1000004: case 1000005:
  308. case 1000006: // super symmetric partners to quars
  309. color = 1; break;
  310. case 1000021: // ~g
  311. color = 2; break;
  312. case 2000001: case 2000002: case 2000003: case 2000004: case 2000005:
  313. case 2000006: // R hadrons
  314. color = 1; break;
  315. case 3000331: case 3100021: case 3200111: case 3100113: case 3200113:
  316. case 3300113: case 3400113:
  317. // Technicolor
  318. color = 2; break;
  319. case 4000001: case 4000002:
  320. color = 1; break;
  321. case 9900443: case 9900441: case 9910441: case 9900553: case 9900551:
  322. case 9910551:
  323. color = 2; break;
  324. }
  325. //std::cout << std::right
  326. fout << std::right
  327. << " " << std::setw(9) << pdg->PdgCode()
  328. << " " << std::left << std::setw(16) << pdg->GetName()
  329. << " " << std::setw(16) << antiName
  330. << std::right
  331. << std::setw(3) << Int_t(pdg->Charge())
  332. << std::setw(3) << color
  333. << std::setw(3) << (anti ? 1 : 0)
  334. << std::fixed << std::setprecision(5)
  335. << std::setw(12) << pdg->Mass()
  336. << std::setw(12) << pdg->Width()
  337. //<< std::setw(12) << 0.0 // Broad
  338. << std::setw(12) << 4*pdg->Width() // Broad
  339. //<< std::setw(12) << 8*pdg->Width() // Broad
  340. << std::scientific
  341. //<< " " << std::setw(13) << pdg->Lifetime()
  342. << std::setw(13) << pdg->Lifetime()
  343. << std::setw(3) << 0 // MWID
  344. << std::setw(3) << !pdg->Stable()
  345. << std::endl;
  346. }
  347. // ===================================================================
  348. void MpdDecayer::MakeDecayList(Int_t pdgCode, std::ofstream &fout)
  349. {
  350. TDatabasePDG* pdgDB = TDatabasePDG::Instance();
  351. if (!pdgDB->ParticleList()) pdgDB->ReadPDGTable();
  352. const THashList* pdgs = pdgDB->ParticleList();
  353. TParticlePDG* pdg = 0;
  354. TIter nextPDG(pdgs);
  355. while ((pdg = static_cast<TParticlePDG*>(nextPDG()))) {
  356. // std::cout << "Processing " << pdg->GetName() << std::endl;
  357. if (pdg->PdgCode() != pdgCode) continue;
  358. PrintPDG(pdg,fout);
  359. TObjArray* decays = pdg->DecayList();
  360. TDecayChannel* decay = 0;
  361. TIter nextDecay(decays);
  362. while ((decay = static_cast<TDecayChannel*>(nextDecay()))) {
  363. // std::cout << "Processing decay number " << decay->Number() << std::endl;
  364. Int_t ndecay = decay->NDaughters();
  365. fout << std::right << std::setw(15) << 1
  366. << std::setw(5) << decay->MatrixElementCode()
  367. << std::fixed << std::setw(12) << std::setprecision(6) << decay->BranchingRatio();
  368. for (Int_t id = 0; id < 5; ++id) {
  369. if (id < ndecay) fout << std::setw(10) << decay->DaughterPdgCode(id);
  370. else fout << std::setw(10) << 0;
  371. }
  372. fout << "\n";
  373. // " %5d%5d%12.5f%10d%10d%10d%10d%10d\n"
  374. }
  375. }
  376. }
  377. #endif
  378. // END COMMENT
  379. // ===================================================================
  380. ////////////////////////////////////////////////////////////////////////////////
  381. /// write particle data to an ASCII file. The file name must
  382. /// previously have been set using the member function
  383. /// SetDecayTableFile.
  384. ///
  385. /// Users can use this function to make an initial decay list file,
  386. /// which then can be edited by hand, and re-loaded into the decayer
  387. /// using ReadDecayTable.
  388. ///
  389. /// The file syntax is
  390. ///
  391. /// particle_list : partcle_data
  392. /// | particle_list particle_data
  393. /// ;
  394. /// particle_data : particle_info
  395. /// | particle_info '\n' decay_list
  396. /// ;
  397. /// particle_info : See below
  398. /// ;
  399. /// decay_list : decay_entry
  400. /// | decay_list decay_entry
  401. /// ;
  402. /// decay_entry : See below
  403. ///
  404. /// The particle_info consists of 13 fields:
  405. ///
  406. /// PDG code int
  407. /// Name string
  408. /// Anti-particle name string if there's no anti-particle,
  409. /// then this field must be the
  410. /// empty string
  411. /// Electic charge int in units of |e|/3
  412. /// Color charge int in units of quark color charges
  413. /// Have anti-particle int 1 of there's an anti-particle
  414. /// to this particle, or 0
  415. /// otherwise
  416. /// Mass float in units of GeV
  417. /// Resonance width float
  418. /// Max broadning float
  419. /// Lifetime float
  420. /// MWID int ??? (some sort of flag)
  421. /// Decay int 1 if it decays. 0 otherwise
  422. ///
  423. /// The format to write these entries in are
  424. ///
  425. /// " %9 %-16s %-16s%3d%3d%3d%12.5f%12.5f%12.5f%13.gf%3d%d\n"
  426. ///
  427. /// The decay_entry consists of 8 fields:
  428. ///
  429. /// On/Off int 1 for on, -1 for off
  430. /// Matrix element type int
  431. /// Branching ratio float
  432. /// Product 1 int PDG code of decay product 1
  433. /// Product 2 int PDG code of decay product 2
  434. /// Product 3 int PDG code of decay product 3
  435. /// Product 4 int PDG code of decay product 4
  436. /// Product 5 int PDG code of decay product 5
  437. ///
  438. /// The format for these lines are
  439. ///
  440. /// " %5d%5d%12.5f%10d%10d%10d%10d%10d\n"
  441. ///
  442. void MpdDecayer::WriteDecayTable()
  443. {
  444. if (fDecayTableFile.IsNull()) {
  445. Warning("ReadDecayTable", "No file set");
  446. return;
  447. }
  448. Int_t lun = 15;
  449. TPythia6::Instance()->OpenFortranFile(lun,
  450. const_cast<char*>(fDecayTableFile.Data()));
  451. TPythia6::Instance()->Pyupda(1,lun);
  452. TPythia6::Instance()->CloseFortranFile(lun);
  453. }
  454. ////////////////////////////////////////////////////////////////////////////////
  455. /// Count number of decay products
  456. Int_t MpdDecayer::CountProducts(Int_t channel, Int_t particle)
  457. {
  458. /*
  459. Int_t np = 0;
  460. for (Int_t i = 1; i <= 5; i++)
  461. if (TMath::Abs(TPythia6::Instance()->GetKFDP(channel,i)) == particle) np++;
  462. return np;
  463. */
  464. return -1;
  465. }
  466. ////////////////////////////////////////////////////////////////////////////////
  467. /// Simulate two-body decay process of Lambda-hyperon to p+\pi-
  468. void MpdDecayer::Gdecay(Int_t idpart, TLorentzVector* p)
  469. {
  470. Double_t pcm[2][4] = {{0},{0}};
  471. Double_t xm0 = p->M(), xm1 = TDatabasePDG::Instance()->GetParticle(2212)->Mass();
  472. Double_t xm2 = TDatabasePDG::Instance()->GetParticle(-211)->Mass();
  473. Gdeca2 (xm0, xm1, xm2, pcm);
  474. // First, second decay products.
  475. TLorentzVector daughter1, daughter2;
  476. daughter1.SetPxPyPzE (pcm[0][0], pcm[0][1], pcm[0][2], pcm[0][3]);
  477. daughter2.SetPxPyPzE (pcm[1][0], pcm[1][1], pcm[1][2], pcm[1][3]);
  478. //Boost to lab frame
  479. //TVector3 boost = p->BoostVector(); //AZ - this is wrong !!!!
  480. TVector3 boost (0, 0, p->P()/p->E());
  481. daughter1.Boost(boost);
  482. daughter2.Boost(boost);
  483. daughter1.Transform(fRotation);
  484. daughter2.Transform(fRotation);
  485. TLorentzVector pos;
  486. gMC->TrackPosition(pos);
  487. //pos.Print();
  488. Int_t npart = fParticles->GetEntriesFast();
  489. new ((*fParticles)[npart++]) TParticle(2212, 1, 1, -1, 0, 0, daughter1, pos);
  490. new ((*fParticles)[npart]) TParticle(-211, 1, 1, -1, 0, 0, daughter2, pos);
  491. fSourceFlag = kCustom;
  492. // Modify mother particle parameters
  493. TParticle *mother = (TParticle*) fParticles->UncheckedAt(0);
  494. mother->SetFirstMother(1);
  495. mother->SetLastMother(-1);
  496. mother->SetFirstDaughter(2);
  497. mother->SetLastDaughter(3);
  498. mother->SetStatusCode(11);
  499. }
  500. ////////////////////////////////////////////////////////////////////////////////
  501. /// Simulate two-body decay process with anisotropic angular distribution in CMS.
  502. void MpdDecayer::Gdeca2(Double_t xm0, Double_t xm1, Double_t xm2, Double_t pcm[2][4])
  503. {
  504. Double_t random[2], pvert[3], costh = 0.0, sinth = 0.0, phi = 0.0;
  505. Double_t e1 = (xm0 * xm0 + xm1 * xm1 - xm2 * xm2) / (2. * xm0);
  506. Double_t p1 = TMath::Sqrt (TMath::Abs ((e1 - xm1) * (e1 + xm1)));
  507. static TRandom *myRandom = NULL;
  508. if (!myRandom) { myRandom = new TRandom(); myRandom->SetSeed(0); }
  509. TRandom *oldRandom = gRandom;
  510. gRandom = myRandom;
  511. gRandom->RndmArray (2, random);
  512. // Sanity check - should not happen (legacy code)
  513. TParticle *part = gMC->GetStack()->GetCurrentTrack();
  514. if (part->GetPdgCode() != 3122) { cout << " ??? Not Lambda - exit " << part->GetPdgCode() << endl; exit(0); }
  515. TLorentzVector lorV;
  516. part->Momentum(lorV);
  517. lorV.Vect().GetXYZ(pvert);
  518. TVector3 polar;
  519. part->GetPolarisation(polar);
  520. if (TMath::Abs(polar.X()) < 0.0001 && TMath::Abs(polar.Z()) < 0.0001) {
  521. costh = 2. * random[0] - 1.0;
  522. phi = TMath::TwoPi() * random[1];
  523. } else {
  524. // Anisotropic decay angular distribution (0.5*(1+alpha*P*cos)).
  525. Anisotropy (pvert, random, polar.X(), phi, costh);
  526. }
  527. if (TMath::Abs(costh) >= 1.0) costh = TMath::Sign (1.0, costh);
  528. else sinth = TMath::Sqrt ((1. - costh) * (1. + costh));
  529. // Polar co-ordinates to momentum components.
  530. pcm[0][0] = p1 * sinth * TMath::Cos(phi);
  531. pcm[0][1] = p1 * sinth * TMath::Sin(phi);
  532. pcm[0][2] = p1 * costh;
  533. pcm[0][3] = e1;
  534. // Generate second decay product.
  535. pcm[1][0] = -pcm[0][0];
  536. pcm[1][1] = -pcm[0][1];
  537. pcm[1][2] = -pcm[0][2];
  538. pcm[1][3] = TMath::Sqrt (pcm[1][0] * pcm[1][0] + pcm[1][1] * pcm[1][1] + pcm[1][2] * pcm[1][2] + xm2 * xm2);
  539. gRandom = oldRandom;
  540. }
  541. ////////////////////////////////////////////////////////////////////////////////
  542. /// Simulate anisotropic (due to polarization) decay of lambda
  543. void MpdDecayer::Anisotropy (Double_t *pvert, Double_t *rndm, Double_t polar, Double_t &phi, Double_t &costh)
  544. {
  545. // Simulate anisotropic (due to polarization) decay of lambda
  546. //exit(0);
  547. //std::ofstream outfile ("costh.txt");
  548. //freopen ("costh.txt", "w", stdout);
  549. //const Double_t alpha = 0.642;
  550. static TF1 f("f","1+0.642*[0]*x",-1,1);
  551. f.SetParameter(0,polar);
  552. f.SetNpx(1000);
  553. Double_t costhe = f.GetRandom();
  554. Double_t sinth = 0.0;
  555. if (TMath::Abs(costhe) >= 1.0) costhe = TMath::Sign (1.0,costhe);
  556. else sinth = TMath::Sqrt (1.0 - costhe * costhe);
  557. phi = TMath::TwoPi() * rndm[1];
  558. //std::cout << " Cos: " << costhe << " " << phi << " " << pvert[0] << " " << pvert[1] << " " << pvert[2] << std::endl;
  559. // Compute normal vector
  560. TVector3 beam(0.0, 0.0, 1.0);
  561. TVector3 lambda(pvert[0], pvert[1], pvert[2]);
  562. TVector3 norm = beam.Cross(lambda).Unit();
  563. // Unit vector with theta and phi w.r.t. normal
  564. TVector3 unit(sinth*TMath::Cos(phi), sinth*TMath::Sin(phi), costhe);
  565. // Coordinate system transformation
  566. // (from lambda to lab)
  567. TVector3 lambU = lambda.Unit(), lambZ(0,0,1), lambX(1,0,0), lambY(0,1,0);
  568. lambZ.RotateUz(lambU);
  569. lambX.RotateUz(lambU);
  570. lambY.RotateUz(lambU);
  571. TRotation rotL;
  572. rotL.RotateAxes(lambX,lambY,lambZ); // transformation to lab. system
  573. fRotation = rotL;
  574. rotL.Invert(); // from lab. to lambda
  575. unit.RotateUz(norm); // to lab. system
  576. unit.Transform(rotL); // to lambda system
  577. /*
  578. TVector3 normZ(0,0,1), normX(1,0,0), normY(0,1,0);
  579. normZ.RotateUz(norm);
  580. normX.RotateUz(norm);
  581. normY.RotateUz(norm);
  582. TRotation rotN;
  583. rotN.RotateAxes(normX,normY,normZ); // transformation to lab. system
  584. unit.Transform(rotN); // to lab. system
  585. unit.Transform(rotL); // to lambda system
  586. */
  587. costh = TMath::Cos(unit.Theta());
  588. phi = unit.Phi();
  589. //std::cout << costh << " " << phi << std::endl;
  590. //outfile << i << " " << costh << endl;
  591. //outfile.close();
  592. }
  593. ///////////////////////////////////////////////////////////////////////////////
  594. /// Add PDG of particle to be decayed by this package
  595. void MpdDecayer::AddMotherPdg(Int_t pdg)
  596. {
  597. if (fMothersPdg.find(pdg) != fMothersPdg.end()) return;
  598. fMothersPdg.insert(pdg);
  599. /*
  600. G4DecayTable *decayTbl = NULL;
  601. if (TString(gMC->GetName()) == "TGeant4") {
  602. // Get decay table to feed it in again - bypass TGeant4 feature of zeroing decay table
  603. G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
  604. G4ParticleDefinition* particleDefinition = particleTable->FindParticle(pdg);
  605. decayTbl = particleDefinition->GetDecayTable();
  606. }
  607. */
  608. gMC->SetUserDecay(pdg);
  609. /*
  610. if (TString(gMC->GetName()) == "TGeant4") {
  611. G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
  612. G4ParticleDefinition* particleDefinition = particleTable->FindParticle(pdg);
  613. particleDefinition->SetDecayTable(decayTbl);
  614. }
  615. */
  616. }
  617. ////////////////////////////////////////////////////////////////////////////////
  618. /// Decay a TParticle
  619. void MpdDecayer::Decay(TParticle* p)
  620. {
  621. // Decay function
  622. //{ cout << " !!! Out !!! " << endl; exit(0); }
  623. // Reset internal particle container
  624. fParticles->Delete();
  625. fSourceFlag = kPythia;
  626. if (!p) return;
  627. p->ProductionVertex(fMother);
  628. TPythia6 *pythia = TPythia6::Instance();
  629. pythia->Py1ent(0, p->GetPdgCode(), p->Energy(), p->Theta(), p->Phi());
  630. pythia->GetPrimaries();
  631. }