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