EVGEN/AliGenLightNuclei.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2009-2014, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 ////////////////////////////////////////////////////////////////////
  17 //
  18 // Afterburner to generate light nuclei for event generators
  19 // such as PYTHIA and PHOJET
  20 //
  21 // Light nuclei are generated whenever a cluster of nucleons is found
  22 // within a sphere of radius p0 (coalescence momentum), i.e. have the
  23 // same momentum.
  24 //
  25 // By default it starts with He4 nuclei which are the most stable,
  26 // then He3 nuclei, tritons and finally deuterons. It can also generate
  27 // a single nucleus species by disabling the others.
  28 //
  29 // Sample code for PYTHIA:
  30 //
  31 //    AliGenLightNuclei* gener = new AliGenLightNuclei();
  32 //
  33 //    AliGenPythia* pythia = new AliGenPythia(-1);
  34 //    pythia->SetCollisionSystem("p+", "p+");
  35 //    pythia->SetEnergyCMS(7000);
  36 //
  37 //    gener->UsePerEventRates();
  38 //    gener->AddGenerator(pythia, "PYTHIA", 1);
  39 //    gener->SetCoalescenceMomentum(0.100); // default (GeV/c)
  40 //
  41 //////////////////////////////////////////////////////////////////////
  42
  43 //
  44 // Author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
  45 //
  46
  47 #include "TMath.h"
  48 #include "TPDGCode.h"
  49 #include "TMCProcess.h"
  50 #include "TList.h"
  51 #include "TVector3.h"
  52 #include "TParticle.h"
  53 #include "AliStack.h"
  54 #include "AliRun.h"
  55 #include "AliLog.h"
  56 #include "AliMC.h"
  57 #include "AliGenLightNuclei.h"
  58
  59 ClassImp(AliGenLightNuclei)
  60
  61 AliGenLightNuclei::AliGenLightNuclei()
  62 :AliGenCocktail()
  63 ,fP0(0.100)
  64 ,fGenDeuterons(kTRUE)
  65 ,fGenTritons(kTRUE)
  66 ,fGenHe3Nuclei(kTRUE)
  67 ,fGenHe4Nuclei(kTRUE)
  68 {
  69 //
  70 // default constructor
  71 //
  72 }
  73
  74 AliGenLightNuclei::~AliGenLightNuclei()
  75 {
  76 //
  77 // default destructor
  78 //
  79 }
  80
  81 void AliGenLightNuclei::Generate()
  82 {
  83 //
  84 // delegate the particle generation to the cocktail
  85 // and modify the stack adding the light nuclei
  86 //
  87         AliGenCocktail::Generate();
  88
  89         // find the nucleons and anti-nucleons
  90
  91         TList* protons      = new TList();
  92         TList* neutrons     = new TList();
  93         TList* antiprotons  = new TList();
  94         TList* antineutrons = new TList();
  95
  96         for (Int_t i=0; i < fStack->GetNprimary(); ++i)
  97         {
  98                 TParticle* iParticle = fStack->Particle(i);
  99
 100                 if(iParticle->GetStatusCode() != 1) continue;
 101
 102                 switch(iParticle->GetPdgCode())
 103                 {
 104                         case kProton:
 105                                 protons->Add(iParticle);
 106                                 break;
 107                         case kNeutron:
 108                                 neutrons->Add(iParticle);
 109                                 break;
 110                         case kProtonBar:
 111                                 antiprotons->Add(iParticle);
 112                                 break;
 113                         case kNeutronBar:
 114                                 antineutrons->Add(iParticle);
 115                                 break;
 116                         default:
 117                                 break;
 118                 }
 119         }
 120
 121         // do not delete content
 122         protons->SetOwner(kFALSE);
 123         neutrons->SetOwner(kFALSE);
 124         antiprotons->SetOwner(kFALSE);
 125         antineutrons->SetOwner(kFALSE);
 126
 127         // first try with He4 nuclei which are the most stable
 128
 129         if(fGenHe4Nuclei)
 130         {
 131                 this->GenerateHe4Nuclei(protons, neutrons);
 132                 this->GenerateHe4Nuclei(antiprotons, antineutrons);
 133         }
 134
 135         // then He3 nuclei
 136
 137         if(fGenHe3Nuclei)
 138         {
 139                 this->GenerateHe3Nuclei(protons, neutrons);
 140                 this->GenerateHe3Nuclei(antiprotons, antineutrons);
 141         }
 142
 143         // then tritons
 144
 145         if(fGenTritons)
 146         {
 147                 this->GenerateTritons(protons, neutrons);
 148                 this->GenerateTritons(antiprotons, antineutrons);
 149         }
 150
 151         // and finally deuterons
 152
 153         if(fGenDeuterons)
 154         {
 155                 this->GenerateDeuterons(protons, neutrons);
 156                 this->GenerateDeuterons(antiprotons, antineutrons);
 157         }
 158
 159         protons->Clear("nodelete");
 160         neutrons->Clear("nodelete");
 161         antiprotons->Clear("nodelete");
 162         antineutrons->Clear("nodelete");
 163
 164         delete protons;
 165         delete neutrons;
 166         delete antiprotons;
 167         delete antineutrons;
 168 }
 169
 170 Bool_t AliGenLightNuclei::Coalescence(const TParticle* n1, const TParticle* n2) const
 171 {
 172 //
 173 // returns true if the nucleons are inside of an sphere of radius p0
 174 // (assume the nucleons are in the same place e.g. PYTHIA, PHOJET,...)
 175 //
 176         Double_t deltaP = this->GetPcm(n1->Px(), n1->Py(), n1->Pz(), n1->GetMass(),
 177                                        n2->Px(), n2->Py(), n2->Pz(), n2->GetMass());
 178
 179         return ( deltaP < fP0);
 180 }
 181
 182 TParticle* AliGenLightNuclei::FindPartner(const TParticle* n0, const TList* nucleons, const TParticle* nx) const
 183 {
 184 //
 185 // find the first nucleon partner within a sphere of radius p0
 186 // centered at n0 and exclude nucleon nx
 187 //
 188         TIter n_next(nucleons);
 189         while(TParticle* n1 = dynamic_cast<TParticle*>( n_next()) )
 190         {
 191                 if(n1 == 0) continue;
 192                 if(n1 == nx) continue;
 193                 if(n1->GetStatusCode() == kCluster) continue;
 194                 if( !this->Coalescence(n0, n1) ) continue;
 195                 return n1;
 196         }
 197
 198         return 0;
 199 }
 200
 201 Int_t AliGenLightNuclei::GenerateDeuterons(const TList* protons, const TList* neutrons)
 202 {
 203 //
 204 // a deuteron is generated from a pair of p-n nucleons
 205 // (the center of the sphere is one of the nucleons)
 206 //
 207         Int_t npart = 0;
 208
 209         TIter p_next(protons);
 210         while(TParticle* n0 = dynamic_cast<TParticle*>(p_next()) )
 211         {
 212                 if(n0 == 0) continue;
 213                 if(n0->GetStatusCode() == kCluster) continue;
 214
 215                 TParticle* partner = this->FindPartner(n0, neutrons, 0);
 216
 217                 if(partner == 0) continue;
 218
 219                 this->PushDeuteron(n0, partner);
 220
 221                 ++npart;
 222         }
 223
 224         return npart;
 225 }
 226
 227 Int_t AliGenLightNuclei::GenerateTritons(const TList* protons, const TList* neutrons)
 228 {
 229 //
 230 // a triton is generated from a cluster of p-n-n nucleons with same momentum
 231 // (triangular configuration)
 232 //
 233         Int_t npart = 0;
 234
 235         TIter p_next(protons);
 236         while(TParticle* n0 = dynamic_cast<TParticle*>(p_next()) )
 237         {
 238                 if(n0 == 0) continue;
 239                 if(n0->GetStatusCode() == kCluster) continue;
 240
 241                 TParticle* partner1 = this->FindPartner(n0, neutrons, 0);
 242
 243                 if(partner1 == 0) continue;
 244
 245                 TParticle* partner2 = this->FindPartner(n0, neutrons, partner1);
 246
 247                 if(partner2 == 0) continue;
 248
 249                 // check that the partners coalesce between themselves
 250
 251                 if(!this->Coalescence(partner1, partner2)) continue;
 252
 253                 this->PushTriton(n0, partner1, partner2);
 254
 255                 ++npart;
 256         }
 257
 258         return npart;
 259 }
 260
 261 Int_t AliGenLightNuclei::GenerateHe3Nuclei(const TList* protons, const TList* neutrons)
 262 {
 263 //
 264 // a He3 nucleus is generated from a cluster of p-n-p nucleons with same momentum
 265 // (triangular configuration)
 266 //
 267         Int_t npart = 0;
 268
 269         TIter p_next(protons); // center of the sphere
 270         while(TParticle* n0 = dynamic_cast<TParticle*>(p_next()) )
 271         {
 272                 if(n0 == 0) continue;
 273                 if(n0->GetStatusCode() == kCluster) continue;
 274
 275                 TParticle* partner1 = this->FindPartner(n0, neutrons, 0);
 276
 277                 if(partner1 == 0) continue;
 278
 279                 TParticle* partner2 = this->FindPartner(n0, protons, n0);
 280
 281                 if(partner2 == 0) continue;
 282
 283                 // check that the partners coalesce between themselves
 284
 285                 if(!this->Coalescence(partner1, partner2)) continue;
 286
 287                 this->PushHe3Nucleus(n0, partner1, partner2);
 288
 289                 ++npart;
 290         }
 291
 292         return npart;
 293 }
 294
 295 Int_t AliGenLightNuclei::GenerateHe4Nuclei(const TList* protons, const TList* neutrons)
 296 {
 297 //
 298 // a He4 nucleus is generated from a cluster of p-n-p-n nucleons with same momentum
 299 // (tetrahedron configuration)
 300 //
 301         Int_t npart = 0;
 302
 303         TIter p_next(protons); // center of the sphere
 304         while(TParticle* n0 = dynamic_cast<TParticle*>(p_next()) )
 305         {
 306                 if(n0 == 0) continue;
 307                 if(n0->GetStatusCode() == kCluster) continue;
 308
 309                 TParticle* partner1 = this->FindPartner(n0, neutrons, 0);
 310
 311                 if(partner1 == 0) continue;
 312
 313                 TParticle* partner2 = this->FindPartner(n0, protons, n0);
 314
 315                 if(partner2 == 0) continue;
 316
 317                 TParticle* partner3 = this->FindPartner(n0, neutrons, partner1);
 318
 319                 if(partner3 == 0) continue;
 320
 321                 // check that the partners coalesce between themselves
 322
 323                 if(!this->Coalescence(partner1, partner2)) continue;
 324                 if(!this->Coalescence(partner1, partner3)) continue;
 325                 if(!this->Coalescence(partner2, partner3)) continue;
 326
 327                 this->PushHe4Nucleus(n0, partner1, partner2, partner3 );
 328
 329                 ++npart;
 330         }
 331
 332         return npart;
 333 }
 334
 335 void AliGenLightNuclei::PushDeuteron(TParticle* parent1, TParticle* parent2)
 336 {
 337 //
 338 // push a deuteron to the particle stack
 339 //
 340         Int_t pdg = ( parent1->GetPdgCode() > 0 ) ? kDeuteron : kAntiDeuteron;
 341         this->PushNucleus(pdg, 1.87561282, parent1, parent2);
 342 }
 343
 344 void AliGenLightNuclei::PushTriton(TParticle* parent1, TParticle* parent2, TParticle* parent3)
 345 {
 346 //
 347 // push a triton to the particle stack
 348 //
 349         Int_t pdg = ( parent1->GetPdgCode() > 0 ) ? kTriton : kAntiTriton;
 350         this->PushNucleus(pdg, 2.80925, parent1, parent2, parent3);
 351 }
 352
 353 void AliGenLightNuclei::PushHe3Nucleus(TParticle* parent1, TParticle* parent2, TParticle* parent3)
 354 {
 355 //
 356 // push a He3 nucleus to the particle stack
 357 //
 358         Int_t pdg = ( parent1->GetPdgCode() > 0 ) ? kHe3Nucleus : kAntiHe3Nucleus;
 359         this->PushNucleus(pdg, 2.80923, parent1, parent2, parent3);
 360 }
 361
 362 void AliGenLightNuclei::PushHe4Nucleus(TParticle* parent1, TParticle* parent2, TParticle* parent3, TParticle* parent4)
 363 {
 364 //
 365 // push a He4 nucleus to the particle stack
 366 //
 367         Int_t pdg = ( parent1->GetPdgCode() > 0 ) ? kAlpha : kAntiAlpha;
 368         this->PushNucleus(pdg, 3.727417, parent1, parent2, parent3, parent4);
 369 }
 370
 371 void AliGenLightNuclei::PushNucleus(Int_t pdg, Double_t mass, TParticle* parent1, TParticle* parent2, TParticle* parent3, TParticle* parent4)
 372 {
 373 //
 374 // push a nucleus to the stack and tag the parents with the kCluster status code
 375 //
 376         Int_t ntrk;
 377
 378         // momentum
 379         TVector3 p1(parent1->Px(), parent1->Py(), parent1->Pz());
 380         TVector3 p2(parent2->Px(), parent2->Py(), parent2->Pz());
 381         TVector3 p3(0, 0, 0);
 382         TVector3 p4(0, 0, 0);
 383         if(parent3 != 0) p3.SetXYZ(parent3->Px(), parent3->Py(), parent3->Pz());
 384         if(parent4 != 0) p4.SetXYZ(parent4->Px(), parent4->Py(), parent4->Pz());
 385
 386         // momentum
 387         TVector3 pN = p1+p2+p3+p4;
 388
 389         // E^2 = p^2 + m^2
 390         Double_t energy = TMath::Sqrt(pN.Mag2() + mass*mass);
 391
 392         // production vertex same as the parent1'
 393         TVector3 vN(parent1->Vx(), parent1->Vy(), parent1->Vz());
 394
 395         Double_t weight = 1;
 396         Int_t is = 1; // final state particle
 397
 398         // add a new nucleus to current event stack
 399         fStack->PushTrack(1, -1, pdg,
 400                          pN.X(), pN.Y(), pN.Z(), energy,
 401                          vN.X(), vN.Y(), vN.Z(), parent1->T(),
 402                          0., 0., 0., kPNCapture, ntrk, weight, is);
 403
 404         // change the status code of the parents
 405         parent1->SetStatusCode(kCluster);
 406         parent2->SetStatusCode(kCluster);
 407         if(parent3 != 0) parent3->SetStatusCode(kCluster);
 408         if(parent4 != 0) parent4->SetStatusCode(kCluster);
 409
 410         // set no transport for the parents
 411         parent1->SetBit(kDoneBit);
 412         parent2->SetBit(kDoneBit);
 413         if(parent3 != 0) parent3->SetBit(kDoneBit);
 414         if(parent4 != 0) parent4->SetBit(kDoneBit);
 415 }
 416
 417 Double_t AliGenLightNuclei::GetS(Double_t p1x, Double_t p1y, Double_t p1z, Double_t m1, Double_t p2x, Double_t p2y, Double_t p2z, Double_t m2) const
 418 {
 419 //
 420 // square of the center of mass energy
 421 //
 422         Double_t E1 = TMath::Sqrt( p1x*p1x + p1y*p1y + p1z*p1z + m1*m1);
 423         Double_t E2 = TMath::Sqrt( p2x*p2x + p2y*p2y + p2z*p2z + m2*m2);
 424
 425         return (E1+E2)*(E1+E2) - ((p1x+p2x)*(p1x+p2x) + (p1y+p2y)*(p1y+p2y) + (p1z+p2z)*(p1z+p2z));
 426 }
 427
 428 Double_t AliGenLightNuclei::GetPcm(Double_t p1x, Double_t p1y, Double_t p1z, Double_t m1, Double_t p2x, Double_t p2y, Double_t p2z, Double_t m2) const
 429 {
 430 //
 431 // momentum in the CM frame for 2 particles
 432 //
 433         Double_t s = this->GetS(p1x, p1y, p1z, m1, p2x, p2y, p2z, m2);
 434
 435         return TMath::Sqrt( (s-(m1-m2)*(m1-m2))*(s-(m1+m2)*(m1+m2)) )/(2.*TMath::Sqrt(s));
 436 }