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  ,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                 fCurStack = gener->GetStack(ns);
 130                 if(!fCurStack)
 131                 {
 132                         AliWarning("no event stack");
 133                         return;
 134                 }
 135
 136                 TList* protons = new TList();
 137                 protons->SetOwner(kFALSE);
 138                 TList* neutrons = new TList();
 139                 neutrons->SetOwner(kFALSE);
 140
 141                 // particles produced by the generator
 142                 for (Int_t i=0; i < fCurStack->GetNprimary(); ++i)
 143                 {
 144                         TParticle* iParticle = fCurStack->Particle(i);
 145
 146                         if(iParticle->GetStatusCode() != 1) continue;
 147
 148                         Int_t pdgCode = iParticle->GetPdgCode();
 149                         if(pdgCode == fSign*2212)// (anti)proton
 150                         {
 151                                 FixProductionVertex(iParticle);
 152                                 protons->Add(iParticle);
 153                         }
 154                         else if(pdgCode == fSign*2112) // (anti)neutron
 155                         {
 156                                 FixProductionVertex(iParticle);
 157                                 neutrons->Add(iParticle);
 158                         }
 159                 }
 160
 161                 // look for clusters
 162                 if(fClusterType==kFirstPartner)
 163                 {
 164                         FirstPartner(protons, neutrons);
 165                 }
 166                 else
 167                 {
 168                         WeightMatrix(protons, neutrons);
 169                 }
 170
 171                 protons->Clear("nodelete");
 172                 neutrons->Clear("nodelete");
 173
 174                 delete protons;
 175                 delete neutrons;
 176         }
 177
 178         AliInfo("DONE");
 179 }
 180
 181 Double_t AliGenDeuteron::GetCoalescenceProbability(const TParticle* nucleon1, const TParticle* nucleon2) const
 182 {
 183 //
 184 // Coalescence conditions as in
 185 // A. J. Baltz et al., Phys. lett B 325(1994)7
 186 //
 187 // returns < 0 if coalescence is not possible
 188 // otherwise returns a coalescence probability
 189 //
 190         const Double_t kProtonMass  = 0.938272013;
 191         const Double_t kNeutronMass = 0.939565378;
 192
 193         TVector3 v1(nucleon1->Vx(), nucleon1->Vy(), nucleon1->Vz());
 194         TVector3 p1(nucleon1->Px(), nucleon1->Py(), nucleon1->Pz());
 195
 196         TVector3 v2(nucleon2->Vx(), nucleon2->Vy(), nucleon2->Vz());
 197         TVector3 p2(nucleon2->Px(), nucleon2->Py(), nucleon2->Pz());
 198
 199         Double_t deltaP = 2.*this->GetPcm(p1, kProtonMass, p2, kNeutronMass); // relative momentum in CM frame
 200         if( deltaP >= fPmax) return -1.;
 201
 202         Double_t deltaR = (v2-v1).Mag();       // relative distance (cm)
 203         if(deltaR >= fRmax*1.e-13) return -1.;
 204
 205         if(Rndm() > fSpinProb)    return -1.;  // spin
 206
 207         if(fClusterType == kLowestMomentum) return 1. - deltaP/fPmax;
 208         if(fClusterType == kLowestDistance) return 1. - 1.e+13*deltaR/fRmax;
 209
 210         return 1. - 1.e+13*(deltaP*deltaR)/(fRmax*fPmax);
 211 }
 212
 213 void AliGenDeuteron::FirstPartner(const TList* protons, TList* neutrons)
 214 {
 215 //
 216 // Clusters are made with the first nucleon pair that fulfill
 217 // the coalescence conditions, starting with the protons
 218 //
 219         TIter p_next(protons);
 220         while(TParticle* n0 = (TParticle*) p_next())
 221         {
 222                 TParticle* partner = 0;
 223                 TIter n_next(neutrons);
 224                 while(TParticle* n1 = (TParticle*) n_next() )
 225                 {
 226                         if(GetCoalescenceProbability(n0, n1) < 0 ) continue; // with next neutron
 227                         partner = n1;
 228                         break;
 229                 }
 230
 231                 if(partner == 0) continue; // with next proton
 232
 233                 PushDeuteron(n0, partner);
 234
 235                 // Remove from the list for the next iteration
 236                 neutrons->Remove(partner);
 237         }
 238 }
 239
 240 void AliGenDeuteron::WeightMatrix(const TList* protons, const TList* neutrons)
 241 {
 242 //
 243 // Build all possible nucleon pairs with their own probability
 244 // and select only those with the highest coalescence probability
 245 //
 246         Int_t nMaxPairs = protons->GetSize()*neutrons->GetSize();
 247
 248         TParticle** cProton = new TParticle*[nMaxPairs];
 249         TParticle** cNeutron = new TParticle*[nMaxPairs];
 250         Double_t*  cWeight = new Double_t[nMaxPairs];
 251
 252         // build all pairs with probability > 0
 253         Int_t cIdx = -1;
 254         TIter p_next(protons);
 255         while(TParticle* n1 = (TParticle*) p_next())
 256         {
 257                 TIter n_next(neutrons);
 258                 while(TParticle* n2 = (TParticle*) n_next() )
 259                 {
 260                         Double_t weight = this->GetCoalescenceProbability(n1,n2);
 261                         if(weight < 0) continue;
 262                         ++cIdx;
 263                         cProton[cIdx]  = n1;
 264                         cNeutron[cIdx] = n2;
 265                         cWeight[cIdx]  = weight;
 266                 }
 267                 n_next.Reset();
 268         }
 269         p_next.Reset();
 270
 271         Int_t nPairs = cIdx + 1;
 272
 273         // find the interacting pairs:
 274         // remove repeated pairs and select only
 275         // the pair with the highest coalescence probability
 276
 277         Int_t nMaxIntPair = TMath::Min(protons->GetSize(), neutrons->GetSize());
 278
 279         TParticle** iProton = new TParticle*[nMaxIntPair];
 280         TParticle** iNeutron = new TParticle*[nMaxIntPair];
 281         Double_t* iWeight = new Double_t[nMaxIntPair];
 282
 283         Int_t iIdx = -1;
 284         while(kTRUE)
 285         {
 286                 Int_t j = -1;
 287                 Double_t wMax = 0;
 288                 for(Int_t i=0; i < nPairs; ++i)
 289                 {
 290                         if(cWeight[i] > wMax)
 291                         {
 292                                 wMax=cWeight[i];
 293                                 j = i;
 294                         }
 295                 }
 296
 297                 if(j == -1 ) break; // end
 298
 299                 // Save the interacting pair
 300                 ++iIdx;
 301                 iProton[iIdx]  = cProton[j];
 302                 iNeutron[iIdx] = cNeutron[j];
 303                 iWeight[iIdx]  = cWeight[j];
 304
 305                 // invalidate all combinations with these pairs for the next iteration
 306                 for(Int_t i=0; i < nPairs; ++i)
 307                 {
 308                         if(cProton[i] == iProton[iIdx]) cWeight[i] = -1.;
 309                         if(cNeutron[i] == iNeutron[iIdx]) cWeight[i] = -1.;
 310                 }
 311         }
 312
 313         Int_t nIntPairs = iIdx + 1;
 314
 315         delete[] cProton;
 316         delete[] cNeutron;
 317         delete[] cWeight;
 318
 319         // Add the (anti)deuterons to the current event stack
 320         for(Int_t i=0; i<nIntPairs; ++i)
 321         {
 322                 TParticle* n1 = iProton[i];
 323                 TParticle* n2 = iNeutron[i];
 324                 PushDeuteron(n1,n2);
 325         }
 326
 327         delete[] iProton;
 328         delete[] iNeutron;
 329         delete[] iWeight;
 330 }
 331
 332 void AliGenDeuteron::PushDeuteron(TParticle* parent1, TParticle* parent2)
 333 {
 334 //
 335 // Create an (anti)deuteron from parent1 and parent2,
 336 // add to the current stack and set kDoneBit for the parents
 337 //
 338         const Double_t kDeuteronMass = 1.87561282;
 339         const Int_t kDeuteronPdg = 1000010020;
 340
 341         // momentum
 342         TVector3 p1(parent1->Px(), parent1->Py(), parent1->Pz());
 343         TVector3 p2(parent2->Px(), parent2->Py(), parent2->Pz());
 344         TVector3 pN = p1+p2;
 345
 346         // production vertex same as the parent1's
 347         TVector3 vN(parent1->Vx(), parent1->Vy(), parent1->Vz());
 348
 349         // E^2 = p^2 + m^2
 350         Double_t energy = TMath::Sqrt(pN.Mag2() + kDeuteronMass*kDeuteronMass);
 351
 352         // Add a new (anti)deuteron to current event stack
 353         fCurStack->PushTrack(1, -1, fSign*kDeuteronPdg,
 354                          pN.X(), pN.Y(), pN.Z(), energy,
 355                          vN.X(), vN.Y(), vN.Z(), parent1->T(),
 356                          0., 0., 0., kPNCapture, fNtrk, 1., 0);
 357
 358         // Set kDoneBit for the parents
 359         parent1->SetBit(kDoneBit);
 360         parent2->SetBit(kDoneBit);
 361 }
 362
 363 void AliGenDeuteron::FixProductionVertex(TParticle* i)
 364 {
 365 //
 366 // Replace current generator nucleon spatial distribution
 367 // with a custom distribution according to the selected model
 368 //
 369         if(fModel == kNone || fModel > kExpansion) return;
 370
 371         // semi-axis from collision geometry (fm)
 372         Double_t a = fTimeLength + TMath::Sqrt(fR*fR - fB*fB/4.);
 373         Double_t b = fTimeLength + fR - fB/2.;
 374         Double_t c = fTimeLength;
 375
 376         Double_t xx = 0;
 377         Double_t yy = 0;
 378         Double_t zz = 0;
 379
 380         if(fModel == kThermal)
 381         {
 382                 // uniformly ditributed in the volume on an ellipsoid
 383                 // random (r,theta,phi) unit sphere
 384                 Double_t r = TMath::Power(Rndm(),1./3.);
 385                 Double_t theta = TMath::ACos(2.*Rndm()-1.);
 386                 Double_t phi = 2.*TMath::Pi()*Rndm();
 387
 388                 // transform coordenates
 389                 xx = a*r*TMath::Sin(theta)*TMath::Cos(phi);
 390                 yy = b*r*TMath::Sin(theta)*TMath::Sin(phi);
 391                 zz = c*r*TMath::Cos(theta);
 392         }
 393         else if(fModel == kExpansion)
 394         {
 395                 // project into the surface of an ellipsoid
 396                 xx = a*TMath::Sin(i->Theta())*TMath::Cos(i->Phi());
 397                 yy = b*TMath::Sin(i->Theta())*TMath::Sin(i->Phi());
 398                 zz = c*TMath::Cos(i->Theta());
 399         }
 400
 401         // rotate by the reaction plane angle
 402         Double_t x = xx*TMath::Cos(fPsiR)+yy*TMath::Sin(fPsiR);
 403         Double_t y = -xx*TMath::Sin(fPsiR)+yy*TMath::Cos(fPsiR);
 404         Double_t z = zz;
 405
 406         // translate by the production vertex (cm)
 407         i->SetProductionVertex(i->Vx() + 1.e-13*x, i->Vy() + 1.e-13*y, i->Vz() + 1.e-13*z, i->T());
 408 }
 409
 410 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
 411 {
 412 //
 413 // square of the center of mass energy
 414 //
 415         Double_t E1 = TMath::Sqrt( p1x*p1x + p1y*p1y + p1z*p1z + m1*m1);
 416         Double_t E2 = TMath::Sqrt( p2x*p2x + p2y*p2y + p2z*p2z + m2*m2);
 417
 418         return (E1+E2)*(E1+E2) - ((p1x+p2x)*(p1x+p2x) + (p1y+p2y)*(p1y+p2y) + (p1z+p2z)*(p1z+p2z));
 419 }
 420
 421 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
 422 {
 423 //
 424 // momentum in the CM frame
 425 //
 426         Double_t s = this->GetS(p1x, p1y, p1z, m1, p2x, p2y, p2z, m2);
 427
 428         return TMath::Sqrt( (s-(m1-m2)*(m1-m2))*(s-(m1+m2)*(m1+m2)) )/(2.*TMath::Sqrt(s));
 429 }
 430
 431 Double_t AliGenDeuteron::GetPcm(const TVector3& p1, Double_t m1, const TVector3& p2, Double_t m2) const
 432 {
 433 //
 434 // momentum in the CM frame
 435 //
 436         return this->GetPcm(p1.Px(),p1.Py(),p1.Pz(),m1,p2.Px(),p2.Py(),p2.Pz(),m2);
 437 }
 438