123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726 |
- // -------------------------------------------------------------------------
- // ----- MpdDecayer source file -----
- // ----- Created 22/04/2020 -----
- // ----- R. Akhat, A. Zinchenko -----
- // ----- External decayer for MPD -----
- // -------------------------------------------------------------------------
- #include "MpdDecayer.h"
- #include "TClonesArray.h"
- #include "TDatabasePDG.h"
- #include <TDecayChannel.h>
- #include "TF1.h"
- //#include "TGeant3.h"
- #include "TH1.h"
- #include <THashList.h>
- #include "TLorentzVector.h"
- #include "TParticle.h"
- #include <TParticlePDG.h>
- #include "TPDGCode.h"
- #include "TPythia6.h"
- #include "TVirtualMC.h"
- //#include "G4ParticleDefinition.hh"
- //#include "G4ParticleTable.hh"
- //#include "G4DecayTable.hh"
- #include <fstream>
- #include <iostream>
- #include <iomanip>
- using std::cout;
- using std::endl;
- using std::set;
- ClassImp(MpdDecayer);
-
- MpdDecayer* MpdDecayer::fgInstance = 0;
-
- ////////////////////////////////////////////////////////////////////////////////
- /// Get the singleton object.
- MpdDecayer* MpdDecayer::Instance()
- {
- if (!fgInstance) fgInstance = new MpdDecayer;
- return fgInstance;
- }
- ////////////////////////////////////////////////////////////////////////////////
- /// Constructor
- MpdDecayer::MpdDecayer()
- // : fDecay(kMaxDecay),
- : fDecay(501),
- fBraPart(501)
- {
- fBraPart.Reset(1);
- Init();
- }
- ////////////////////////////////////////////////////////////////////////////////
- /// Initialize the decayer
- void MpdDecayer::Init()
- {
- static Bool_t init = kFALSE;
- if (init) return;
- init = kTRUE;
- //ForceDecay();
- fParticles = new TClonesArray("TParticle");
-
- // Get branching of Lmabda to p + \pi-
- TPythia6 *pythia = TPythia6::Instance();
- Int_t ccLamb = pythia->Pycomp(3122); // compressed code for Lambda
- Int_t idcb = pythia->GetMDCY(ccLamb,2);
- Int_t idce = idcb + pythia->GetMDCY(ccLamb,3);;
- //for (Int_t idc = idcb; idc < idce; ++idc) cout << pythia->GetBRAT(idc) << " ";
- //cout << endl;
- fBranch = pythia->GetBRAT(idcb);
-
- // Define additional particles
- DefineParticles();
- // Since the pi0 lifetime is less than 1.E-15 sec, if pi0 is produced as a decay
- // product in pythia6, e.g. KL0 -> pi0 pi+ pi-, the pi0 will be immediately
- // decayed by pythia6 to 2 gammas, and the KL0 decay product list passed
- // back to the transport mechanism will be "gamma gamma pi+ pi-", i.e.
- // the pi0 will not be visible in the list of secondaries passed back to
- // the transport mechanism and will not be pushed to the stack for possible
- // storage to the stdhep output array.
- // To avoid this, the pi0 is set to stable in pythia6, and its decay
- // will be handled by Geant.
- pythia->SetMDCY(pythia->Pycomp(111),1,0);
- if (TString(gMC->GetName()).Contains("TGeant4")) return;
- /*
- TGeant3 *g3 = (TGeant3*) gMC;
- // Get decay modes
- Int_t ipart = g3->IdFromPDG(3122);
- Int_t jpart = g3->Gclink()->jpart;
- Int_t jpa = g3->Lq()[jpart-ipart];
- Int_t jpa1 = g3->Lq()[jpa-1];
- Int_t jpa2 = g3->Lq()[jpa-2];
- for (Int_t i = 1; i <= 6; ++i) {
- fBratio[i-1] = g3->Q()[jpa1+i];
- fMode[i-1] = g3->Iq()[jpa2+i];
- std::cout << i << " " << fBratio[i-1] << " " << fMode[i-1] << std::endl;
- }
- */
- }
- ////////////////////////////////////////////////////////////////////////////////
- /// Decay a particle of type IDPART (PDG code) and momentum P.
- void MpdDecayer::Decay(Int_t idpart, TLorentzVector* p)
- {
- // Decay function
-
- //{ cout << " !!! Out !!! " << endl; exit(0); }
- // Reset internal particle container
- fParticles->Delete();
- fSourceFlag = kPythia;
- fMother.SetXYZT(0,0,0,0);
-
- if (!p) return;
- TPythia6 *pythia = TPythia6::Instance();
-
- TParticle *part = gMC->GetStack()->GetCurrentTrack();
- //part->Print();
- if (fMothersPdg.find(idpart) == fMothersPdg.end()) {
- // Particle is not defined
- new ((*fParticles)[0]) TParticle(*part); // store mother particle
- fSourceFlag = kCustom;
- return;
- }
- TVector3 polar;
- part->GetPolarisation(polar);
- //polar.SetX(1.0); // just for test
- //polar.Print();
- //p->Print();
- //exit(0);
- if (TMath::Abs(polar.X()) < 0.0001 && TMath::Abs(polar.Z()) < 0.0001) {
- //if (polar.X() < -0.0001 && polar.Z() < 0.0001) {
- // No polarisation - check this !!!
- pythia->Py1ent(0, idpart, p->Energy(), p->Theta(), p->Phi());
- pythia->GetPrimaries();
- /*
- Int_t kcLamb = pythia->Pycomp(3122); // compressed code for Lambda
- cout << " ----------------- " << pythia->GetMDCY(kcLamb,1) << " " << pythia->GetMDCY(kcLamb,2) << " "
- << pythia->GetMDCY(kcLamb,3) << endl;
- Int_t idcb = pythia->GetMDCY(kcLamb,2);
- Int_t idce = idcb + pythia->GetMDCY(kcLamb,3);;
- for (Int_t idc = idcb; idc < idce; ++idc) cout << pythia->GetBRAT(idc) << " ";
- cout << endl;
- */
- } else {
- // Polarized lambda - simulate anysotropic decay only to p + \pi-
- if (gRandom->Rndm() < fBranch) {
- // Anysotropic decay
- //cout << " ----------------- " << fBranch << endl;
- new ((*fParticles)[0]) TParticle(*part); // store mother particle
- Gdecay (idpart, p);
- } else {
- // Force the other channels
- Int_t idcb = pythia->GetMDCY(pythia->Pycomp(3122),2);
- pythia->SetMDME(idcb,1,0); // disable decay to p + \pi-
- pythia->Py1ent(0, idpart, p->Energy(), p->Theta(), p->Phi());
- pythia->GetPrimaries();
- pythia->SetMDME(idcb,1,1); // enable decay to p + \pi-
- }
- }
- /*
- // Restore decay modes - to do the decay in Geant3
- TGeant3 *g3 = (TGeant3*) gMC;
- g3->Gsdk(g3->IdFromPDG(3122), fBratio, fMode);
- //g3decay_ ();
- g3decayaz_();
- g3->SetUserDecay(3122);
- */
- }
- ////////////////////////////////////////////////////////////////////////////////
- /// Get the decay products into the passed PARTICLES TClonesArray of
- /// TParticles
- Int_t MpdDecayer::ImportParticles(TClonesArray *particles)
- {
- Int_t npart = 0;
-
- if (fSourceFlag == kPythia) {
- fParticles->Delete();
- //npart = TPythia6::Instance()->ImportParticles(particles,"All");
- npart = TPythia6::Instance()->ImportParticles(fParticles,"All");
- TLorentzVector lvDaugh;
-
- for (Int_t j = 0; j < npart; ++j) {
- TParticle *part = (TParticle*) fParticles->UncheckedAt(j);
- part->ProductionVertex(lvDaugh);
- lvDaugh += fMother;
- part->SetProductionVertex(lvDaugh);
- new ((*particles)[j]) TParticle(*part);
- }
- } else {
- npart = fParticles->GetEntriesFast();
- for (Int_t j = 0; j < npart; ++j) new ((*particles)[j]) TParticle(*((TParticle*)fParticles->UncheckedAt(j)));
- }
- /*
- cout << " Number of particles: " << npart << endl;
- for (Int_t j = 0; j < npart; ++j) {
- TParticle *part = (TParticle*) particles->UncheckedAt(j);
- part->Print();
- cout << " " << part->GetFirstMother() << " " << part->GetSecondMother() << " "
- << part->GetFirstDaughter() << " " << part->GetLastDaughter() << " " << part->GetNDaughters()
- << " " << part->GetStatusCode() << endl;
- }
- */
- return npart;
- }
- ////////////////////////////////////////////////////////////////////////////////
- /// Force a particular decay type
- void MpdDecayer::SetForceDecay(Int_t type)
- {
- //if (type > kMaxDecay) {
- if (type > 501) {
- Warning("SetForceDecay", "Invalid decay mode: %d", type);
- return;
- }
- //fDecay = EDecayType(type);
- }
- ////////////////////////////////////////////////////////////////////////////////
- /// Force a particle decay mode
- void MpdDecayer::ForceDecay()
- {
- }
- ////////////////////////////////////////////////////////////////////////////////
- /// Get the partial branching ratio for a particle of type IPART (a
- /// PDG code).
- Float_t MpdDecayer::GetPartialBranchingRatio(Int_t ipart)
- {
- /*
- Int_t kc = TPythia6::Instance()->Pycomp(TMath::Abs(ipart));
- // return TPythia6::Instance()->GetBRAT(kc);
- return fBraPart[kc];
- */
- return -1;
- }
- ////////////////////////////////////////////////////////////////////////////////
- /// Get the life-time of a particle of type KF (a PDG code).
- Float_t MpdDecayer::GetLifetime(Int_t kf)
- {
- /*
- Int_t kc=TPythia6::Instance()->Pycomp(TMath::Abs(kf));
- return TPythia6::Instance()->GetPMAS(kc,4) * 3.3333e-12;
- */
- return -1;
- }
- ////////////////////////////////////////////////////////////////////////////////
- /// Define additional particles
- void MpdDecayer::DefineParticles()
- {
- SetDecayTableFile("pythia6.dat");
- std::ofstream fout(fDecayTableFile);
- // Lambda(1520)
- //PrintPDG(TDatabasePDG::Instance()->GetParticle(3124),fout);
- MakeDecayList(3124, fout);
- fout.close();
- ReadDecayTable();
- TPythia6 *pythia = TPythia6::Instance();
- Int_t ccLamb = pythia->Pycomp(2214); // compressed code for phi
- Int_t idcb = pythia->GetMDCY(ccLamb,2);
- Int_t idce = idcb + pythia->GetMDCY(ccLamb,3);;
- for (Int_t idc = idcb; idc < idce; ++idc) cout << pythia->GetBRAT(idc) << " ";
- cout << endl;
- for (Int_t idc = idcb; idc < idce; ++idc) cout << pythia->GetMDME(idc,2) << " ";
- cout << endl;
- }
- ////////////////////////////////////////////////////////////////////////////////
- /// Read in particle data from an ASCII file. The file name must
- /// previously have been set using the member function
- /// SetDecayTableFile.
- void MpdDecayer::ReadDecayTable()
- {
- if (fDecayTableFile.IsNull()) {
- Warning("ReadDecayTable", "No file set");
- return;
- }
- Int_t lun = 15;
- TPythia6::Instance()->OpenFortranFile(lun,
- const_cast<char*>(fDecayTableFile.Data()));
- TPythia6::Instance()->Pyupda(3,lun);
- TPythia6::Instance()->CloseFortranFile(lun);
- }
- // ===================================================================
- // BEGIN COMMENT
- //
- // It would be better if the particle and decay information could be
- // read from the current TDatabasePDG instance.
- //
- // However, it seems to me that some information is missing. In
- // particular
- //
- // - The broadning cut-off,
- // - Resonance width
- // - Color charge
- // - MWID (?)
- //
- // Further more, it's not clear to me at least, what all the
- // parameters Pythia needs are.
- //
- // Code like the below could be used to make a temporary file that
- // Pythia could then read in. Ofcourse, one could also manipulate
- // the data structures directly, but that's propably more dangerous.
- //
- #if 1
- void MpdDecayer::PrintPDG(TParticlePDG* pdg,std::ofstream &fout)
- {
- TParticlePDG* anti = pdg->AntiParticle();
- const char* antiName = (anti ? anti->GetName() : "");
- Int_t color = 0;
- switch (TMath::Abs(pdg->PdgCode())) {
- case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8: // Quarks
- color = 1; break;
- case 21: // Gluon
- color = 2; break;
- case 1103:
- case 2101: case 2103: case 2203:
- case 3101: case 3103: case 3201: case 3203: case 3303:
- case 4101: case 4103: case 4201: case 4203: case 4301: case 4303: case 4403:
- case 5101: case 5103: case 5201: case 5203: case 5301: case 5303: case 5401:
- case 5403: case 5503:
- // Quark combinations
- color = -1; break;
- case 1000001: case 1000002: case 1000003: case 1000004: case 1000005:
- case 1000006: // super symmetric partners to quars
- color = 1; break;
- case 1000021: // ~g
- color = 2; break;
- case 2000001: case 2000002: case 2000003: case 2000004: case 2000005:
- case 2000006: // R hadrons
- color = 1; break;
- case 3000331: case 3100021: case 3200111: case 3100113: case 3200113:
- case 3300113: case 3400113:
- // Technicolor
- color = 2; break;
- case 4000001: case 4000002:
- color = 1; break;
- case 9900443: case 9900441: case 9910441: case 9900553: case 9900551:
- case 9910551:
- color = 2; break;
- }
- //std::cout << std::right
- fout << std::right
- << " " << std::setw(9) << pdg->PdgCode()
- << " " << std::left << std::setw(16) << pdg->GetName()
- << " " << std::setw(16) << antiName
- << std::right
- << std::setw(3) << Int_t(pdg->Charge())
- << std::setw(3) << color
- << std::setw(3) << (anti ? 1 : 0)
- << std::fixed << std::setprecision(5)
- << std::setw(12) << pdg->Mass()
- << std::setw(12) << pdg->Width()
- //<< std::setw(12) << 0.0 // Broad
- << std::setw(12) << 4*pdg->Width() // Broad
- //<< std::setw(12) << 8*pdg->Width() // Broad
- << std::scientific
- //<< " " << std::setw(13) << pdg->Lifetime()
- << std::setw(13) << pdg->Lifetime()
- << std::setw(3) << 0 // MWID
- << std::setw(3) << !pdg->Stable()
- << std::endl;
- }
- // ===================================================================
- void MpdDecayer::MakeDecayList(Int_t pdgCode, std::ofstream &fout)
- {
- TDatabasePDG* pdgDB = TDatabasePDG::Instance();
- if (!pdgDB->ParticleList()) pdgDB->ReadPDGTable();
- const THashList* pdgs = pdgDB->ParticleList();
- TParticlePDG* pdg = 0;
- TIter nextPDG(pdgs);
- while ((pdg = static_cast<TParticlePDG*>(nextPDG()))) {
- // std::cout << "Processing " << pdg->GetName() << std::endl;
- if (pdg->PdgCode() != pdgCode) continue;
- PrintPDG(pdg,fout);
-
- TObjArray* decays = pdg->DecayList();
- TDecayChannel* decay = 0;
- TIter nextDecay(decays);
- while ((decay = static_cast<TDecayChannel*>(nextDecay()))) {
- // std::cout << "Processing decay number " << decay->Number() << std::endl;
- Int_t ndecay = decay->NDaughters();
- fout << std::right << std::setw(15) << 1
- << std::setw(5) << decay->MatrixElementCode()
- << std::fixed << std::setw(12) << std::setprecision(6) << decay->BranchingRatio();
-
- for (Int_t id = 0; id < 5; ++id) {
- if (id < ndecay) fout << std::setw(10) << decay->DaughterPdgCode(id);
- else fout << std::setw(10) << 0;
- }
- fout << "\n";
- // " %5d%5d%12.5f%10d%10d%10d%10d%10d\n"
- }
- }
- }
- #endif
- // END COMMENT
- // ===================================================================
- ////////////////////////////////////////////////////////////////////////////////
- /// write particle data to an ASCII file. The file name must
- /// previously have been set using the member function
- /// SetDecayTableFile.
- ///
- /// Users can use this function to make an initial decay list file,
- /// which then can be edited by hand, and re-loaded into the decayer
- /// using ReadDecayTable.
- ///
- /// The file syntax is
- ///
- /// particle_list : partcle_data
- /// | particle_list particle_data
- /// ;
- /// particle_data : particle_info
- /// | particle_info '\n' decay_list
- /// ;
- /// particle_info : See below
- /// ;
- /// decay_list : decay_entry
- /// | decay_list decay_entry
- /// ;
- /// decay_entry : See below
- ///
- /// The particle_info consists of 13 fields:
- ///
- /// PDG code int
- /// Name string
- /// Anti-particle name string if there's no anti-particle,
- /// then this field must be the
- /// empty string
- /// Electic charge int in units of |e|/3
- /// Color charge int in units of quark color charges
- /// Have anti-particle int 1 of there's an anti-particle
- /// to this particle, or 0
- /// otherwise
- /// Mass float in units of GeV
- /// Resonance width float
- /// Max broadning float
- /// Lifetime float
- /// MWID int ??? (some sort of flag)
- /// Decay int 1 if it decays. 0 otherwise
- ///
- /// The format to write these entries in are
- ///
- /// " %9 %-16s %-16s%3d%3d%3d%12.5f%12.5f%12.5f%13.gf%3d%d\n"
- ///
- /// The decay_entry consists of 8 fields:
- ///
- /// On/Off int 1 for on, -1 for off
- /// Matrix element type int
- /// Branching ratio float
- /// Product 1 int PDG code of decay product 1
- /// Product 2 int PDG code of decay product 2
- /// Product 3 int PDG code of decay product 3
- /// Product 4 int PDG code of decay product 4
- /// Product 5 int PDG code of decay product 5
- ///
- /// The format for these lines are
- ///
- /// " %5d%5d%12.5f%10d%10d%10d%10d%10d\n"
- ///
- void MpdDecayer::WriteDecayTable()
- {
- if (fDecayTableFile.IsNull()) {
- Warning("ReadDecayTable", "No file set");
- return;
- }
- Int_t lun = 15;
- TPythia6::Instance()->OpenFortranFile(lun,
- const_cast<char*>(fDecayTableFile.Data()));
- TPythia6::Instance()->Pyupda(1,lun);
- TPythia6::Instance()->CloseFortranFile(lun);
- }
- ////////////////////////////////////////////////////////////////////////////////
- /// Count number of decay products
- Int_t MpdDecayer::CountProducts(Int_t channel, Int_t particle)
- {
- /*
- Int_t np = 0;
- for (Int_t i = 1; i <= 5; i++)
- if (TMath::Abs(TPythia6::Instance()->GetKFDP(channel,i)) == particle) np++;
- return np;
- */
- return -1;
- }
- ////////////////////////////////////////////////////////////////////////////////
- /// Simulate two-body decay process of Lambda-hyperon to p+\pi-
- void MpdDecayer::Gdecay(Int_t idpart, TLorentzVector* p)
- {
-
- Double_t pcm[2][4] = {{0},{0}};
- Double_t xm0 = p->M(), xm1 = TDatabasePDG::Instance()->GetParticle(2212)->Mass();
- Double_t xm2 = TDatabasePDG::Instance()->GetParticle(-211)->Mass();
-
- Gdeca2 (xm0, xm1, xm2, pcm);
- // First, second decay products.
-
- TLorentzVector daughter1, daughter2;
- daughter1.SetPxPyPzE (pcm[0][0], pcm[0][1], pcm[0][2], pcm[0][3]);
- daughter2.SetPxPyPzE (pcm[1][0], pcm[1][1], pcm[1][2], pcm[1][3]);
- //Boost to lab frame
- //TVector3 boost = p->BoostVector(); //AZ - this is wrong !!!!
- TVector3 boost (0, 0, p->P()/p->E());
- daughter1.Boost(boost);
- daughter2.Boost(boost);
- daughter1.Transform(fRotation);
- daughter2.Transform(fRotation);
- TLorentzVector pos;
- gMC->TrackPosition(pos);
- //pos.Print();
- Int_t npart = fParticles->GetEntriesFast();
- new ((*fParticles)[npart++]) TParticle(2212, 1, 1, -1, 0, 0, daughter1, pos);
- new ((*fParticles)[npart]) TParticle(-211, 1, 1, -1, 0, 0, daughter2, pos);
- fSourceFlag = kCustom;
- // Modify mother particle parameters
- TParticle *mother = (TParticle*) fParticles->UncheckedAt(0);
- mother->SetFirstMother(1);
- mother->SetLastMother(-1);
- mother->SetFirstDaughter(2);
- mother->SetLastDaughter(3);
- mother->SetStatusCode(11);
- }
- ////////////////////////////////////////////////////////////////////////////////
- /// Simulate two-body decay process with anisotropic angular distribution in CMS.
- void MpdDecayer::Gdeca2(Double_t xm0, Double_t xm1, Double_t xm2, Double_t pcm[2][4])
- {
- Double_t random[2], pvert[3], costh = 0.0, sinth = 0.0, phi = 0.0;
-
- Double_t e1 = (xm0 * xm0 + xm1 * xm1 - xm2 * xm2) / (2. * xm0);
- Double_t p1 = TMath::Sqrt (TMath::Abs ((e1 - xm1) * (e1 + xm1)));
- static TRandom *myRandom = NULL;
- if (!myRandom) { myRandom = new TRandom(); myRandom->SetSeed(0); }
- TRandom *oldRandom = gRandom;
- gRandom = myRandom;
- gRandom->RndmArray (2, random);
-
- // Sanity check - should not happen (legacy code)
- TParticle *part = gMC->GetStack()->GetCurrentTrack();
- if (part->GetPdgCode() != 3122) { cout << " ??? Not Lambda - exit " << part->GetPdgCode() << endl; exit(0); }
- TLorentzVector lorV;
- part->Momentum(lorV);
- lorV.Vect().GetXYZ(pvert);
-
- TVector3 polar;
- part->GetPolarisation(polar);
-
- if (TMath::Abs(polar.X()) < 0.0001 && TMath::Abs(polar.Z()) < 0.0001) {
- costh = 2. * random[0] - 1.0;
- phi = TMath::TwoPi() * random[1];
- } else {
- // Anisotropic decay angular distribution (0.5*(1+alpha*P*cos)).
- Anisotropy (pvert, random, polar.X(), phi, costh);
- }
- if (TMath::Abs(costh) >= 1.0) costh = TMath::Sign (1.0, costh);
- else sinth = TMath::Sqrt ((1. - costh) * (1. + costh));
-
- // Polar co-ordinates to momentum components.
- pcm[0][0] = p1 * sinth * TMath::Cos(phi);
- pcm[0][1] = p1 * sinth * TMath::Sin(phi);
- pcm[0][2] = p1 * costh;
- pcm[0][3] = e1;
- // Generate second decay product.
- pcm[1][0] = -pcm[0][0];
- pcm[1][1] = -pcm[0][1];
- pcm[1][2] = -pcm[0][2];
- 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);
- gRandom = oldRandom;
- }
- ////////////////////////////////////////////////////////////////////////////////
- /// Simulate anisotropic (due to polarization) decay of lambda
- void MpdDecayer::Anisotropy (Double_t *pvert, Double_t *rndm, Double_t polar, Double_t &phi, Double_t &costh)
- {
- // Simulate anisotropic (due to polarization) decay of lambda
- //exit(0);
- //std::ofstream outfile ("costh.txt");
- //freopen ("costh.txt", "w", stdout);
- //const Double_t alpha = 0.642;
- static TF1 f("f","1+0.642*[0]*x",-1,1);
- f.SetParameter(0,polar);
- f.SetNpx(1000);
- Double_t costhe = f.GetRandom();
- Double_t sinth = 0.0;
- if (TMath::Abs(costhe) >= 1.0) costhe = TMath::Sign (1.0,costhe);
- else sinth = TMath::Sqrt (1.0 - costhe * costhe);
- phi = TMath::TwoPi() * rndm[1];
- //std::cout << " Cos: " << costhe << " " << phi << " " << pvert[0] << " " << pvert[1] << " " << pvert[2] << std::endl;
- // Compute normal vector
- TVector3 beam(0.0, 0.0, 1.0);
- TVector3 lambda(pvert[0], pvert[1], pvert[2]);
- TVector3 norm = beam.Cross(lambda).Unit();
- // Unit vector with theta and phi w.r.t. normal
- TVector3 unit(sinth*TMath::Cos(phi), sinth*TMath::Sin(phi), costhe);
- // Coordinate system transformation
- // (from lambda to lab)
- TVector3 lambU = lambda.Unit(), lambZ(0,0,1), lambX(1,0,0), lambY(0,1,0);
- lambZ.RotateUz(lambU);
- lambX.RotateUz(lambU);
- lambY.RotateUz(lambU);
- TRotation rotL;
- rotL.RotateAxes(lambX,lambY,lambZ); // transformation to lab. system
- fRotation = rotL;
- rotL.Invert(); // from lab. to lambda
- unit.RotateUz(norm); // to lab. system
- unit.Transform(rotL); // to lambda system
- /*
- TVector3 normZ(0,0,1), normX(1,0,0), normY(0,1,0);
- normZ.RotateUz(norm);
- normX.RotateUz(norm);
- normY.RotateUz(norm);
- TRotation rotN;
- rotN.RotateAxes(normX,normY,normZ); // transformation to lab. system
- unit.Transform(rotN); // to lab. system
- unit.Transform(rotL); // to lambda system
- */
- costh = TMath::Cos(unit.Theta());
- phi = unit.Phi();
- //std::cout << costh << " " << phi << std::endl;
- //outfile << i << " " << costh << endl;
- //outfile.close();
- }
- ///////////////////////////////////////////////////////////////////////////////
- /// Add PDG of particle to be decayed by this package
- void MpdDecayer::AddMotherPdg(Int_t pdg)
- {
- if (fMothersPdg.find(pdg) != fMothersPdg.end()) return;
- fMothersPdg.insert(pdg);
- /*
- G4DecayTable *decayTbl = NULL;
- if (TString(gMC->GetName()) == "TGeant4") {
- // Get decay table to feed it in again - bypass TGeant4 feature of zeroing decay table
- G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
- G4ParticleDefinition* particleDefinition = particleTable->FindParticle(pdg);
- decayTbl = particleDefinition->GetDecayTable();
- }
- */
- gMC->SetUserDecay(pdg);
- /*
- if (TString(gMC->GetName()) == "TGeant4") {
- G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
- G4ParticleDefinition* particleDefinition = particleTable->FindParticle(pdg);
- particleDefinition->SetDecayTable(decayTbl);
- }
- */
- }
- ////////////////////////////////////////////////////////////////////////////////
- /// Decay a TParticle
- void MpdDecayer::Decay(TParticle* p)
- {
- // Decay function
-
- //{ cout << " !!! Out !!! " << endl; exit(0); }
- // Reset internal particle container
- fParticles->Delete();
- fSourceFlag = kPythia;
-
- if (!p) return;
- p->ProductionVertex(fMother);
- TPythia6 *pythia = TPythia6::Instance();
-
- pythia->Py1ent(0, p->GetPdgCode(), p->Energy(), p->Theta(), p->Phi());
- pythia->GetPrimaries();
- }
|