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
  50   // if particle is set to be stable then return 0
  51   if(pDef->GetStableStatus())
  52     return 0.;
  53
  54   if(width > weakDecayLimit && weakDecayLimit>=0.0) {
  55     const Double_t slope =  parent.E() * 0.1973 / (pDef->GetMass() * width);
  56     return -slope * TMath::Log(gRandom->Rndm());//in fm/c
  57   }
  58
  59   return 0.;
  60 }
  61
  62
  63 extern "C" void mydelta_();
  64 extern SERVICEEVCommon SERVICEEV;
  65
  66 void Decay(List_t &output, Particle &parent, ParticleAllocator &allocator, DatabasePDG* database) {
  67   // check if the parent particle has been decayed already
  68   //  std::cout << "HadronDecayer::Decay() IN" << std::endl;
  69   Int_t daughters = parent.GetNDaughters();
  70
  71   //cout<<"in Decay pdg"<< parent.Def()->GetPDG()<<"daughters "<<daughters<<endl;
  72
  73   if(daughters>0)  // particle decayed already
  74     return;
  75
  76   // Get the PDG properties of the particle
  77   ParticlePDG *pDef = parent.Def();
  78
  79   // Get the number of posible decay channels
  80   Int_t nDecayChannel = pDef->GetNDecayChannels();
  81
  82   // check the full branching of this specie
  83   Double_t fullBranching = pDef->GetFullBranching(); // Only 3 or less body decays
  84
  85   // return if particle has no branching
  86   if(fullBranching < 0.00001)
  87     return;
  88
  89   // get the PDG mass of the specie
  90   Double_t PDGmass = pDef->GetMass();
  91   int ComprCodePyth=0;
  92   Float_t Delta =0;
  93
  94   Bool_t success = kFALSE;
  95   Int_t iterations = 0;
  96   // Try to decay the particle
  97   while(!success) {
  98     //    std::cout << "HadronDecayer::Decay() iteration #" << iterations << std::endl;
  99     if(iterations>1000) {
 100       std::cout << "HadronDecayer::Decay() more than 1000 iterations to decay particle with code "
 101                 << pDef->GetPDG() << std::endl;
 102       std::cout << "               Will be left undecayed ... check it out!" << std::endl;
 103       return;
 104     }
 105     // get a random mass using the Breit Wigner distribution
 106     Double_t BWmass = gRandom->BreitWigner(PDGmass, pDef->GetWidth());
 107     //!!!!
 108     //      BWmass = PDGmass;
 109     // Try to cut the Breit Wigner tail of the particle using the cuts from pythia
 110     // The Delta variable is obtained from pythia based on the specie
 111     int encoding =pDef->GetPDG();
 112     SERVICEEV.ipdg = encoding;
 113     mydelta_();
 114     ComprCodePyth=SERVICEEV.KC;
 115     Delta = SERVICEEV.delta;// PYDAT2.PMAS[KC][3];
 116
 117     //if there are no such particle in PYTHIA particle table, we take Delta=0.4
 118     if(ComprCodePyth==0){
 119       BWmass=PDGmass;
 120       Delta=0.0;
 121     }
 122
 123     //bad delta - an exception
 124     if(ComprCodePyth==254){
 125       BWmass=PDGmass;
 126       Delta=0.0;
 127     }
 128
 129     // K0 decay into K0s or K0l
 130     if(TMath::Abs(encoding)==311) {
 131       BWmass=PDGmass;
 132       Delta=0.0;
 133     }
 134
 135     //    std::cout << "HadronDecayer::Decay() BWmass = " << BWmass << std::endl;
 136     //    std::cout << "HadronDecayer::Decay() Delta = " << Delta << std::endl;
 137     //for particles from PYTHIA table only, if the BW mass is outside the cut range then quit this iteration and generate another BW mass
 138     if(ComprCodePyth!=0 && Delta>0 && (BWmass<PDGmass-Delta || BWmass>PDGmass+Delta)){
 139       //      std::cout<<"encoding"<<encoding<<"delta"<<Delta<<"width "<<pDef->GetWidth()<<"mass"<<BWmass<<std::endl;
 140       iterations++;
 141       //      std::cout << "HadronDecayer::Decay() BWmass outside cut, try again" << std::endl;
 142       continue;
 143     }
 144     //----
 145
 146     //    if(BWmass>5)
 147     //      std::cout<<" > 5 encoding"<<encoding<<" pdgmass "<<PDGmass<<" delta "<<Delta<<"width "<<pDef->GetWidth()<<" mass "<<BWmass<<"CC"<<ComprCodePyth<<std::endl;
 148
 149     // check how many decay channels are allowed with the generated mass
 150     Int_t nAllowedChannels = database->GetNAllowedChannels(pDef, BWmass);
 151     // if no decay channels are posible with this mass, then generate another BW mass
 152     if(nAllowedChannels==0) {
 153       iterations++;
 154       //      std::cout << "HadronDecayer::Decay() no decays allowed at this BW mass" << std::endl;
 155       continue;
 156     }
 157
 158     std::vector<Particle> apDaughter;
 159     std::vector<Double_t> dMass; //daughters'mass
 160     std::vector<Double_t> dMom;
 161     std::vector<Double_t> sm;
 162     std::vector<Double_t> rd;
 163
 164     // we need to choose an allowed decay channel
 165     Double_t randValue = gRandom->Rndm() * fullBranching;
 166     Int_t chosenChannel = 1000;
 167     Bool_t found = kFALSE;
 168     Int_t channelIterations = 0;
 169     //    std::cout << "HadronDecayer::Decay() randValue = " << randValue  << std::endl;
 170     while(!found) {
 171       //      std::cout << "HadronDecayer::Decay() channel iteration #" << channelIterations << std::endl;
 172       for(Int_t nChannel = 0; nChannel < nDecayChannel; ++nChannel) {
 173         randValue -= pDef->GetDecayChannel(nChannel)->GetBranching();
 174         //      std::cout << "HadronDecayer::Decay() channel #" << nChannel << " randValue = " << randValue << std::endl;
 175         if(randValue <= 0. && database->IsChannelAllowed(pDef->GetDecayChannel(nChannel), BWmass)) {
 176           //      std::cout << "HadronDecayer::Decay() channel found" << std::endl;
 177           chosenChannel = nChannel;
 178           found = kTRUE;
 179           break;
 180         }
 181       }
 182       channelIterations++;
 183     }
 184
 185     // get the PDG information for the chosen decay channel
 186     DecayChannel *dc = pDef->GetDecayChannel(chosenChannel);
 187     Int_t nSec = dc->GetNDaughters();
 188
 189     // Adjust the parent momentum four-vector for the MC generated Breit-Wigner mass
 190     Particle parentBW(database->GetPDGParticle(parent.Encoding()));
 191     parentBW.Pos(parent.Pos());
 192     Double_t BWenergy = TMath::Sqrt(parent.Mom().X()*parent.Mom().X() +
 193                                     parent.Mom().Y()*parent.Mom().Y() +
 194                                     parent.Mom().Z()*parent.Mom().Z() +
 195                                     BWmass*BWmass);
 196
 197     Int_t NB = (Int_t)parent.GetType(); //particle from jets
 198
 199     TLorentzVector MomparentBW(parent.Mom().X(), parent.Mom().Y(), parent.Mom().Z(), BWenergy);
 200     parentBW.Mom(MomparentBW);
 201     // take into account BW when calculating boost velocity (for wide resonances it matters)
 202     TVector3 velocityBW(parentBW.Mom().BoostVector());
 203
 204     // now we have an allowed decay
 205     // first case: one daughter particle
 206     if(nSec == 1) {
 207       // initialize the daughter particle
 208       Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
 209       p1.Pos(parentBW.Pos());
 210       p1.Mom(parent.Mom());
 211       p1.SetLastMotherPdg(parentBW.Encoding());
 212       p1.SetLastMotherDecayCoor(parentBW.Pos());
 213       p1.SetLastMotherDecayMom(parentBW.Mom());
 214       p1.SetType(NB);
 215
 216       // link the parent and daughters trough their indexes in the list
 217       Int_t parentIndex = -1;
 218       if(parent.GetMother()==-1) parentIndex = parent.SetIndex();   // parents which are primaries don't have yet an index
 219       else parentIndex = parent.GetIndex();                         // parents which are secondaries have an index
 220       Int_t p1Index = p1.SetIndex();                           // set the daughter index
 221       p1.SetMother(parentIndex);                               // set the mother index for this daughter
 222       parent.SetDaughter(p1Index);                             // set p1 as daughter to the parent
 223       if(parent.GetMother()==-1) allocator.AddParticle(parent, output);   // add it only if its a primary particle
 224       allocator.AddParticle(p1, output);
 225       success = kTRUE;
 226     }
 227     // second case: two daughter particles
 228     else if(nSec == 2) {
 229       // initialize the daughter particles
 230       Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
 231       p1.Pos(parentBW.Pos());
 232       Particle p2(database->GetPDGParticle(dc->GetDaughterPDG(1)));
 233       p2.Pos(parentBW.Pos());
 234
 235       // calculate the momenta in rest frame of mother for the two particles (theta and phi are isotropic)
 236       MomAntiMom(p1.Mom(), p1.TableMass(), p2.Mom(), p2.TableMass(), BWmass);
 237
 238       // boost to the laboratory system (to the mother velocity)
 239       p1.Mom().Boost(velocityBW);
 240       p2.Mom().Boost(velocityBW);
 241
 242       //store information about mother
 243       p1.SetLastMotherPdg(parentBW.Encoding());
 244       p1.SetLastMotherDecayCoor(parentBW.Pos());
 245       p1.SetLastMotherDecayMom(parentBW.Mom());
 246       p2.SetLastMotherPdg(parentBW.Encoding());
 247       p2.SetLastMotherDecayCoor(parentBW.Pos());
 248       p2.SetLastMotherDecayMom(parentBW.Mom());
 249       // std::cout<<"2d NB="<<NB<<std::endl;
 250       //set to daughters the same type as has mother
 251       p1.SetType(NB);
 252       p2.SetType(NB);
 253
 254
 255       // check the kinematics in the lab system
 256       Double_t deltaS = TMath::Sqrt((parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X())*(parentBW.Mom().X()-p1.Mom().X()-p2.Mom().X())+
 257                                     (parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y())*(parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y())+
 258                                     (parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z())*(parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z())+
 259                                     (parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E())*(parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()));
 260       // if deltaS is too big then repeat the kinematic procedure
 261
 262
 263       if(deltaS>0.001) {
 264
 265         //      cout << "2-body decay kinematic check in lab system: " << pDef->GetPDG() << " --> " << p1.Encoding() << " + " << p2.Encoding() << endl;
 266         //      cout << "Mother    (e,px,py,pz): " << parentBW.Mom().E() << "\t" << parentBW.Mom().X() << "\t" << parentBW.Mom().Y() << "\t" << parentBW.Mom().Z() << endl;
 267         //      cout << "Mother    (x,y,z,t): " << parentBW.Pos().X() << "\t" << parentBW.Pos().Y() << "\t" << parentBW.Pos().Z() << "\t" << parentBW.Pos().T() << endl;
 268
 269         //      cout << "Daughter1 (e,px,py,pz): " << p1.Mom().E() << "\t" << p1.Mom().X() << "\t" << p1.Mom().Y() << "\t" << p1.Mom().Z() << endl;
 270         //      cout << "Daughter2 (e,px,py,pz): " << p2.Mom().E() << "\t" << p2.Mom().X() << "\t" << p2.Mom().Y() << "\t" << p2.Mom().Z() << endl;
 271         //      cout << "2-body decay delta(sqrtS) = " << deltaS << endl;
 272         //      cout << "Repeating the decay algorithm ..." << endl;
 273
 274         iterations++;
 275         continue;
 276       }
 277       // push particles to the list of secondaries
 278       Int_t parentIndex = -1;
 279       if(parent.GetMother()==-1) parentIndex = parent.SetIndex();   // parents which are primaries don't have yet an index
 280       else parentIndex = parent.GetIndex();                         // parents which are secondaries have an index
 281       p1.SetIndex();
 282       p2.SetIndex();
 283       p1.SetMother(parentIndex);
 284       p2.SetMother(parentIndex);
 285       parent.SetDaughter(p1.GetIndex());
 286       parent.SetDaughter(p2.GetIndex());
 287       if(parent.GetMother()==-1) allocator.AddParticle(parent, output);      // add it only if its a primary
 288       allocator.AddParticle(p1, output);
 289       allocator.AddParticle(p2, output);
 290       success = kTRUE;
 291     }
 292
 293     // third case: three daughter particle
 294     else if(nSec == 3) {
 295       // initialize the daughter particle
 296       Particle p1(database->GetPDGParticle(dc->GetDaughterPDG(0)));
 297       p1.Pos(parentBW.Pos());
 298       Particle p2(database->GetPDGParticle(dc->GetDaughterPDG(1)));
 299       p2.Pos(parentBW.Pos());
 300       Particle p3(database->GetPDGParticle(dc->GetDaughterPDG(2)));
 301       p3.Pos(parentBW.Pos());
 302       // calculate the momenta in the rest frame of the mother particle
 303       Double_t pAbs1 = 0., pAbs2 = 0., pAbs3 = 0., sumPabs = 0., maxPabs = 0.;
 304       Double_t mass1 = p1.TableMass(), mass2 = p2.TableMass(), mass3 = p3.TableMass();
 305       TLorentzVector &mom1 = p1.Mom(), &mom2 = p2.Mom(), &mom3 = p3.Mom();
 306       Double_t deltaMass = BWmass - mass1 - mass2 - mass3;
 307
 308       do {
 309         Double_t rd1 = gRandom->Rndm();
 310         Double_t rd2 = gRandom->Rndm();
 311         if (rd2 > rd1)
 312           std::swap(rd1, rd2);
 313         // 1
 314         Double_t e = rd2*deltaMass;
 315         pAbs1 = TMath::Sqrt(e*e + 2*e*mass1);
 316         sumPabs = pAbs1;
 317         maxPabs = sumPabs;
 318         // 2
 319         e = (1-rd1)*deltaMass;
 320         pAbs2 = TMath::Sqrt(e*e + 2*e*mass2);
 321
 322         if(pAbs2 > maxPabs)
 323           maxPabs = pAbs2;
 324
 325         sumPabs += pAbs2;
 326         // 3
 327         e = (rd1-rd2)*deltaMass;
 328         pAbs3 = TMath::Sqrt(e*e + 2*e*mass3);
 329
 330         if (pAbs3 > maxPabs)
 331           maxPabs =  pAbs3;
 332         sumPabs  +=  pAbs3;
 333       } while(maxPabs > sumPabs - maxPabs);
 334
 335       // isotropic sample first particle 3-momentum
 336       Double_t cosTheta = 2*(gRandom->Rndm()) - 1;
 337       Double_t sinTheta = TMath::Sqrt(1 - cosTheta*cosTheta);
 338       Double_t phi      = TMath::TwoPi()*(gRandom->Rndm());
 339       Double_t sinPhi   = TMath::Sin(phi);
 340       Double_t cosPhi   = TMath::Cos(phi);
 341
 342       mom1.SetPxPyPzE(sinTheta*cosPhi, sinTheta*sinPhi, cosTheta, 0);
 343       mom1 *= pAbs1;
 344       // sample rest particle 3-momentum
 345       Double_t cosThetaN = (pAbs2*pAbs2 - pAbs3*pAbs3 - pAbs1*pAbs1)/(2*pAbs1*pAbs3);
 346       Double_t sinThetaN = TMath::Sqrt(1 - cosThetaN*cosThetaN);
 347       Double_t phiN      = TMath::TwoPi()*(gRandom->Rndm());
 348       Double_t sinPhiN   = TMath::Sin(phiN);
 349       Double_t cosPhiN   = TMath::Cos(phiN);
 350
 351       mom3.SetPxPyPzE(sinThetaN*cosPhiN*cosTheta*cosPhi - sinThetaN*sinPhiN*sinPhi + cosThetaN*sinTheta*cosPhi,
 352                       sinThetaN*cosPhiN*cosTheta*sinPhi + sinThetaN*sinPhiN*cosPhi + cosThetaN*sinTheta*sinPhi,
 353                       -sinThetaN*cosPhiN*sinTheta + cosThetaN*cosTheta,
 354                       0.);
 355
 356       mom3 *= pAbs3*mom3.P();
 357       mom2 = mom1;
 358       mom2 += mom3;
 359       mom2 *= -1.;
 360       // calculate energy
 361       mom1.SetE(TMath::Sqrt(mom1.P()*mom1.P() + mass1*mass1));
 362       mom2.SetE(TMath::Sqrt(mom2.P()*mom2.P() + mass2*mass2));
 363       mom3.SetE(TMath::Sqrt(mom3.P()*mom3.P() + mass3*mass3));
 364
 365       // boost to Lab system
 366       mom1.Boost(velocityBW);
 367       mom2.Boost(velocityBW);
 368       mom3.Boost(velocityBW);
 369
 370       p1.SetLastMotherPdg(parentBW.Encoding());
 371       p1.SetLastMotherDecayCoor(parentBW.Pos());
 372       p1.SetLastMotherDecayMom(parentBW.Mom());
 373       p2.SetLastMotherPdg(parentBW.Encoding());
 374       p2.SetLastMotherDecayCoor(parentBW.Pos());
 375       p2.SetLastMotherDecayMom(parentBW.Mom());
 376       p3.SetLastMotherPdg(parentBW.Encoding());
 377       p3.SetLastMotherDecayCoor(parentBW.Pos());
 378       p3.SetLastMotherDecayMom(parentBW.Mom());
 379
 380       //set to daughters the same type as has mother
 381       p1.SetType(NB);
 382       p2.SetType(NB);
 383       p3.SetType(NB);
 384  //      std::cout<<"3d NB="<<NB<<std::endl;
 385
 386
 387
 388       // energy conservation check in the lab system
 389       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()) +
 390                                     (parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y()-p3.Mom().Y())*(parentBW.Mom().Y()-p1.Mom().Y()-p2.Mom().Y()-p3.Mom().Y()) +
 391                                     (parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z()-p3.Mom().Z())*(parentBW.Mom().Z()-p1.Mom().Z()-p2.Mom().Z()-p3.Mom().Z())     +
 392                                     (parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()-p3.Mom().E())*(parentBW.Mom().E()-p1.Mom().E()-p2.Mom().E()-p3.Mom().E()));
 393       // if deltaS is too big then repeat the kinematic procedure
 394       if(deltaS>0.001) {
 395
 396         //      cout << "3-body decay kinematic check in lab system: " << pDef->GetPDG() << " --> " << p1.Encoding() << " + " << p2.Encoding() << " + " << p3.Encoding() << endl;
 397         //      cout << "Mother    (e,px,py,pz): " << parentBW.Mom().E() << "\t" << parentBW.Mom().X() << "\t" << parentBW.Mom().Y() << "\t" << parentBW.Mom().Z() << endl;
 398         //      cout << "Daughter1 (e,px,py,pz): " << p1.Mom().E() << "\t" << p1.Mom().X() << "\t" << p1.Mom().Y() << "\t" << p1.Mom().Z() << endl;
 399         //      cout << "Daughter2 (e,px,py,pz): " << p2.Mom().E() << "\t" << p2.Mom().X() << "\t" << p2.Mom().Y() << "\t" << p2.Mom().Z() << endl;
 400         //      cout << "Daughter3 (e,px,py,pz): " << p3.Mom().E() << "\t" << p3.Mom().X() << "\t" << p3.Mom().Y() << "\t" << p3.Mom().Z() << endl;
 401         //      cout << "3-body decay delta(sqrtS) = " << deltaS << endl;
 402         //      cout << "Repeating the decay algorithm..." << endl;
 403
 404         iterations++;
 405         continue;
 406       }
 407
 408       Int_t parentIndex = -1;
 409       if(parent.GetMother()==-1) parentIndex = parent.SetIndex();   // parents which are primaries don't have yet an index
 410       else parentIndex = parent.GetIndex();                         // parents which are secondaries have an index
 411       p1.SetIndex();
 412       p2.SetIndex();
 413       p3.SetIndex();
 414       p1.SetMother(parentIndex);
 415       p2.SetMother(parentIndex);
 416       p3.SetMother(parentIndex);
 417       parent.SetDaughter(p1.GetIndex());
 418       parent.SetDaughter(p2.GetIndex());
 419       parent.SetDaughter(p3.GetIndex());
 420       if(parent.GetMother()==-1) allocator.AddParticle(parent, output);    // add it only if its a primary
 421       allocator.AddParticle(p1, output);
 422       allocator.AddParticle(p2, output);
 423       allocator.AddParticle(p3, output);
 424       success = kTRUE;
 425     }
 426   }
 427   return;
 428 }