TUHKMgen/UHKM/HadronDecayer.cxx

   1 /*
   2
   3        July 2008 BW mass is limited by "PYTHIA method", by I. Lokhtin and L. Malinina
   4
   5
   6         Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
   7       amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru
   8                            November. 2, 2005
   9
  10 */
  11
  12 #include <functional>
  13 #include <algorithm>
  14 #include <vector>
  15 #include <iostream>
  16 #include <TRandom.h>
  17 #include <TError.h>
  18 #include <TMath.h>
  19
  20 #ifndef DATABASE_PDG
  21 #include "DatabasePDG.h"
  22 #endif
  23 #ifndef PARTICLE_PDG
  24 #include "ParticlePDG.h"
  25 #endif
  26 #ifndef DECAY_CHANNEL
  27 #include "DecayChannel.h"
  28 #endif
  29 #ifndef HADRONDECAYER_INCLUDED
  30 #include "HadronDecayer.h"
  31 #endif
  32 #ifndef UKUTILITY_INCLUDED
  33 #include "UKUtility.h"
  34 #endif
  35 #ifndef PARTICLE_INCLUDED
  36 #include "Particle.h"
  37 #endif
  38 #include "HYJET_COMMONS.h"
  39
  40 //calculates decay time in fm/c
  41 //calculates 1,2 and 3 body decays
  42
  43 using std::cout;
  44 using std::endl;
  45
  46 Double_t GetDecayTime(const Particle &parent, Double_t weakDecayLimit) {
  47   ParticlePDG *pDef = parent.Def();
  48   Double_t width = pDef->GetWidth(); //GeV
  49   if(width > weakDecayLimit) {
  50     const Double_t slope =  parent.E() * 0.1973 / (pDef->GetMass() * width);
  51 //    cout<<"get decay time"<<pDef->GetPDG()<<" "<<slope<<endl;
  52     return -slope * TMath::Log(gRandom->Rndm());//in fm/c
  53   }
  54
  55   return 0.;
  56 }
  57
  58
  59 extern "C" void mydelta_();
  60 extern SERVICEEVCommon SERVICEEV;
  61
  62 void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, DatabasePDG* database) {
  63   // check if the parent particle has been decayed already
  64
  65   Int_t daughters = parent.GetNDaughters();
  66
  67 //cout<<"in Decay pdg"<< parent.Def()->GetPDG()<<"daughters "<<daughters<<endl;
  68
  69   if(daughters>0)  // particle decayed already
  70     return;
  71
  72   // Get the PDG properties of the particle
  73   ParticlePDG *pDef = parent.Def();
  74
  75   // Get the number of posible decay channels
  76   Int_t nDecayChannel = pDef->GetNDecayChannels();
  77
  78   // check the full branching of this specie
  79   Double_t fullBranching = pDef->GetFullBranching(); // Only 3 or less body decays
  80
  81   // return if particle has no branching
  82   if(fullBranching < 0.00001)
  83     return;
  84
  85   // get the PDG mass of the specie
  86   Double_t PDGmass = pDef->GetMass();
  87   int ComprCodePyth=0;
  88   Float_t Delta =0;
  89
  90   Bool_t success = kFALSE;
  91   Int_t iterations = 0;
  92   // Try to decay the particle
  93   while(!success) {
  94     // get a random mass using the Breit Wigner distribution
  95     Double_t BWmass = gRandom->BreitWigner(PDGmass, pDef->GetWidth());
  96     //!!!!
  97     //      BWmass = PDGmass;
  98     // Try to cut the Breit Wigner tail of the particle using the cuts from pythia
  99     // The Delta variable is obtained from pythia based on the specie
 100     int encoding =pDef->GetPDG();
 101     SERVICEEV.ipdg = encoding;
 102     mydelta_();
 103     ComprCodePyth=SERVICEEV.KC;
 104     Delta = SERVICEEV.delta;// PYDAT2.PMAS[KC][3];
 105
 106     //if there are no such particle in PYTHIA particle table, we take Delta=0.4
 107     if(ComprCodePyth==0){
 108       BWmass=PDGmass;
 109       Delta=0.0;
 110     }
 111
 112     //bad delta - an exception
 113     if(ComprCodePyth==254){
 114       BWmass=PDGmass;
 115       Delta=0.0;
 116     }
 117
 118     //for particles from PYTHIA table only, if the BW mass is outside the cut range then quit this iteration and generate another BW mass
 119     if(ComprCodePyth!=0 && Delta>0 && (BWmass<PDGmass-Delta || BWmass>PDGmass+Delta)){
 120       //      std::cout<<"encoding"<<encoding<<"delta"<<Delta<<"width "<<pDef->GetWidth()<<"mass"<<BWmass<<std::endl;
 121       continue;
 122     }
 123     //----
 124
 125     if(BWmass>5)
 126       std::cout<<" > 5 encoding"<<encoding<<" pdgmass "<<PDGmass<<" delta "<<Delta<<"width "<<pDef->GetWidth()<<" mass "<<BWmass<<"CC"<<ComprCodePyth<<std::endl;
 127
 128     // check how many decay channels are allowed with the generated mass
 129     Int_t nAllowedChannels = database->GetNAllowedChannels(pDef, BWmass);
 130     // if no decay channels are posible with this mass, then generate another BW mass
 131     if(nAllowedChannels==0) {
 132       iterations++;
 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.SetLastMotherPdg(parentBW.Encoding());
 184       p1.SetLastMotherDecayCoor(parentBW.Pos());
 185       p1.SetLastMotherDecayMom(parentBW.Mom());
 186
 187       // link the parent and daughters trough their indexes in the list
 188       Int_t parentIndex = -1;
 189       if(parent.GetMother()==-1) parentIndex = parent.SetIndex();   // parents which are primaries don't have yet an index
 190       else parentIndex = parent.GetIndex();                         // parents which are secondaries have an index
 191       Int_t p1Index = p1.SetIndex();                           // set the daughter index
 192       p1.SetMother(parentIndex);                               // set the mother index for this daughter
 193       parent.SetDaughter(p1Index);                             // set p1 as daughter to the parent
 194       if(parent.GetMother()==-1) allocator.AddParticle(parent, output);   // add it only if its a primary particle
 195       allocator.AddParticle(p1, output);
 196       success = kTRUE;
 197     }
 198     // second case: two daughter particles
 199     else if(nSec == 2) {
 200       // initialize the daughter particles
 201       Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
 202       p1.Pos(parentBW.Pos());
 203       Particle p2(database->GetPDGParticle(dc->GetDaughterPDG(1)));
 204       p2.Pos(parentBW.Pos());
 205
 206       // calculate the momenta in rest frame of mother for the two particles (theta and phi are isotropic)
 207       MomAntiMom(p1.Mom(), p1.TableMass(), p2.Mom(), p2.TableMass(), BWmass);
 208
 209       // boost to the laboratory system (to the mother velocity)
 210       p1.Mom().Boost(velocityBW);
 211       p2.Mom().Boost(velocityBW);
 212
 213       //store information about mother
 214       p1.SetLastMotherPdg(parentBW.Encoding());
 215       p1.SetLastMotherDecayCoor(parentBW.Pos());
 216       p1.SetLastMotherDecayMom(parentBW.Mom());
 217       p2.SetLastMotherPdg(parentBW.Encoding());
 218       p2.SetLastMotherDecayCoor(parentBW.Pos());
 219       p2.SetLastMotherDecayMom(parentBW.Mom());
 220       // std::cout<<"2d NB="<<NB<<std::endl;
 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
 234       if(deltaS>0.001) {
 235
 236         cout << "2-body decay kinematic check in lab system: " << pDef->GetPDG() << " --> " << p1.Encoding() << " + " << p2.Encoding() << endl;
 237         cout << "Mother    (e,px,py,pz): " << parentBW.Mom().E() << "\t" << parentBW.Mom().X() << "\t" << parentBW.Mom().Y() << "\t" << parentBW.Mom().Z() << endl;
 238         cout << "Mother    (x,y,z,t): " << parentBW.Pos().X() << "\t" << parentBW.Pos().Y() << "\t" << parentBW.Pos().Z() << "\t" << parentBW.Pos().T() << endl;
 239
 240         cout << "Daughter1 (e,px,py,pz): " << p1.Mom().E() << "\t" << p1.Mom().X() << "\t" << p1.Mom().Y() << "\t" << p1.Mom().Z() << endl;
 241         cout << "Daughter2 (e,px,py,pz): " << p2.Mom().E() << "\t" << p2.Mom().X() << "\t" << p2.Mom().Y() << "\t" << p2.Mom().Z() << endl;
 242         cout << "2-body decay delta(sqrtS) = " << deltaS << endl;
 243         cout << "Repeating the decay algorithm ..." << endl;
 244
 245         iterations++;
 246         continue;
 247       }
 248       // push particles to the list of secondaries
 249       Int_t parentIndex = -1;
 250       if(parent.GetMother()==-1) parentIndex = parent.SetIndex();   // parents which are primaries don't have yet an index
 251       else parentIndex = parent.GetIndex();                         // parents which are secondaries have an index
 252       p1.SetIndex();
 253       p2.SetIndex();
 254       p1.SetMother(parentIndex);
 255       p2.SetMother(parentIndex);
 256       parent.SetDaughter(p1.GetIndex());
 257       parent.SetDaughter(p2.GetIndex());
 258       if(parent.GetMother()==-1) allocator.AddParticle(parent, output);      // add it only if its a primary
 259       allocator.AddParticle(p1, output);
 260       allocator.AddParticle(p2, output);
 261       success = kTRUE;
 262     }
 263
 264     // third case: three daughter particle
 265     else if(nSec == 3) {
 266       // initialize the daughter particle
 267       Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
 268       p1.Pos(parentBW.Pos());
 269       Particle p2(database->GetPDGParticle(dc->GetDaughterPDG(1)));
 270       p2.Pos(parentBW.Pos());
 271       Particle p3(database->GetPDGParticle(dc->GetDaughterPDG(2)));
 272       p3.Pos(parentBW.Pos());
 273       // calculate the momenta in the rest frame of the mother particle
 274       Double_t pAbs1 = 0., pAbs2 = 0., pAbs3 = 0., sumPabs = 0., maxPabs = 0.;
 275       Double_t mass1 = p1.TableMass(), mass2 = p2.TableMass(), mass3 = p3.TableMass();
 276       TLorentzVector &mom1 = p1.Mom(), &mom2 = p2.Mom(), &mom3 = p3.Mom();
 277       Double_t deltaMass = BWmass - mass1 - mass2 - mass3;
 278
 279       do {
 280         Double_t rd1 = gRandom->Rndm();
 281         Double_t rd2 = gRandom->Rndm();
 282         if (rd2 > rd1)
 283           std::swap(rd1, rd2);
 284         // 1
 285         Double_t e = rd2*deltaMass;
 286         pAbs1 = TMath::Sqrt(e*e + 2*e*mass1);
 287         sumPabs = pAbs1;
 288         maxPabs = sumPabs;
 289         // 2
 290         e = (1-rd1)*deltaMass;
 291         pAbs2 = TMath::Sqrt(e*e + 2*e*mass2);
 292
 293         if(pAbs2 > maxPabs)
 294           maxPabs = pAbs2;
 295
 296         sumPabs += pAbs2;
 297         // 3
 298         e = (rd1-rd2)*deltaMass;
 299         pAbs3 = TMath::Sqrt(e*e + 2*e*mass3);
 300
 301         if (pAbs3 > maxPabs)
 302           maxPabs =  pAbs3;
 303         sumPabs  +=  pAbs3;
 304       } while(maxPabs > sumPabs - maxPabs);
 305
 306       // isotropic sample first particle 3-momentum
 307       Double_t cosTheta = 2*(gRandom->Rndm()) - 1;
 308       Double_t sinTheta = TMath::Sqrt(1 - cosTheta*cosTheta);
 309       Double_t phi      = TMath::TwoPi()*(gRandom->Rndm());
 310       Double_t sinPhi   = TMath::Sin(phi);
 311       Double_t cosPhi   = TMath::Cos(phi);
 312
 313       mom1.SetPxPyPzE(sinTheta*cosPhi, sinTheta*sinPhi, cosTheta, 0);
 314       mom1 *= pAbs1;
 315       // sample rest particle 3-momentum
 316       Double_t cosThetaN = (pAbs2*pAbs2 - pAbs3*pAbs3 - pAbs1*pAbs1)/(2*pAbs1*pAbs3);
 317       Double_t sinThetaN = TMath::Sqrt(1 - cosThetaN*cosThetaN);
 318       Double_t phiN      = TMath::TwoPi()*(gRandom->Rndm());
 319       Double_t sinPhiN   = TMath::Sin(phiN);
 320       Double_t cosPhiN   = TMath::Cos(phiN);
 321
 322       mom3.SetPxPyPzE(sinThetaN*cosPhiN*cosTheta*cosPhi - sinThetaN*sinPhiN*sinPhi + cosThetaN*sinTheta*cosPhi,
 323                       sinThetaN*cosPhiN*cosTheta*sinPhi + sinThetaN*sinPhiN*cosPhi + cosThetaN*sinTheta*sinPhi,
 324                       -sinThetaN*cosPhiN*sinTheta + cosThetaN*cosTheta,
 325                       0.);
 326
 327       mom3 *= pAbs3*mom3.P();
 328       mom2 = mom1;
 329       mom2 += mom3;
 330       mom2 *= -1.;
 331       // calculate energy
 332       mom1.SetE(TMath::Sqrt(mom1.P()*mom1.P() + mass1*mass1));
 333       mom2.SetE(TMath::Sqrt(mom2.P()*mom2.P() + mass2*mass2));
 334       mom3.SetE(TMath::Sqrt(mom3.P()*mom3.P() + mass3*mass3));
 335
 336       // boost to Lab system
 337       mom1.Boost(velocityBW);
 338       mom2.Boost(velocityBW);
 339       mom3.Boost(velocityBW);
 340
 341       p1.SetLastMotherPdg(parentBW.Encoding());
 342       p1.SetLastMotherDecayCoor(parentBW.Pos());
 343       p1.SetLastMotherDecayMom(parentBW.Mom());
 344       p2.SetLastMotherPdg(parentBW.Encoding());
 345       p2.SetLastMotherDecayCoor(parentBW.Pos());
 346       p2.SetLastMotherDecayMom(parentBW.Mom());
 347       p3.SetLastMotherPdg(parentBW.Encoding());
 348       p3.SetLastMotherDecayCoor(parentBW.Pos());
 349       p3.SetLastMotherDecayMom(parentBW.Mom());
 350
 351       //set to daughters the same type as has mother
 352       p1.SetType(NB);
 353       p2.SetType(NB);
 354       p3.SetType(NB);
 355  //      std::cout<<"3d NB="<<NB<<std::endl;
 356
 357
 358
 359       // energy conservation check in the lab system
 360       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()) +
 361                                     (parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y()-p3.Mom().Y())*(parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y()-p3.Mom().Y()) +
 362                                     (parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z()-p3.Mom().Z())*(parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z()-p3.Mom().Z())     +
 363                                     (parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()-p3.Mom().E())*(parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()-p3.Mom().E()));
 364       // if deltaS is too big then repeat the kinematic procedure
 365       if(deltaS>0.001) {
 366
 367         cout << "3-body decay kinematic check in lab system: " << pDef->GetPDG() << " --> " << p1.Encoding() << " + " << p2.Encoding() << " + " << p3.Encoding() << endl;
 368         cout << "Mother    (e,px,py,pz): " << parentBW.Mom().E() << "\t" << parentBW.Mom().X() << "\t" << parentBW.Mom().Y() << "\t" << parentBW.Mom().Z() << endl;
 369         cout << "Daughter1 (e,px,py,pz): " << p1.Mom().E() << "\t" << p1.Mom().X() << "\t" << p1.Mom().Y() << "\t" << p1.Mom().Z() << endl;
 370         cout << "Daughter2 (e,px,py,pz): " << p2.Mom().E() << "\t" << p2.Mom().X() << "\t" << p2.Mom().Y() << "\t" << p2.Mom().Z() << endl;
 371         cout << "Daughter3 (e,px,py,pz): " << p3.Mom().E() << "\t" << p3.Mom().X() << "\t" << p3.Mom().Y() << "\t" << p3.Mom().Z() << endl;
 372         cout << "3-body decay delta(sqrtS) = " << deltaS << endl;
 373         cout << "Repeating the decay algorithm..." << endl;
 374
 375         iterations++;
 376         continue;
 377       }
 378
 379       Int_t parentIndex = -1;
 380       if(parent.GetMother()==-1) parentIndex = parent.SetIndex();   // parents which are primaries don't have yet an index
 381       else parentIndex = parent.GetIndex();                         // parents which are secondaries have an index
 382       p1.SetIndex();
 383       p2.SetIndex();
 384       p3.SetIndex();
 385       p1.SetMother(parentIndex);
 386       p2.SetMother(parentIndex);
 387       p3.SetMother(parentIndex);
 388       parent.SetDaughter(p1.GetIndex());
 389       parent.SetDaughter(p2.GetIndex());
 390       parent.SetDaughter(p3.GetIndex());
 391       if(parent.GetMother()==-1) allocator.AddParticle(parent, output);    // add it only if its a primary
 392       allocator.AddParticle(p1, output);
 393       allocator.AddParticle(p2, output);
 394       allocator.AddParticle(p3, output);
 395       success = kTRUE;
 396     }
 397   }
 398   return;
 399 }