EVGEN/AliGenDeuteron.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 2009-2010, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 \r
  16 //\r
  17 // A very simple model based on the article "Physics Letters B325 (1994) 7-12".\r
  18 // The only addition is to model the freeze-out at which the\r
  19 // light nuclei can form for the energies of the LHC, since PYTHIA does not give\r
  20 // any detailed spatial description of the generated particles at the level of a\r
  21 // few fermis.\r
  22 //\r
  23 // There are two options to place the nucleons: a thermal picture where all\r
  24 // nucleons are placed randomly and homogeneously in an spherical volume of a\r
  25 // few fermis, and expansion one where they are projected on its surface.\r
  26 //\r
  27 // A (anti)deuteron will form if there is a pair of (anti)proton-(anti)neutron\r
  28 // with momentum difference less than ~ 300MeV and relative distance less than a\r
  29 // ~ 2.1fm. Only 3/4 of this clusters are accepted due to the deuteron spin.\r
  30 //\r
  31 // Author: Eulogio Serradilla <eulogio.serradilla@ciemat.es>,\r
  32 //         Arturo Menchaca <menchaca@fisica.unam.mx>\r
  33 //\r
  34 \r
  35 #include "Riostream.h"\r
  36 #include "TParticle.h"\r
  37 #include "AliStack.h"\r
  38 #include "AliGenDeuteron.h"\r
  39 #include "AliGenCocktailAfterBurner.h"\r
  40 #include "TMath.h"\r
  41 #include "TMCProcess.h"\r
  42 #include "TList.h"\r
  43 #include "TVector3.h"\r
  44 #include "AliMC.h"\r
  45 #include "TArrayF.h"\r
  46 #include "AliGenCocktailEventHeader.h"\r
  47 \r
  48 ClassImp(AliGenDeuteron)\r
  49 \r
  50 AliGenDeuteron::AliGenDeuteron(Int_t sign, Double_t pmax, Double_t rmax, Double_t rsrc, Int_t model):\r
  51 fDeuteronMass(1.87561282), //pdg\r
  52 fPmax(pmax), // GeV/c\r
  53 fRmax(rmax*1.e-13), //cm\r
  54 fRsrc(rsrc*1.e-13), // cm\r
  55 fSpinProb(0.75),\r
  56 fModel(model),\r
  57 fSign(1)\r
  58 {\r
  59         fSign = sign > 0 ? 1:-1;\r
  60 }\r
  61 \r
  62 AliGenDeuteron::~AliGenDeuteron()\r
  63 {\r
  64 }\r
  65 \r
  66 void AliGenDeuteron::Init()\r
  67 {\r
  68         // Standard AliGenerator Initializer\r
  69 }\r
  70 \r
  71 void AliGenDeuteron::Generate()\r
  72 {\r
  73 // Modify the stack of each event, adding the new particles at the end and\r
  74 // not transporting the parents.\r
  75 \r
  76         Info("Generate","freeze-out model : %d (0 expansion, 1 thermal)",fModel);\r
  77         Info("Generate","relative momentum: %g GeV/c",fPmax);\r
  78         Info("Generate","relative distance: %g cm",fRmax);\r
  79         Info("Generate","source radius    : %g cm",fRsrc);\r
  80         Info("Generate","spin probability : %g ",fSpinProb);\r
  81         Info("Generate","sign             : %d ",fSign);\r
  82         \r
  83         Int_t ntr;\r
  84         \r
  85         // Get the cocktail generator\r
  86         AliGenCocktailAfterBurner* gener = (AliGenCocktailAfterBurner*)gAlice->GetMCApp()->Generator();\r
  87         \r
  88         // Loop over events\r
  89         for(Int_t ns = 0; ns < gener->GetNumberOfEvents(); ++ns)\r
  90         {\r
  91                 gener->SetActiveEventNumber(ns);\r
  92                 \r
  93                 AliStack* stack = gener->GetStack(ns); // Stack of event ns\r
  94                 if(!stack)\r
  95                 {\r
  96                         Info("Generate", "no stack, exiting");\r
  97                         return;\r
  98                 }\r
  99                 \r
 100                 // find emerging nucleons from the collision of current event\r
 101                 \r
 102                 TList* protons = new TList();\r
 103                 protons->SetOwner(kFALSE);\r
 104                 TList* neutrons = new TList();\r
 105                 neutrons->SetOwner(kFALSE);\r
 106                 \r
 107                 // primary vertex of current event\r
 108                 TArrayF primVtx;\r
 109                 gener->GetActiveEventHeader()->PrimaryVertex(primVtx);\r
 110                 TVector3 r0(primVtx[0],primVtx[1],primVtx[2]);\r
 111                 \r
 112                 for (Int_t i=0; i<stack->GetNtrack(); ++i)\r
 113                 {\r
 114                         TParticle* iParticle = stack->Particle(i);\r
 115                         \r
 116                         if(iParticle->TestBit(kDoneBit)) continue;\r
 117                         // select only nucleons within the freeze-out volume\r
 118                         TVector3 r(iParticle->Vx(),iParticle->Vy(),iParticle->Vz());\r
 119                         if((r-r0).Mag()>fRsrc) continue;\r
 120                         \r
 121                         Int_t pdgCode = iParticle->GetPdgCode();\r
 122                         if(pdgCode == fSign*2212)// (anti)proton\r
 123                         {\r
 124                                 FixProductionVertex(iParticle);\r
 125                                 protons->Add(iParticle);\r
 126                         }\r
 127                         else if(pdgCode == fSign*2112) // (anti)neutron\r
 128                         {\r
 129                                 FixProductionVertex(iParticle);\r
 130                                 neutrons->Add(iParticle);\r
 131                         }\r
 132                 }\r
 133                 \r
 134                 // coalescence conditions\r
 135                 // A.J. Baltz et al., Phys. lett B 325(1994)7\r
 136                 \r
 137                 TIter p_next(protons);\r
 138                 while(TParticle* n0 = (TParticle*) p_next())\r
 139                 {\r
 140                         TParticle* iProton = 0;\r
 141                         TParticle* jNeutron = 0;\r
 142                         \r
 143                         // Select the first nucleon\r
 144                         iProton = n0;\r
 145                         \r
 146                         TVector3 v0(n0->Vx(), n0->Vy(), n0->Vz());\r
 147                         TVector3 p0(n0->Px(), n0->Py(), n0->Pz()), p1(0,0,0);\r
 148                         \r
 149                         // See if there is a neibourgh...\r
 150                         TIter n_next(neutrons);\r
 151                         while(TParticle* n1 = (TParticle*) n_next() )\r
 152                         {\r
 153                                 TVector3 v(n1->Vx(), n1->Vy(), n1->Vz());\r
 154                                 TVector3 p(n1->Px(), n1->Py(), n1->Pz());\r
 155                                 \r
 156                                 if((p-p0).Mag() > fPmax) continue;\r
 157                                 if((v-v0).Mag() > fRmax) continue;\r
 158                                 \r
 159                                 jNeutron = n1;\r
 160                                 p1 = p;\r
 161                                 break;\r
 162                         }\r
 163                         \r
 164                         if(jNeutron == 0) continue; // with next proton\r
 165                         \r
 166                         if(Rndm() > fSpinProb) continue;\r
 167                         \r
 168                         // neutron captured!\r
 169                         \r
 170                         TVector3 pN = p0+p1;\r
 171                         TVector3 vN(iProton->Vx(), iProton->Vy(), iProton->Vz()); // production vertex = p = n = collision vertex\r
 172                         \r
 173                         // E^2 = p^2 + m^2\r
 174                         Double_t EN = TMath::Sqrt(pN.Mag2()+fDeuteronMass*fDeuteronMass);\r
 175                         \r
 176                         // Add a new (anti)deuteron to the current event stack\r
 177                         stack->PushTrack(1, stack->Particles()->IndexOf(iProton), fSign*1000010020,\r
 178                                          pN.X(), pN.Y(), pN.Z(), EN,\r
 179                                          vN.X(), vN.Y(), vN.Z(), iProton->T(),\r
 180                                          0., 0., 0., kPNCapture, ntr,1.,0);\r
 181                         \r
 182                         //Info("Generate","neutron capture (NEW DEUTERON)");\r
 183                         \r
 184                         // Set bit not to transport for the nucleons\r
 185                         iProton->SetBit(kDoneBit);\r
 186                         jNeutron->SetBit(kDoneBit);\r
 187                         \r
 188                         // Remove from the list\r
 189                         neutrons->Remove(jNeutron);\r
 190                 }\r
 191                 \r
 192                 protons->Clear("nodelete");\r
 193                 neutrons->Clear("nodelete");\r
 194                 \r
 195                 delete protons;\r
 196                 delete neutrons;\r
 197         }\r
 198         \r
 199         Info("Generate","Coalescence afterburner: DONE");\r
 200 }\r
 201 \r
 202 // create the freeze-out nucleon distribution around the collision vertex\r
 203 void AliGenDeuteron::FixProductionVertex(TParticle* i)\r
 204 {\r
 205         Double_t x,y,z;\r
 206         \r
 207         if(fModel == kThermal) // homogeneous volume\r
 208         {\r
 209                 // random (r,theta,phi)\r
 210                 Double_t r = fRsrc*TMath::Power(Rndm(),1./3.);\r
 211                 Double_t theta = TMath::ACos(2.*Rndm()-1.);\r
 212                 Double_t phi = 2.*TMath::Pi()*Rndm();\r
 213         \r
 214                 // transform coordenates\r
 215                 x = r*TMath::Sin(theta)*TMath::Cos(phi);\r
 216                 y = r*TMath::Sin(theta)*TMath::Sin(phi);\r
 217                 z = r*TMath::Cos(theta);\r
 218         }\r
 219         else // projection into the surface of an sphere\r
 220         {\r
 221                 x = fRsrc*TMath::Sin(i->Theta())*TMath::Cos(i->Phi());\r
 222                 y = fRsrc*TMath::Sin(i->Theta())*TMath::Sin(i->Phi());\r
 223                 z = fRsrc*TMath::Cos(i->Theta());\r
 224         }\r
 225         \r
 226         // assume nucleons at the freeze-out have the same production vertex\r
 227         i->SetProductionVertex(i->Vx()+x, i->Vy()+y, i->Vz()+z, i->T());\r
 228 }\r