TUHKMgen/UHKM/HadronDecayer.cxx

   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 /////////////////////////////////////////////////////////////////////////////////////
  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
  19 #include "DatabasePDG.h"
  20 #include "ParticlePDG.h"
  21 #include "DecayChannel.h"
  22 #include "UKUtility.h"
  23 #include "Particle.h"
  24 #include "HYJET_COMMONS.h"
  25 #include "HadronDecayer.h"
  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) {
  34   //
  35   // return a random decay time according to the particle's width
  36   //
  37
  38   ParticlePDG *pDef = parent.Def();
  39   Double_t width = pDef->GetWidth(); //GeV
  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) {
  46     const Double_t slope =  parent.E() * 0.1973 / (pDef->GetMass() * width);
  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
  57 void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, const DatabasePDG* database) {
  58   //
  59   // perform the decay
  60   //
  61
  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
  76   Double_t pdgMass = pDef->GetMass();
  77   Int_t comprCodePyth=0;
  78   Float_t delta =0;
  79
  80   Bool_t success = kFALSE;
  81   Int_t iterations = 0;
  82   // Try to decay the particle
  83   while(!success) {
  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     }
  91
  92     // get a random mass using the Breit Wigner distribution
  93     Double_t bwMass = gRandom->BreitWigner(pdgMass, pDef->GetWidth());
  94     //      bwMass = pdgMass;
  95     // Try to cut the Breit Wigner tail of the particle using the cuts from pythia
  96     // The delta variable is obtained from pythia based on the specie
  97     int encoding =pDef->GetPDG();
  98     SERVICEEV.ipdg = encoding;
  99     mydelta_();
 100     comprCodePyth=SERVICEEV.KC;
 101     delta = SERVICEEV.delta;// PYDAT2.PMAS[KC][3];
 102
 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;
 107     }
 108
 109     //bad delta - an exception
 110     if(comprCodePyth==254){
 111       bwMass=pdgMass;
 112       delta=0.0;
 113     }
 114
 115     // K0 decay into K0s or K0l
 116     if(TMath::Abs(encoding)==311) {
 117       bwMass=pdgMass;
 118       delta=0.0;
 119     }
 120
 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
 122     if(comprCodePyth!=0 && delta>0 && (bwMass<pdgMass-delta || bwMass>pdgMass+delta)){
 123       iterations++;
 124       continue;
 125     }
 126
 127     // check how many decay channels are allowed with the generated mass
 128     Int_t nAllowedChannels = database->GetNAllowedChannels(pDef, bwMass);
 129     // if no decay channels are posible with this mass, then generate another BW mass
 130     if(nAllowedChannels==0) {
 131       iterations++;
 132             std::cout << "HadronDecayer::Decay() no decays allowed at this BW mass" << std::endl;
 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) {
 148       for(Int_t nChannel = 0; nChannel < nDecayChannel; nChannel++) {
 149         randValue -= pDef->GetDecayChannel(nChannel)->GetBranching();
 150         if(randValue <= 0. && database->IsChannelAllowed(pDef->GetDecayChannel(nChannel), bwMass)) {
 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());
 166     Double_t bwEnergy = TMath::Sqrt(parent.Mom().X()*parent.Mom().X() +
 167                                     parent.Mom().Y()*parent.Mom().Y() +
 168                                     parent.Mom().Z()*parent.Mom().Z() +
 169                                     bwMass*bwMass);
 170
 171     Int_t nb = (Int_t)parent.GetType(); //particle from jets
 172
 173     TLorentzVector MomparentBW(parent.Mom().X(), parent.Mom().Y(), parent.Mom().Z(), bwEnergy);
 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)));
 183       p1.Pos(parentBW.Pos());
 184       p1.Mom(parent.Mom());
 185       p1.SetLastMotherPdg(parentBW.Encoding());
 186       p1.SetLastMotherDecayCoor(parentBW.Pos());
 187       p1.SetLastMotherDecayMom(parentBW.Mom());
 188       p1.SetType(nb);
 189
 190       // add the daughter to the list of secondaries
 191       Int_t parentIndex = parent.GetIndex();
 192       Int_t p1Index = p1.SetIndex();                           // set the daughter index
 193       p1.SetMother(parentIndex);                               // set the mother index for this daughter
 194       parent.SetFirstDaughterIndex(p1Index);
 195       parent.SetLastDaughterIndex(p1Index);
 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)
 208       MomAntiMom(p1.Mom(), p1.TableMass(), p2.Mom(), p2.TableMass(), bwMass);
 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());
 221       //set to daughters the same type as has mother
 222       p1.SetType(nb);
 223       p2.SetType(nb);
 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
 233       if(deltaS>0.001) {
 234         iterations++;
 235         continue;
 236       }
 237       // push particles to the list of secondaries
 238       Int_t parentIndex = parent.GetIndex();
 239       p1.SetIndex();
 240       p2.SetIndex();
 241       p1.SetMother(parentIndex);
 242       p2.SetMother(parentIndex);
 243       parent.SetFirstDaughterIndex(p1.GetIndex());
 244       parent.SetLastDaughterIndex(p2.GetIndex());
 245       allocator.AddParticle(p1, output);
 246       allocator.AddParticle(p2, output);
 247       success = kTRUE;
 248
 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();
 264       Double_t deltaMass = bwMass - mass1 - mass2 - mass3;
 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(e*e + 2*e*mass1);
 274         sumPabs = pAbs1;
 275         maxPabs = sumPabs;
 276         // 2
 277         e = (1-rd1)*deltaMass;
 278         pAbs2 = TMath::Sqrt(e*e + 2*e*mass2);
 279
 280         if(pAbs2 > maxPabs)
 281           maxPabs = pAbs2;
 282
 283         sumPabs += pAbs2;
 284         // 3
 285         e = (rd1-rd2)*deltaMass;
 286         pAbs3 = TMath::Sqrt(e*e + 2*e*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(sinTheta*cosPhi, sinTheta*sinPhi, cosTheta, 0);
 301       mom1 *= pAbs1;
 302       // sample rest particle 3-momentum
 303       Double_t cosThetaN = (pAbs2*pAbs2 - pAbs3*pAbs3 - pAbs1*pAbs1)/(2*pAbs1*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(sinThetaN*cosPhiN*cosTheta*cosPhi - sinThetaN*sinPhiN*sinPhi + cosThetaN*sinTheta*cosPhi,
 310                       sinThetaN*cosPhiN*cosTheta*sinPhi + sinThetaN*sinPhiN*cosPhi + cosThetaN*sinTheta*sinPhi,
 311                       -sinThetaN*cosPhiN*sinTheta + cosThetaN*cosTheta,
 312                       0.);
 313
 314       mom3 *= pAbs3*mom3.P();
 315       mom2 = mom1;
 316       mom2 += mom3;
 317       mom2 *= -1.;
 318       // calculate energy
 319       mom1.SetE(TMath::Sqrt(mom1.P()*mom1.P() + mass1*mass1));
 320       mom2.SetE(TMath::Sqrt(mom2.P()*mom2.P() + mass2*mass2));
 321       mom3.SetE(TMath::Sqrt(mom3.P()*mom3.P() + mass3*mass3));
 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
 339       p1.SetType(nb);
 340       p2.SetType(nb);
 341       p3.SetType(nb);
 342
 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) {
 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;
 357
 358         iterations++;
 359         continue;
 360       }
 361
 362       Int_t parentIndex = parent.GetIndex();
 363       p1.SetIndex();
 364       p2.SetIndex();
 365       p3.SetIndex();
 366       p1.SetMother(parentIndex);
 367       p2.SetMother(parentIndex);
 368       p3.SetMother(parentIndex);
 369       parent.SetFirstDaughterIndex(p1.GetIndex());
 370       parent.SetLastDaughterIndex(p3.GetIndex());
 371       allocator.AddParticle(p1, output);
 372       allocator.AddParticle(p2, output);
 373       allocator.AddParticle(p3, output);
 374       success = kTRUE;
 375     }
 376   }
 377   return;
 378 }