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