[u/mrichter/AliRoot.git] / TUHKMgen / UHKM / HadronDecayer.cxx

///////////////////////////////////////////////////////////////////////////////////// 
//                                                                                 //
//  July 2008 BW mass is limited by "PYTHIA method", by I. Lokhtin and L. Malinina //
//                                                                                 //
//        Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna        //
//      amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru     //
//                           November. 2, 2005                                     //
//                                                                                 //
/////////////////////////////////////////////////////////////////////////////////////

#include <functional>
#include <algorithm>
#include <vector>
#include <iostream>
#include <TRandom.h>
#include <TError.h>
#include <TMath.h>

#include "DatabasePDG.h"
#include "ParticlePDG.h"
#include "DecayChannel.h"
#include "UKUtility.h"
#include "Particle.h"
#include "HYJET_COMMONS.h"
#include "HadronDecayer.h"

//calculates decay time in fm/c
//calculates 1,2 and 3 body decays

using std::cout;
using std::endl;

Double_t GetDecayTime(const Particle &parent, Double_t weakDecayLimit) {
  //
  // return a random decay time according to the particle's width
  //
  
  ParticlePDG *pDef = parent.Def(); 
  Double_t width = pDef->GetWidth(); //GeV

  // if particle is set to be stable then return 0
  if(pDef->GetStableStatus())
    return 0.;

  if(width > weakDecayLimit && weakDecayLimit>=0.0) {
    const Double_t slope =  parent.E() * 0.1973 / (pDef->GetMass() * width);
    return -slope * TMath::Log(gRandom->Rndm());//in fm/c
  }

  return 0.;
}


extern "C" void mydelta_();
extern SERVICEEVCommon SERVICEEV;

void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, const DatabasePDG* database) {
  //
  // perform the decay
  //
  
  // Get the PDG properties of the particle
  ParticlePDG *pDef = parent.Def();

  // Get the number of posible decay channels
  Int_t nDecayChannel = pDef->GetNDecayChannels();

  // check the full branching of this specie
  Double_t fullBranching = pDef->GetFullBranching(); // Only 3 or less body decays

  // return if particle has no branching
  if(fullBranching < 0.00001)
    return;

  // get the PDG mass of the specie  
  Double_t pdgMass = pDef->GetMass();
  Int_t comprCodePyth=0;
  Float_t delta =0;

  Bool_t success = kFALSE;
  Int_t iterations = 0;
  // Try to decay the particle
  while(!success) {
    //    std::cout << "HadronDecayer::Decay() iteration #" << iterations << std::endl;
    if(iterations>1000) {
      std::cout << "HadronDecayer::Decay() more than 1000 iterations to decay particle with code " 
		<< pDef->GetPDG() << std::endl;
      std::cout << "               Will be left undecayed ... check it out!" << std::endl;
      return;
    }
    
    // get a random mass using the Breit Wigner distribution 
    Double_t bwMass = gRandom->BreitWigner(pdgMass, pDef->GetWidth());     
    //      bwMass = pdgMass;
    // Try to cut the Breit Wigner tail of the particle using the cuts from pythia
    // The delta variable is obtained from pythia based on the specie
    int encoding =pDef->GetPDG();
    SERVICEEV.ipdg = encoding;
    mydelta_();      
    comprCodePyth=SERVICEEV.KC;
    delta = SERVICEEV.delta;// PYDAT2.PMAS[KC][3];

    //if there are no such particle in PYTHIA particle table, we take delta=0.4
    if(comprCodePyth==0){
      bwMass=pdgMass; 
      delta=0.0;
    } 

    //bad delta - an exception
    if(comprCodePyth==254){
      bwMass=pdgMass; 
      delta=0.0;
    } 
      
    // K0 decay into K0s or K0l
    if(TMath::Abs(encoding)==311) {
      bwMass=pdgMass;
      delta=0.0;
    }

    //for particles from PYTHIA table only, if the BW mass is outside the cut range then quit this iteration and generate another BW mass
    if(comprCodePyth!=0 && delta>0 && (bwMass<pdgMass-delta || bwMass>pdgMass+delta)){
      iterations++;
      continue;
    }    
    
    // check how many decay channels are allowed with the generated mass
    Int_t nAllowedChannels = database->GetNAllowedChannels(pDef, bwMass);
    // if no decay channels are posible with this mass, then generate another BW mass
    if(nAllowedChannels==0) {    
      iterations++;
            std::cout << "HadronDecayer::Decay() no decays allowed at this BW mass" << std::endl;
      continue;
    }

    std::vector<Particle> apDaughter;
    std::vector<Double_t> dMass; //daughters'mass
    std::vector<Double_t> dMom;
    std::vector<Double_t> sm;
    std::vector<Double_t> rd;

    // we need to choose an allowed decay channel
    Double_t randValue = gRandom->Rndm() * fullBranching;
    Int_t chosenChannel = 1000;
    Bool_t found = kFALSE;
    Int_t channelIterations = 0;
    while(!found) {
      for(Int_t nChannel = 0; nChannel < nDecayChannel; nChannel++) {
	randValue -= pDef->GetDecayChannel(nChannel)->GetBranching();
	if(randValue <= 0. && database->IsChannelAllowed(pDef->GetDecayChannel(nChannel), bwMass)) {
	  chosenChannel = nChannel;
	  found = kTRUE;
	  break;
	}
      }
      channelIterations++;
    }

    // get the PDG information for the chosen decay channel
    DecayChannel *dc = pDef->GetDecayChannel(chosenChannel);
    Int_t nSec = dc->GetNDaughters();

    // Adjust the parent momentum four-vector for the MC generated Breit-Wigner mass
    Particle parentBW(database->GetPDGParticle(parent.Encoding()));
    parentBW.Pos(parent.Pos());
    Double_t bwEnergy = TMath::Sqrt(parent.Mom().X()*parent.Mom().X() + 
				    parent.Mom().Y()*parent.Mom().Y() +
				    parent.Mom().Z()*parent.Mom().Z() +
				    bwMass*bwMass);

    Int_t nb = (Int_t)parent.GetType(); //particle from jets

    TLorentzVector MomparentBW(parent.Mom().X(), parent.Mom().Y(), parent.Mom().Z(), bwEnergy); 
    parentBW.Mom(MomparentBW);
    // take into account BW when calculating boost velocity (for wide resonances it matters)
    TVector3 velocityBW(parentBW.Mom().BoostVector());

    // now we have an allowed decay
    // first case: one daughter particle
    if(nSec == 1) {
      // initialize the daughter particle
      Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
      p1.Pos(parentBW.Pos());
      p1.Mom(parent.Mom());
      p1.SetLastMotherPdg(parentBW.Encoding());
      p1.SetLastMotherDecayCoor(parentBW.Pos());
      p1.SetLastMotherDecayMom(parentBW.Mom());
      p1.SetType(nb);

      // add the daughter to the list of secondaries
      Int_t parentIndex = parent.GetIndex();
      Int_t p1Index = p1.SetIndex();                           // set the daughter index
      p1.SetMother(parentIndex);                               // set the mother index for this daughter 
      parent.SetFirstDaughterIndex(p1Index);
      parent.SetLastDaughterIndex(p1Index);
      allocator.AddParticle(p1, output);
      success = kTRUE;  
    }
    // second case: two daughter particles
    else if(nSec == 2) {
      // initialize the daughter particles
      Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
      p1.Pos(parentBW.Pos());
      Particle p2(database->GetPDGParticle(dc->GetDaughterPDG(1)));
      p2.Pos(parentBW.Pos());
      
      // calculate the momenta in rest frame of mother for the two particles (theta and phi are isotropic)
      MomAntiMom(p1.Mom(), p1.TableMass(), p2.Mom(), p2.TableMass(), bwMass);
    
      // boost to the laboratory system (to the mother velocity)
      p1.Mom().Boost(velocityBW);
      p2.Mom().Boost(velocityBW);

      //store information about mother
      p1.SetLastMotherPdg(parentBW.Encoding());
      p1.SetLastMotherDecayCoor(parentBW.Pos());
      p1.SetLastMotherDecayMom(parentBW.Mom());
      p2.SetLastMotherPdg(parentBW.Encoding());
      p2.SetLastMotherDecayCoor(parentBW.Pos());
      p2.SetLastMotherDecayMom(parentBW.Mom());
      //set to daughters the same type as has mother
      p1.SetType(nb);
      p2.SetType(nb);


      // check the kinematics in the lab system
      Double_t deltaS = TMath::Sqrt((parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X())*(parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X())+
				    (parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y())*(parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y())+
				    (parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z())*(parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z())+
				    (parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E())*(parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()));
      // if deltaS is too big then repeat the kinematic procedure
 
      if(deltaS>0.001) {
	iterations++;
	continue;
      }
      // push particles to the list of secondaries
      Int_t parentIndex = parent.GetIndex();
      p1.SetIndex(); 
      p2.SetIndex();
      p1.SetMother(parentIndex); 
      p2.SetMother(parentIndex);
      parent.SetFirstDaughterIndex(p1.GetIndex());
      parent.SetLastDaughterIndex(p2.GetIndex());
      allocator.AddParticle(p1, output);
      allocator.AddParticle(p2, output);
      success = kTRUE;
    
    }

    // third case: three daughter particle
    else if(nSec == 3) {
      // initialize the daughter particle
      Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
      p1.Pos(parentBW.Pos());
      Particle p2(database->GetPDGParticle(dc->GetDaughterPDG(1)));
      p2.Pos(parentBW.Pos());
      Particle p3(database->GetPDGParticle(dc->GetDaughterPDG(2)));
      p3.Pos(parentBW.Pos());
      // calculate the momenta in the rest frame of the mother particle
      Double_t pAbs1 = 0., pAbs2 = 0., pAbs3 = 0., sumPabs = 0., maxPabs = 0.;
      Double_t mass1 = p1.TableMass(), mass2 = p2.TableMass(), mass3 = p3.TableMass();
      TLorentzVector &mom1 = p1.Mom(), &mom2 = p2.Mom(), &mom3 = p3.Mom(); 
      Double_t deltaMass = bwMass - mass1 - mass2 - mass3;

      do {
	Double_t rd1 = gRandom->Rndm();
	Double_t rd2 = gRandom->Rndm();
	if (rd2 > rd1)
	  std::swap(rd1, rd2);
	// 1
	Double_t e = rd2*deltaMass;
	pAbs1 = TMath::Sqrt(e*e + 2*e*mass1);
	sumPabs = pAbs1;
	maxPabs = sumPabs;
	// 2
	e = (1-rd1)*deltaMass;
	pAbs2 = TMath::Sqrt(e*e + 2*e*mass2);
	
	if(pAbs2 > maxPabs)
	  maxPabs = pAbs2;
	
	sumPabs += pAbs2;
	// 3
	e = (rd1-rd2)*deltaMass;
	pAbs3 = TMath::Sqrt(e*e + 2*e*mass3);
	
	if (pAbs3 > maxPabs)
	  maxPabs =  pAbs3;
	sumPabs  +=  pAbs3;
      } while(maxPabs > sumPabs - maxPabs);
      
      // isotropic sample first particle 3-momentum
      Double_t cosTheta = 2*(gRandom->Rndm()) - 1;
      Double_t sinTheta = TMath::Sqrt(1 - cosTheta*cosTheta);
      Double_t phi      = TMath::TwoPi()*(gRandom->Rndm());
      Double_t sinPhi   = TMath::Sin(phi);
      Double_t cosPhi   = TMath::Cos(phi);
      
      mom1.SetPxPyPzE(sinTheta*cosPhi, sinTheta*sinPhi, cosTheta, 0);
      mom1 *= pAbs1;
      // sample rest particle 3-momentum
      Double_t cosThetaN = (pAbs2*pAbs2 - pAbs3*pAbs3 - pAbs1*pAbs1)/(2*pAbs1*pAbs3);
      Double_t sinThetaN = TMath::Sqrt(1 - cosThetaN*cosThetaN);
      Double_t phiN      = TMath::TwoPi()*(gRandom->Rndm());
      Double_t sinPhiN   = TMath::Sin(phiN);
      Double_t cosPhiN   = TMath::Cos(phiN);
      
      mom3.SetPxPyPzE(sinThetaN*cosPhiN*cosTheta*cosPhi - sinThetaN*sinPhiN*sinPhi + cosThetaN*sinTheta*cosPhi,
		      sinThetaN*cosPhiN*cosTheta*sinPhi + sinThetaN*sinPhiN*cosPhi + cosThetaN*sinTheta*sinPhi,
		      -sinThetaN*cosPhiN*sinTheta + cosThetaN*cosTheta,
		      0.);
      
      mom3 *= pAbs3*mom3.P();
      mom2 = mom1;
      mom2 += mom3;
      mom2 *= -1.;
      // calculate energy
      mom1.SetE(TMath::Sqrt(mom1.P()*mom1.P() + mass1*mass1));
      mom2.SetE(TMath::Sqrt(mom2.P()*mom2.P() + mass2*mass2));
      mom3.SetE(TMath::Sqrt(mom3.P()*mom3.P() + mass3*mass3));
      
      // boost to Lab system
      mom1.Boost(velocityBW);
      mom2.Boost(velocityBW);
      mom3.Boost(velocityBW);
      
      p1.SetLastMotherPdg(parentBW.Encoding());
      p1.SetLastMotherDecayCoor(parentBW.Pos());
      p1.SetLastMotherDecayMom(parentBW.Mom());
      p2.SetLastMotherPdg(parentBW.Encoding());
      p2.SetLastMotherDecayCoor(parentBW.Pos());
      p2.SetLastMotherDecayMom(parentBW.Mom());
      p3.SetLastMotherPdg(parentBW.Encoding());
      p3.SetLastMotherDecayCoor(parentBW.Pos());
      p3.SetLastMotherDecayMom(parentBW.Mom());

      //set to daughters the same type as has mother  
      p1.SetType(nb);
      p2.SetType(nb);
      p3.SetType(nb);

      // energy conservation check in the lab system
      Double_t deltaS = TMath::Sqrt((parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X()-p3.Mom().X())*(parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X()-p3.Mom().X()) +
				    (parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y()-p3.Mom().Y())*(parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y()-p3.Mom().Y()) +
				    (parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z()-p3.Mom().Z())*(parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z()-p3.Mom().Z())	+
				    (parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()-p3.Mom().E())*(parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()-p3.Mom().E()));
      // if deltaS is too big then repeat the kinematic procedure
      if(deltaS>0.001) {
	//	cout << "3-body decay kinematic check in lab system: " << pDef->GetPDG() << " --> " << p1.Encoding() << " + " << p2.Encoding() << " + " << p3.Encoding() << endl;
	//	cout << "Mother    (e,px,py,pz): " << parentBW.Mom().E() << "\t" << parentBW.Mom().X() << "\t" << parentBW.Mom().Y() << "\t" << parentBW.Mom().Z() << endl;
	//	cout << "Daughter1 (e,px,py,pz): " << p1.Mom().E() << "\t" << p1.Mom().X() << "\t" << p1.Mom().Y() << "\t" << p1.Mom().Z() << endl;
	//	cout << "Daughter2 (e,px,py,pz): " << p2.Mom().E() << "\t" << p2.Mom().X() << "\t" << p2.Mom().Y() << "\t" << p2.Mom().Z() << endl;
	//	cout << "Daughter3 (e,px,py,pz): " << p3.Mom().E() << "\t" << p3.Mom().X() << "\t" << p3.Mom().Y() << "\t" << p3.Mom().Z() << endl;
	//	cout << "3-body decay delta(sqrtS) = " << deltaS << endl;
	//	cout << "Repeating the decay algorithm..." << endl;

	iterations++;
	continue;
      }

      Int_t parentIndex = parent.GetIndex();
      p1.SetIndex();
      p2.SetIndex();
      p3.SetIndex();
      p1.SetMother(parentIndex); 
      p2.SetMother(parentIndex);
      p3.SetMother(parentIndex);
      parent.SetFirstDaughterIndex(p1.GetIndex());
      parent.SetLastDaughterIndex(p3.GetIndex());
      allocator.AddParticle(p1, output);
      allocator.AddParticle(p2, output);
      allocator.AddParticle(p3, output);
      success = kTRUE;
    }
  }
  return;
}
Commit	Line	Data
03896fc4	1	/////////////////////////////////////////////////////////////////////////////////////
	2	// //
	3	// July 2008 BW mass is limited by "PYTHIA method", by I. Lokhtin and L. Malinina //
	4	// //
	5	// Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna //
	6	// amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru //
	7	// November. 2, 2005 //
	8	// //
	9	/////////////////////////////////////////////////////////////////////////////////////
b1c2e580	10
	11	#include <functional>
	12	#include <algorithm>
	13	#include <vector>
	14	#include <iostream>
	15	#include <TRandom.h>
	16	#include <TError.h>
	17	#include <TMath.h>
	18
b1c2e580	19	#include "DatabasePDG.h"
b1c2e580	20	#include "ParticlePDG.h"
b1c2e580	21	#include "DecayChannel.h"
b1c2e580	22	#include "UKUtility.h"
b1c2e580	23	#include "Particle.h"
b1c2e580	24	#include "HYJET_COMMONS.h"
03896fc4	25	#include "HadronDecayer.h"
b1c2e580	26
	27	//calculates decay time in fm/c
	28	//calculates 1,2 and 3 body decays
	29
	30	using std::cout;
	31	using std::endl;
	32
	33	Double_t GetDecayTime(const Particle &parent, Double_t weakDecayLimit) {
03896fc4	34	//
	35	// return a random decay time according to the particle's width
	36	//
	37
b1c2e580	38	ParticlePDG *pDef = parent.Def();
b1c2e580	39	Double_t width = pDef->GetWidth(); //GeV
7b7936e9	40
	41	// if particle is set to be stable then return 0
	42	if(pDef->GetStableStatus())
	43	return 0.;
	44
	45	if(width > weakDecayLimit && weakDecayLimit>=0.0) {
b1c2e580	46	const Double_t slope = parent.E() * 0.1973 / (pDef->GetMass() * width);
b1c2e580	47	return -slope * TMath::Log(gRandom->Rndm());//in fm/c
	48	}
	49
	50	return 0.;
	51	}
	52
	53
	54	extern "C" void mydelta_();
	55	extern SERVICEEVCommon SERVICEEV;
	56
03896fc4	57	void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, const DatabasePDG* database) {
	58	//
	59	// perform the decay
	60	//
	61
b1c2e580	62	// Get the PDG properties of the particle
	63	ParticlePDG *pDef = parent.Def();
	64
	65	// Get the number of posible decay channels
	66	Int_t nDecayChannel = pDef->GetNDecayChannels();
	67
	68	// check the full branching of this specie
	69	Double_t fullBranching = pDef->GetFullBranching(); // Only 3 or less body decays
	70
	71	// return if particle has no branching
	72	if(fullBranching < 0.00001)
	73	return;
	74
	75	// get the PDG mass of the specie
03896fc4	76	Double_t pdgMass = pDef->GetMass();
	77	Int_t comprCodePyth=0;
	78	Float_t delta =0;
b1c2e580	79
	80	Bool_t success = kFALSE;
	81	Int_t iterations = 0;
	82	// Try to decay the particle
	83	while(!success) {
7b7936e9	84	// std::cout << "HadronDecayer::Decay() iteration #" << iterations << std::endl;
	85	if(iterations>1000) {
	86	std::cout << "HadronDecayer::Decay() more than 1000 iterations to decay particle with code "
	87	<< pDef->GetPDG() << std::endl;
	88	std::cout << " Will be left undecayed ... check it out!" << std::endl;
	89	return;
	90	}
3fa37a65	91
b1c2e580	92	// get a random mass using the Breit Wigner distribution
03896fc4	93	Double_t bwMass = gRandom->BreitWigner(pdgMass, pDef->GetWidth());
03896fc4	94	// bwMass = pdgMass;
b1c2e580	95	// Try to cut the Breit Wigner tail of the particle using the cuts from pythia
03896fc4	96	// The delta variable is obtained from pythia based on the specie
b1c2e580	97	int encoding =pDef->GetPDG();
	98	SERVICEEV.ipdg = encoding;
	99	mydelta_();
03896fc4	100	comprCodePyth=SERVICEEV.KC;
03896fc4	101	delta = SERVICEEV.delta;// PYDAT2.PMAS[KC][3];
b1c2e580	102
03896fc4	103	//if there are no such particle in PYTHIA particle table, we take delta=0.4
	104	if(comprCodePyth==0){
	105	bwMass=pdgMass;
	106	delta=0.0;
b1c2e580	107	}
	108
	109	//bad delta - an exception
03896fc4	110	if(comprCodePyth==254){
	111	bwMass=pdgMass;
	112	delta=0.0;
b1c2e580	113	}
b1c2e580	114
7b7936e9	115	// K0 decay into K0s or K0l
7b7936e9	116	if(TMath::Abs(encoding)==311) {
03896fc4	117	bwMass=pdgMass;
03896fc4	118	delta=0.0;
7b7936e9	119	}
7b7936e9	120
b1c2e580	121	//for particles from PYTHIA table only, if the BW mass is outside the cut range then quit this iteration and generate another BW mass
03896fc4	122	if(comprCodePyth!=0 && delta>0 && (bwMass<pdgMass-delta \|\| bwMass>pdgMass+delta)){
7b7936e9	123	iterations++;
b1c2e580	124	continue;
b1c2e580	125	}
b1c2e580	126
b1c2e580	127	// check how many decay channels are allowed with the generated mass
03896fc4	128	Int_t nAllowedChannels = database->GetNAllowedChannels(pDef, bwMass);
b1c2e580	129	// if no decay channels are posible with this mass, then generate another BW mass
	130	if(nAllowedChannels==0) {
	131	iterations++;
3fa37a65	132	std::cout << "HadronDecayer::Decay() no decays allowed at this BW mass" << std::endl;
b1c2e580	133	continue;
	134	}
	135
	136	std::vector<Particle> apDaughter;
	137	std::vector<Double_t> dMass; //daughters'mass
	138	std::vector<Double_t> dMom;
	139	std::vector<Double_t> sm;
	140	std::vector<Double_t> rd;
	141
	142	// we need to choose an allowed decay channel
	143	Double_t randValue = gRandom->Rndm() * fullBranching;
	144	Int_t chosenChannel = 1000;
	145	Bool_t found = kFALSE;
	146	Int_t channelIterations = 0;
	147	while(!found) {
3fa37a65	148	for(Int_t nChannel = 0; nChannel < nDecayChannel; nChannel++) {
b1c2e580	149	randValue -= pDef->GetDecayChannel(nChannel)->GetBranching();
03896fc4	150	if(randValue <= 0. && database->IsChannelAllowed(pDef->GetDecayChannel(nChannel), bwMass)) {
b1c2e580	151	chosenChannel = nChannel;
	152	found = kTRUE;
	153	break;
	154	}
	155	}
	156	channelIterations++;
	157	}
	158
	159	// get the PDG information for the chosen decay channel
	160	DecayChannel *dc = pDef->GetDecayChannel(chosenChannel);
	161	Int_t nSec = dc->GetNDaughters();
	162
	163	// Adjust the parent momentum four-vector for the MC generated Breit-Wigner mass
	164	Particle parentBW(database->GetPDGParticle(parent.Encoding()));
	165	parentBW.Pos(parent.Pos());
03896fc4	166	Double_t bwEnergy = TMath::Sqrt(parent.Mom().X()*parent.Mom().X() +
b1c2e580	167	parent.Mom().Y()*parent.Mom().Y() +
b1c2e580	168	parent.Mom().Z()*parent.Mom().Z() +
03896fc4	169	bwMass*bwMass);
b1c2e580	170
03896fc4	171	Int_t nb = (Int_t)parent.GetType(); //particle from jets
b1c2e580	172
03896fc4	173	TLorentzVector MomparentBW(parent.Mom().X(), parent.Mom().Y(), parent.Mom().Z(), bwEnergy);
b1c2e580	174	parentBW.Mom(MomparentBW);
	175	// take into account BW when calculating boost velocity (for wide resonances it matters)
	176	TVector3 velocityBW(parentBW.Mom().BoostVector());
	177
	178	// now we have an allowed decay
	179	// first case: one daughter particle
	180	if(nSec == 1) {
	181	// initialize the daughter particle
	182	Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
7b7936e9	183	p1.Pos(parentBW.Pos());
7b7936e9	184	p1.Mom(parent.Mom());
b1c2e580	185	p1.SetLastMotherPdg(parentBW.Encoding());
	186	p1.SetLastMotherDecayCoor(parentBW.Pos());
	187	p1.SetLastMotherDecayMom(parentBW.Mom());
03896fc4	188	p1.SetType(nb);
b1c2e580	189
3fa37a65	190	// add the daughter to the list of secondaries
3fa37a65	191	Int_t parentIndex = parent.GetIndex();
b1c2e580	192	Int_t p1Index = p1.SetIndex(); // set the daughter index
b1c2e580	193	p1.SetMother(parentIndex); // set the mother index for this daughter
3fa37a65	194	parent.SetFirstDaughterIndex(p1Index);
3fa37a65	195	parent.SetLastDaughterIndex(p1Index);
b1c2e580	196	allocator.AddParticle(p1, output);
	197	success = kTRUE;
	198	}
	199	// second case: two daughter particles
	200	else if(nSec == 2) {
	201	// initialize the daughter particles
	202	Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
	203	p1.Pos(parentBW.Pos());
	204	Particle p2(database->GetPDGParticle(dc->GetDaughterPDG(1)));
	205	p2.Pos(parentBW.Pos());
	206
	207	// calculate the momenta in rest frame of mother for the two particles (theta and phi are isotropic)
03896fc4	208	MomAntiMom(p1.Mom(), p1.TableMass(), p2.Mom(), p2.TableMass(), bwMass);
b1c2e580	209
	210	// boost to the laboratory system (to the mother velocity)
	211	p1.Mom().Boost(velocityBW);
	212	p2.Mom().Boost(velocityBW);
	213
	214	//store information about mother
	215	p1.SetLastMotherPdg(parentBW.Encoding());
	216	p1.SetLastMotherDecayCoor(parentBW.Pos());
	217	p1.SetLastMotherDecayMom(parentBW.Mom());
	218	p2.SetLastMotherPdg(parentBW.Encoding());
	219	p2.SetLastMotherDecayCoor(parentBW.Pos());
	220	p2.SetLastMotherDecayMom(parentBW.Mom());
b1c2e580	221	//set to daughters the same type as has mother
03896fc4	222	p1.SetType(nb);
03896fc4	223	p2.SetType(nb);
b1c2e580	224
	225
	226	// check the kinematics in the lab system
	227	Double_t deltaS = TMath::Sqrt((parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X())*(parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X())+
	228	(parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y())*(parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y())+
	229	(parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z())*(parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z())+
	230	(parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E())*(parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()));
	231	// if deltaS is too big then repeat the kinematic procedure
	232
b1c2e580	233	if(deltaS>0.001) {
b1c2e580	234	iterations++;
	235	continue;
	236	}
	237	// push particles to the list of secondaries
3fa37a65	238	Int_t parentIndex = parent.GetIndex();
b1c2e580	239	p1.SetIndex();
	240	p2.SetIndex();
	241	p1.SetMother(parentIndex);
	242	p2.SetMother(parentIndex);
3fa37a65	243	parent.SetFirstDaughterIndex(p1.GetIndex());
3fa37a65	244	parent.SetLastDaughterIndex(p2.GetIndex());
b1c2e580	245	allocator.AddParticle(p1, output);
	246	allocator.AddParticle(p2, output);
	247	success = kTRUE;
3fa37a65	248
b1c2e580	249	}
	250
	251	// third case: three daughter particle
	252	else if(nSec == 3) {
	253	// initialize the daughter particle
	254	Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
	255	p1.Pos(parentBW.Pos());
	256	Particle p2(database->GetPDGParticle(dc->GetDaughterPDG(1)));
	257	p2.Pos(parentBW.Pos());
	258	Particle p3(database->GetPDGParticle(dc->GetDaughterPDG(2)));
	259	p3.Pos(parentBW.Pos());
	260	// calculate the momenta in the rest frame of the mother particle
	261	Double_t pAbs1 = 0., pAbs2 = 0., pAbs3 = 0., sumPabs = 0., maxPabs = 0.;
	262	Double_t mass1 = p1.TableMass(), mass2 = p2.TableMass(), mass3 = p3.TableMass();
	263	TLorentzVector &mom1 = p1.Mom(), &mom2 = p2.Mom(), &mom3 = p3.Mom();
03896fc4	264	Double_t deltaMass = bwMass - mass1 - mass2 - mass3;
b1c2e580	265
	266	do {
	267	Double_t rd1 = gRandom->Rndm();
	268	Double_t rd2 = gRandom->Rndm();
	269	if (rd2 > rd1)
	270	std::swap(rd1, rd2);
	271	// 1
	272	Double_t e = rd2*deltaMass;
	273	pAbs1 = TMath::Sqrt(ee + 2e*mass1);
	274	sumPabs = pAbs1;
	275	maxPabs = sumPabs;
	276	// 2
	277	e = (1-rd1)*deltaMass;
	278	pAbs2 = TMath::Sqrt(ee + 2e*mass2);
	279
	280	if(pAbs2 > maxPabs)
	281	maxPabs = pAbs2;
	282
	283	sumPabs += pAbs2;
	284	// 3
	285	e = (rd1-rd2)*deltaMass;
	286	pAbs3 = TMath::Sqrt(ee + 2e*mass3);
	287
	288	if (pAbs3 > maxPabs)
	289	maxPabs = pAbs3;
	290	sumPabs += pAbs3;
	291	} while(maxPabs > sumPabs - maxPabs);
	292
	293	// isotropic sample first particle 3-momentum
	294	Double_t cosTheta = 2*(gRandom->Rndm()) - 1;
	295	Double_t sinTheta = TMath::Sqrt(1 - cosTheta*cosTheta);
	296	Double_t phi = TMath::TwoPi()*(gRandom->Rndm());
	297	Double_t sinPhi = TMath::Sin(phi);
	298	Double_t cosPhi = TMath::Cos(phi);
	299
	300	mom1.SetPxPyPzE(sinThetacosPhi, sinThetasinPhi, cosTheta, 0);
	301	mom1 *= pAbs1;
	302	// sample rest particle 3-momentum
	303	Double_t cosThetaN = (pAbs2pAbs2 - pAbs3pAbs3 - pAbs1pAbs1)/(2pAbs1*pAbs3);
	304	Double_t sinThetaN = TMath::Sqrt(1 - cosThetaN*cosThetaN);
	305	Double_t phiN = TMath::TwoPi()*(gRandom->Rndm());
	306	Double_t sinPhiN = TMath::Sin(phiN);
	307	Double_t cosPhiN = TMath::Cos(phiN);
	308
	309	mom3.SetPxPyPzE(sinThetaNcosPhiNcosThetacosPhi - sinThetaNsinPhiNsinPhi + cosThetaNsinTheta*cosPhi,
	310	sinThetaNcosPhiNcosThetasinPhi + sinThetaNsinPhiNcosPhi + cosThetaNsinTheta*sinPhi,
	311	-sinThetaNcosPhiNsinTheta + cosThetaN*cosTheta,
	312	0.);
	313
	314	mom3 = pAbs3mom3.P();
	315	mom2 = mom1;
	316	mom2 += mom3;
	317	mom2 *= -1.;
	318	// calculate energy
	319	mom1.SetE(TMath::Sqrt(mom1.P()mom1.P() + mass1mass1));
	320	mom2.SetE(TMath::Sqrt(mom2.P()mom2.P() + mass2mass2));
	321	mom3.SetE(TMath::Sqrt(mom3.P()mom3.P() + mass3mass3));
	322
	323	// boost to Lab system
	324	mom1.Boost(velocityBW);
	325	mom2.Boost(velocityBW);
	326	mom3.Boost(velocityBW);
	327
	328	p1.SetLastMotherPdg(parentBW.Encoding());
329	p1.SetLastMotherDecayCoor(parentBW.Pos());
330	p1.SetLastMotherDecayMom(parentBW.Mom());
331	p2.SetLastMotherPdg(parentBW.Encoding());
332	p2.SetLastMotherDecayCoor(parentBW.Pos());
333	p2.SetLastMotherDecayMom(parentBW.Mom());
334	p3.SetLastMotherPdg(parentBW.Encoding());
335	p3.SetLastMotherDecayCoor(parentBW.Pos());
336	p3.SetLastMotherDecayMom(parentBW.Mom());
337
338	//set to daughters the same type as has mother
03896fc4	339	p1.SetType(nb);
	340	p2.SetType(nb);
	341	p3.SetType(nb);
b1c2e580	342
b1c2e580	343	// energy conservation check in the lab system
	344	Double_t deltaS = TMath::Sqrt((parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X()-p3.Mom().X())*(parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X()-p3.Mom().X()) +
	345	(parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y()-p3.Mom().Y())*(parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y()-p3.Mom().Y()) +
	346	(parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z()-p3.Mom().Z())*(parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z()-p3.Mom().Z()) +
	347	(parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()-p3.Mom().E())*(parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()-p3.Mom().E()));
	348	// if deltaS is too big then repeat the kinematic procedure
	349	if(deltaS>0.001) {
7b7936e9	350	// cout << "3-body decay kinematic check in lab system: " << pDef->GetPDG() << " --> " << p1.Encoding() << " + " << p2.Encoding() << " + " << p3.Encoding() << endl;
	351	// cout << "Mother (e,px,py,pz): " << parentBW.Mom().E() << "\t" << parentBW.Mom().X() << "\t" << parentBW.Mom().Y() << "\t" << parentBW.Mom().Z() << endl;
	352	// cout << "Daughter1 (e,px,py,pz): " << p1.Mom().E() << "\t" << p1.Mom().X() << "\t" << p1.Mom().Y() << "\t" << p1.Mom().Z() << endl;
	353	// cout << "Daughter2 (e,px,py,pz): " << p2.Mom().E() << "\t" << p2.Mom().X() << "\t" << p2.Mom().Y() << "\t" << p2.Mom().Z() << endl;
	354	// cout << "Daughter3 (e,px,py,pz): " << p3.Mom().E() << "\t" << p3.Mom().X() << "\t" << p3.Mom().Y() << "\t" << p3.Mom().Z() << endl;
	355	// cout << "3-body decay delta(sqrtS) = " << deltaS << endl;
	356	// cout << "Repeating the decay algorithm..." << endl;
b1c2e580	357
	358	iterations++;
	359	continue;
	360	}
	361
3fa37a65	362	Int_t parentIndex = parent.GetIndex();
b1c2e580	363	p1.SetIndex();
	364	p2.SetIndex();
	365	p3.SetIndex();
	366	p1.SetMother(parentIndex);
	367	p2.SetMother(parentIndex);
	368	p3.SetMother(parentIndex);
3fa37a65	369	parent.SetFirstDaughterIndex(p1.GetIndex());
3fa37a65	370	parent.SetLastDaughterIndex(p3.GetIndex());
b1c2e580	371	allocator.AddParticle(p1, output);
	372	allocator.AddParticle(p2, output);
	373	allocator.AddParticle(p3, output);
	374	success = kTRUE;
	375	}
	376	}
	377	return;
	378	}