EVGEN/AliGenDeuteron.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2009-2010, 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 // Afterburner to generate (anti)deuterons simulating coalescence of
  18 // (anti)nucleons as in A. J. Baltz et al., Phys. lett B 325(1994)7.
  19 // If the nucleon generator does not provide a spatial description of
  20 // the source the afterburner can provide one.
  21 //
  22 // There two options for the source: a thermal picture where all nucleons
  23 // are placed randomly and homogeneously in an spherical volume, or
  24 // an expansion one where they are projected onto its surface.
  25 //
  26 // An (anti)deuteron will form if there is a pair of (anti)proton-(anti)neutron
  27 // with momentum difference less than ~ 300MeV and relative distance less than
  28 // ~ 2.1fm. Only 3/4 of these clusters are accepted due to the deuteron spin.
  29 //
  30 // When there are more than one pair fulfilling the coalescence conditions,
  31 // the selected pair can be the one with the first partner, with
  32 // the lowest relative momentum, the lowest relative distance or both.
  33 //////////////////////////////////////////////////////////////////////////
  34
  35 // Author: Eulogio Serradilla <eulogio.serradilla@ciemat.es>,
  36 //         Arturo Menchaca <menchaca@fisica.unam.mx>
  37 //
  38
  39 #include "Riostream.h"
  40 #include "TParticle.h"
  41 #include "AliStack.h"
  42 #include "AliRun.h"
  43 #include "TMath.h"
  44 #include "TMCProcess.h"
  45 #include "TList.h"
  46 #include "TVector3.h"
  47 #include "AliMC.h"
  48 #include "TArrayF.h"
  49 #include "AliCollisionGeometry.h"
  50 #include "AliGenCocktailEventHeader.h"
  51 #include "AliGenCocktailEntry.h"
  52 #include "AliGenCocktailAfterBurner.h"
  53 #include "AliGenDeuteron.h"
  54 #include "AliLog.h"
  55
  56 ClassImp(AliGenDeuteron)
  57
  58 AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Int_t cluster)
  59  :fSign(1)
  60  ,fPmax(pmax)
  61  ,fRmax(rmax)
  62  ,fSpinProb(1)
  63  ,fClusterType(cluster)
  64  ,fModel(0)
  65  ,fTimeLength(2.5)
  66  ,fB(0)
  67  ,fR(0)
  68  ,fPsiR(0)
  69  ,fCurStack(0)
  70 {
  71 //
  72 // constructor
  73 //
  74         fSign = sign > 0 ? 1:-1;
  75
  76 }
  77
  78 AliGenDeuteron::~AliGenDeuteron()
  79 {
  80 //
  81 // Default destructor
  82 //
  83 }
  84
  85 void AliGenDeuteron::Init()
  86 {
  87 //
  88 // Standard AliGenerator initializer
  89 //
  90 }
  91
  92 void AliGenDeuteron::Generate()
  93 {
  94 //
  95 // Look for coalescence of (anti)nucleons in the nucleon source
  96 // provided by the generator cocktail
  97 //
  98         AliInfo(Form("sign: %d, Pmax: %g GeV/c, Rmax: %g fm, cluster: %d",fSign, fPmax, fRmax, fClusterType));
  99         if(fModel!=kNone) AliInfo(Form("model: %d, time: %g fm/c", fModel, fTimeLength));
 100
 101         AliGenCocktailAfterBurner* gener = (AliGenCocktailAfterBurner*)gAlice->GetMCApp()->Generator();
 102
 103         // find nuclear radius, only for first generator and projectile=target
 104         Bool_t collisionGeometry=0;
 105         if(fModel != kNone && gener->FirstGenerator()->Generator()->ProvidesCollisionGeometry())
 106         {
 107                 TString name;
 108                 Int_t ap, zp, at, zt;
 109                 gener->FirstGenerator()->Generator()->GetProjectile(name,ap,zp);
 110                 gener->FirstGenerator()->Generator()->GetTarget(name,at,zt);
 111                 if(ap != at) AliWarning("projectile and target have different size");
 112                 fR = 1.29*TMath::Power(at, 1./3.);
 113                 collisionGeometry = 1;
 114         }
 115
 116         for(Int_t ns = 0; ns < gener->GetNumberOfEvents(); ++ns)
 117         {
 118                 gener->SetActiveEventNumber(ns);
 119
 120                 if(fModel != kNone && collisionGeometry)
 121                 {
 122                         fPsiR = (gener->GetCollisionGeometry(ns))->ReactionPlaneAngle();
 123                         fB = (gener->GetCollisionGeometry(ns))->ImpactParameter();
 124
 125                         if(fB >= 2.*fR) continue; // no collision
 126                 }
 127
 128                 fCurStack = gener->GetStack(ns);
 129                 if(!fCurStack)
 130                 {
 131                         AliWarning("no event stack");
 132                         return;
 133                 }
 134
 135                 TList* protons = new TList();
 136                 protons->SetOwner(kFALSE);
 137                 TList* neutrons = new TList();
 138                 neutrons->SetOwner(kFALSE);
 139
 140                 // particles produced by the generator
 141                 for (Int_t i=0; i < fCurStack->GetNprimary(); ++i)
 142                 {
 143                         TParticle* iParticle = fCurStack->Particle(i);
 144
 145                         if(iParticle->GetStatusCode() != 1) continue;
 146
 147                         Int_t pdgCode = iParticle->GetPdgCode();
 148                         if(pdgCode == fSign*2212)// (anti)proton
 149                         {
 150                                 FixProductionVertex(iParticle);
 151                                 protons->Add(iParticle);
 152                         }
 153                         else if(pdgCode == fSign*2112) // (anti)neutron
 154                         {
 155                                 FixProductionVertex(iParticle);
 156                                 neutrons->Add(iParticle);
 157                         }
 158                 }
 159
 160                 // look for clusters
 161                 if(fClusterType==kFirstPartner)
 162                 {
 163                         FirstPartner(protons, neutrons);
 164                 }
 165                 else
 166                 {
 167                         WeightMatrix(protons, neutrons);
 168                 }
 169
 170                 protons->Clear("nodelete");
 171                 neutrons->Clear("nodelete");
 172
 173                 delete protons;
 174                 delete neutrons;
 175         }
 176
 177         AliInfo("DONE");
 178 }
 179
 180 Double_t AliGenDeuteron::GetCoalescenceProbability(const TParticle* nucleon1, const TParticle* nucleon2) const
 181 {
 182 //
 183 // Coalescence conditions as in
 184 // A. J. Baltz et al., Phys. lett B 325(1994)7
 185 //
 186 // returns < 0 if coalescence is not possible
 187 // otherwise returns a coalescence probability
 188 //
 189         const Double_t kProtonMass  = 0.938272013;
 190         const Double_t kNeutronMass = 0.939565378;
 191
 192         TVector3 v1(nucleon1->Vx(), nucleon1->Vy(), nucleon1->Vz());
 193         TVector3 p1(nucleon1->Px(), nucleon1->Py(), nucleon1->Pz());
 194
 195         TVector3 v2(nucleon2->Vx(), nucleon2->Vy(), nucleon2->Vz());
 196         TVector3 p2(nucleon2->Px(), nucleon2->Py(), nucleon2->Pz());
 197
 198         Double_t deltaP = this->GetPcm(p1, kProtonMass, p2, kNeutronMass); // relative momentum in CM frame
 199         if( deltaP >= fPmax) return -1.;
 200
 201         Double_t deltaR = (v2-v1).Mag();       // relative distance (cm)
 202         if(deltaR >= fRmax*1.e-13) return -1.;
 203
 204         if(Rndm() > fSpinProb)    return -1.;  // spin
 205
 206         if(fClusterType == kLowestMomentum) return 1. - deltaP/fPmax;
 207         if(fClusterType == kLowestDistance) return 1. - 1.e+13*deltaR/fRmax;
 208
 209         return 1. - 1.e+13*(deltaP*deltaR)/(fRmax*fPmax);
 210 }
 211
 212 void AliGenDeuteron::FirstPartner(const TList* protons, TList* neutrons)
 213 {
 214 //
 215 // Clusters are made with the first nucleon pair that fulfill
 216 // the coalescence conditions, starting with the protons
 217 //
 218         TIter p_next(protons);
 219         while(TParticle* n0 = (TParticle*) p_next())
 220         {
 221                 TParticle* partner = 0;
 222                 TIter n_next(neutrons);
 223                 while(TParticle* n1 = (TParticle*) n_next() )
 224                 {
 225                         if(GetCoalescenceProbability(n0, n1) < 0 ) continue; // with next neutron
 226                         partner = n1;
 227                         break;
 228                 }
 229
 230                 if(partner == 0) continue; // with next proton
 231
 232                 PushDeuteron(n0, partner);
 233
 234                 // Remove from the list for the next iteration
 235                 neutrons->Remove(partner);
 236         }
 237 }
 238
 239 void AliGenDeuteron::WeightMatrix(const TList* protons, const TList* neutrons)
 240 {
 241 //
 242 // Build all possible nucleon pairs with their own probability
 243 // and select only those with the highest coalescence probability
 244 //
 245         Int_t nMaxPairs = protons->GetSize()*neutrons->GetSize();
 246
 247         TParticle** cProton = new TParticle*[nMaxPairs];
 248         TParticle** cNeutron = new TParticle*[nMaxPairs];
 249         Double_t*  cWeight = new Double_t[nMaxPairs];
 250
 251         // build all pairs with probability > 0
 252         Int_t cIdx = -1;
 253         TIter p_next(protons);
 254         while(TParticle* n1 = (TParticle*) p_next())
 255         {
 256                 TIter n_next(neutrons);
 257                 while(TParticle* n2 = (TParticle*) n_next() )
 258                 {
 259                         Double_t weight = this->GetCoalescenceProbability(n1,n2);
 260                         if(weight < 0) continue;
 261                         ++cIdx;
 262                         cProton[cIdx]  = n1;
 263                         cNeutron[cIdx] = n2;
 264                         cWeight[cIdx]  = weight;
 265                 }
 266                 n_next.Reset();
 267         }
 268         p_next.Reset();
 269
 270         Int_t nPairs = cIdx + 1;
 271
 272         // find the interacting pairs:
 273         // remove repeated pairs and select only
 274         // the pair with the highest coalescence probability
 275
 276         Int_t nMaxIntPair = TMath::Min(protons->GetSize(), neutrons->GetSize());
 277
 278         TParticle** iProton = new TParticle*[nMaxIntPair];
 279         TParticle** iNeutron = new TParticle*[nMaxIntPair];
 280         Double_t* iWeight = new Double_t[nMaxIntPair];
 281
 282         Int_t iIdx = -1;
 283         while(kTRUE)
 284         {
 285                 Int_t j = -1;
 286                 Double_t wMax = 0;
 287                 for(Int_t i=0; i < nPairs; ++i)
 288                 {
 289                         if(cWeight[i] > wMax)
 290                         {
 291                                 wMax=cWeight[i];
 292                                 j = i;
 293                         }
 294                 }
 295
 296                 if(j == -1 ) break; // end
 297
 298                 // Save the interacting pair
 299                 ++iIdx;
 300                 iProton[iIdx]  = cProton[j];
 301                 iNeutron[iIdx] = cNeutron[j];
 302                 iWeight[iIdx]  = cWeight[j];
 303
 304                 // invalidate all combinations with these pairs for the next iteration
 305                 for(Int_t i=0; i < nPairs; ++i)
 306                 {
 307                         if(cProton[i] == iProton[iIdx]) cWeight[i] = -1.;
 308                         if(cNeutron[i] == iNeutron[iIdx]) cWeight[i] = -1.;
 309                 }
 310         }
 311
 312         Int_t nIntPairs = iIdx + 1;
 313
 314         delete[] cProton;
 315         delete[] cNeutron;
 316         delete[] cWeight;
 317
 318         // Add the (anti)deuterons to the current event stack
 319         for(Int_t i=0; i<nIntPairs; ++i)
 320         {
 321                 TParticle* n1 = iProton[i];
 322                 TParticle* n2 = iNeutron[i];
 323                 PushDeuteron(n1,n2);
 324         }
 325
 326         delete[] iProton;
 327         delete[] iNeutron;
 328         delete[] iWeight;
 329 }
 330
 331 void AliGenDeuteron::PushDeuteron(TParticle* parent1, TParticle* parent2)
 332 {
 333 //
 334 // Create an (anti)deuteron from parent1 and parent2,
 335 // add to the current stack and set kDoneBit for the parents
 336 //
 337         const Double_t kDeuteronMass = 1.87561282;
 338         const Int_t kDeuteronPdg = 1000010020;
 339
 340         // momentum
 341         TVector3 p1(parent1->Px(), parent1->Py(), parent1->Pz());
 342         TVector3 p2(parent2->Px(), parent2->Py(), parent2->Pz());
 343         TVector3 pN = p1+p2;
 344
 345         // production vertex same as the parent1's
 346         TVector3 vN(parent1->Vx(), parent1->Vy(), parent1->Vz());
 347
 348         // E^2 = p^2 + m^2
 349         Double_t energy = TMath::Sqrt(pN.Mag2() + kDeuteronMass*kDeuteronMass);
 350
 351         Int_t ntrk = 0;
 352         Double_t weight = 1;
 353         Int_t is = 1; // final state particle
 354
 355         // Add a new (anti)deuteron to current event stack
 356         fCurStack->PushTrack(1, -1, fSign*kDeuteronPdg,
 357                          pN.X(), pN.Y(), pN.Z(), energy,
 358                          vN.X(), vN.Y(), vN.Z(), parent1->T(),
 359                          0., 0., 0., kPNCapture, ntrk, weight, is);
 360
 361         // change the status code of the parents
 362         parent1->SetStatusCode(kCluster);
 363         parent2->SetStatusCode(kCluster);
 364
 365         // Set kDoneBit for the parents
 366         parent1->SetBit(kDoneBit);
 367         parent2->SetBit(kDoneBit);
 368 }
 369
 370 void AliGenDeuteron::FixProductionVertex(TParticle* i)
 371 {
 372 //
 373 // Replace current generator nucleon spatial distribution
 374 // with a custom distribution according to the selected model
 375 //
 376         if(fModel == kNone || fModel > kExpansion) return;
 377
 378         // semi-axis from collision geometry (fm)
 379         Double_t a = fTimeLength + TMath::Sqrt(fR*fR - fB*fB/4.);
 380         Double_t b = fTimeLength + fR - fB/2.;
 381         Double_t c = fTimeLength;
 382
 383         Double_t xx = 0;
 384         Double_t yy = 0;
 385         Double_t zz = 0;
 386
 387         if(fModel == kThermal)
 388         {
 389                 // uniformly ditributed in the volume on an ellipsoid
 390                 // random (r,theta,phi) unit sphere
 391                 Double_t r = TMath::Power(Rndm(),1./3.);
 392                 Double_t theta = TMath::ACos(2.*Rndm()-1.);
 393                 Double_t phi = 2.*TMath::Pi()*Rndm();
 394
 395                 // transform coordenates
 396                 xx = a*r*TMath::Sin(theta)*TMath::Cos(phi);
 397                 yy = b*r*TMath::Sin(theta)*TMath::Sin(phi);
 398                 zz = c*r*TMath::Cos(theta);
 399         }
 400         else if(fModel == kExpansion)
 401         {
 402                 // project into the surface of an ellipsoid
 403                 xx = a*TMath::Sin(i->Theta())*TMath::Cos(i->Phi());
 404                 yy = b*TMath::Sin(i->Theta())*TMath::Sin(i->Phi());
 405                 zz = c*TMath::Cos(i->Theta());
 406         }
 407
 408         // rotate by the reaction plane angle
 409         Double_t x = xx*TMath::Cos(fPsiR)+yy*TMath::Sin(fPsiR);
 410         Double_t y = -xx*TMath::Sin(fPsiR)+yy*TMath::Cos(fPsiR);
 411         Double_t z = zz;
 412
 413         // translate by the production vertex (cm)
 414         i->SetProductionVertex(i->Vx() + 1.e-13*x, i->Vy() + 1.e-13*y, i->Vz() + 1.e-13*z, i->T());
 415 }
 416
 417 Double_t AliGenDeuteron::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 AliGenDeuteron::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
 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 }
 437
 438 Double_t AliGenDeuteron::GetPcm(const TVector3& p1, Double_t m1, const TVector3& p2, Double_t m2) const
 439 {
 440 //
 441 // momentum in the CM frame
 442 //
 443         return this->GetPcm(p1.Px(),p1.Py(),p1.Pz(),m1,p2.Px(),p2.Py(),p2.Pz(),m2);
 444 }
 445