PWG2/FLOW/AliFlowTools/glauberMC/AliGlauberNucleus.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, 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 //
  18 //  AliGlauberNucleus implementation
  19 //  support class for Glauber MC
  20 //
  21 //  origin: PHOBOS experiment
  22 //  alification: Mikolaj Krzewicki, Nikhef, mikolaj.krzewicki@cern.ch
  23 //
  24 ////////////////////////////////////////////////////////////////////////////////
  25
  26 #include <Riostream.h>
  27 #include <TMath.h>
  28 #include <TEllipse.h>
  29 #include <TNamed.h>
  30 #include <TObjArray.h>
  31 #include <TF1.h>
  32 #include <TRandom.h>
  33 #include "AliGlauberNucleon.h"
  34 #include "AliGlauberNucleus.h"
  35
  36 ClassImp(AliGlauberNucleus)
  37
  38 //______________________________________________________________________________
  39 AliGlauberNucleus::AliGlauberNucleus(Option_t* iname, Int_t iN, Double_t iR, Double_t ia, Double_t iw, TF1* ifunc) :
  40    fN(iN),fR(iR),fA(ia),fW(iw),fMinDist(-1),
  41    fF(0),fTrials(0),fFunction(ifunc),
  42    fNucleons(0)
  43 {
  44    if (fN==0) {
  45       cout << "Setting up nucleus " << iname << endl;
  46       Lookup(iname);
  47    }
  48 }
  49
  50 //______________________________________________________________________________
  51 AliGlauberNucleus::~AliGlauberNucleus()
  52 {
  53    if (fNucleons) {
  54       delete fNucleons;
  55    }
  56    delete fFunction;
  57 }
  58
  59 //______________________________________________________________________________
  60 void AliGlauberNucleus::Draw(Double_t xs, Int_t col)
  61 {
  62    Double_t r = 0.5*sqrt(xs/TMath::Pi()/10.);
  63    TEllipse e;
  64    e.SetLineColor(col);
  65    e.SetFillColor(0);
  66    e.SetLineWidth(1);
  67
  68    for (Int_t i = 0;i<fNucleons->GetEntries();++i) {
  69       AliGlauberNucleon* gn = (AliGlauberNucleon*) fNucleons->At(i);
  70       e.SetLineStyle(1);
  71       if (gn->IsSpectator()) e.SetLineStyle(3);
  72       e.DrawEllipse(gn->GetX(),gn->GetY(),r,r,0,360,0,"");
  73    }
  74 }
  75
  76 //______________________________________________________________________________
  77 void AliGlauberNucleus::Lookup(Option_t* name)
  78 {
  79    SetName(name);
  80
  81    if      (TString(name) == "p")    {fN = 1;   fR = 0.6;   fA = 0;      fW =  0;      fF = 0;}
  82    else if (TString(name) == "d")    {fN = 2;   fR = 0.01;  fA = 0.5882; fW =  0;      fF = 1;}
  83    else if (TString(name) == "dh")   {fN = 2;   fR = 0.01;  fA = 0.5882; fW =  0;      fF = 3;}
  84    else if (TString(name) == "dhh")  {fN = 2;   fR = 0.01;  fA = 0.5882; fW =  0;      fF = 4;}
  85    else if (TString(name) == "O")    {fN = 16;  fR = 2.608; fA = 0.513;  fW = -0.051;  fF = 1;}
  86    else if (TString(name) == "Si")   {fN = 28;  fR = 3.34;  fA = 0.580;  fW = -0.233;  fF = 1;}
  87    else if (TString(name) == "S")    {fN = 32;  fR = 2.54;  fA = 2.191;  fW =  0.16;   fF = 2;}
  88    else if (TString(name) == "Ca")   {fN = 40;  fR = 3.766; fA = 0.586;  fW = -0.161;  fF = 1;}
  89    else if (TString(name) == "Ni")   {fN = 58;  fR = 4.309; fA = 0.517;  fW = -0.1308; fF = 1;}
  90    else if (TString(name) == "Cu")   {fN = 63;  fR = 4.2;   fA = 0.596;  fW =  0;      fF = 1;}
  91    else if (TString(name) == "W")    {fN = 186; fR = 6.58;  fA = 0.480;  fW =  0;      fF = 1;}
  92    else if (TString(name) == "Au")   {fN = 197; fR = 6.38;  fA = 0.535;  fW =  0;      fF = 1;}
  93    else if (TString(name) == "Pb")   {fN = 208; fR = 6.62;  fA = 0.546;  fW =  0;      fF = 1;}
  94    else if (TString(name) == "U")    {fN = 238; fR = 6.81;  fA = 0.6;    fW =  0;      fF = 1;}
  95    else {
  96       cout << "Could not find nucleus " << name << endl;
  97       return;
  98    }
  99
 100    switch (fF)
 101    {
 102       case 0: // Proton
 103          fFunction = new TF1("prot","x*x*exp(-x/[0])",0,10);
 104          fFunction->SetParameter(0,fR);
 105          break;
 106       case 1: // 3pF
 107          fFunction = new TF1(name,"x*x*(1+[2]*(x/[0])**2)/(1+exp((x-[0])/[1]))",0,15);
 108          fFunction->SetParameters(fR,fA,fW);
 109          break;
 110       case 2: // 3pG
 111          fFunction = new TF1("3pg","x*x*(1+[2]*(x/[0])**2)/(1+exp((x**2-[0]**2)/[1]**2))",0,15);
 112          fFunction->SetParameters(fR,fA,fW);
 113          break;
 114       case 3: // Hulthen
 115          fFunction = new TF1("f3","x*x*([0]*[1]*([0]+[1]))/(2*pi*(pow([0]-[1],2)))*pow((exp(-[0]*x)-exp(-[1]*x))/x,2)",0,10);
 116          fFunction->SetParameters(1/4.38,1/.85);
 117          break;
 118       case 4: // Hulthen HIJING
 119          fFunction = new TF1("f4","x*x*([0]*[1]*([0]+[1]))/(2*pi*(pow([0]-[1],2)))*pow((exp(-[0]*x)-exp(-[1]*x))/x,2)",0,20);
 120          fFunction->SetParameters(2/4.38,2/.85);
 121          break;
 122       default:
 123          cerr << "Could not find function type " << fF << endl;
 124          return;
 125    }
 126 }
 127
 128 //______________________________________________________________________________
 129 void AliGlauberNucleus::SetR(Double_t ir)
 130 {
 131    fR = ir;
 132    switch (fF)
 133    {
 134       case 0: // Proton
 135          fFunction->SetParameter(0,fR);
 136          break;
 137       case 1: // 3pF
 138          fFunction->SetParameter(0,fR);
 139          break;
 140       case 2: // 3pG
 141          fFunction->SetParameter(0,fR);
 142          break;
 143    }
 144 }
 145
 146 //______________________________________________________________________________
 147 void AliGlauberNucleus::SetA(Double_t ia)
 148 {
 149    fA = ia;
 150    switch (fF)
 151    {
 152       case 0: // Proton
 153          break;
 154       case 1: // 3pF
 155          fFunction->SetParameter(1,fA);
 156          break;
 157       case 2: // 3pG
 158          fFunction->SetParameter(1,fA);
 159          break;
 160    }
 161 }
 162
 163 //______________________________________________________________________________
 164 void AliGlauberNucleus::SetW(Double_t iw)
 165 {
 166    fW = iw;
 167    switch (fF)
 168    {
 169       case 0: // Proton
 170          break;
 171       case 1: // 3pF
 172          fFunction->SetParameter(2,fW);
 173          break;
 174       case 2: // 3pG
 175          fFunction->SetParameter(2,fW);
 176          break;
 177    }
 178 }
 179
 180 //______________________________________________________________________________
 181 void AliGlauberNucleus::ThrowNucleons(Double_t xshift)
 182 {
 183    if (fNucleons==0) {
 184       fNucleons=new TObjArray(fN);
 185       fNucleons->SetOwner();
 186       for(Int_t i=0;i<fN;i++) {
 187          AliGlauberNucleon *nucleon=new AliGlauberNucleon();
 188          fNucleons->Add(nucleon);
 189       }
 190    }
 191
 192    fTrials = 0;
 193
 194    Double_t sumx=0;
 195    Double_t sumy=0;
 196    Double_t sumz=0;
 197
 198    Bool_t hulthen = (TString(GetName())=="dh");
 199    if (fN==2 && hulthen) { //special treatmeant for Hulten
 200
 201       Double_t r = fFunction->GetRandom()/2;
 202       Double_t phi = gRandom->Rndm() * 2 * TMath::Pi() ;
 203       Double_t ctheta = 2*gRandom->Rndm() - 1 ;
 204       Double_t stheta = sqrt(1-ctheta*ctheta);
 205
 206       AliGlauberNucleon *nucleon1=(AliGlauberNucleon*)(fNucleons->At(0));
 207       AliGlauberNucleon *nucleon2=(AliGlauberNucleon*)(fNucleons->At(1));
 208       nucleon1->Reset();
 209       nucleon1->SetXYZ(r * stheta * cos(phi) + xshift,
 210                        r * stheta * sin(phi),
 211                        r * ctheta);
 212       nucleon2->Reset();
 213       nucleon2->SetXYZ(-nucleon1->GetX() + 2*xshift,
 214                        -nucleon1->GetY(),
 215                        -nucleon1->GetZ());
 216       fTrials = 1;
 217       return;
 218    }
 219
 220    for (Int_t i = 0; i<fN; i++) {
 221       AliGlauberNucleon *nucleon=(AliGlauberNucleon*)(fNucleons->At(i));
 222       nucleon->Reset();
 223       while(1) {
 224          fTrials++;
 225          Double_t r = fFunction->GetRandom();
 226          Double_t phi = gRandom->Rndm() * 2 * TMath::Pi() ;
 227          Double_t ctheta = 2*gRandom->Rndm() - 1 ;
 228          Double_t stheta = TMath::Sqrt(1-ctheta*ctheta);
 229          Double_t x = r * stheta * cos(phi) + xshift;
 230          Double_t y = r * stheta * sin(phi);
 231          Double_t z = r * ctheta;
 232          nucleon->SetXYZ(x,y,z);
 233          if(fMinDist<0) break;
 234          Bool_t test=1;
 235          for (Int_t j = 0; j<i; j++) {
 236             AliGlauberNucleon *other=(AliGlauberNucleon*)fNucleons->At(j);
 237             Double_t xo=other->GetX();
 238             Double_t yo=other->GetY();
 239             Double_t zo=other->GetZ();
 240             Double_t dist = TMath::Sqrt((x-xo)*(x-xo)+
 241                                        (y-yo)*(y-yo)+
 242                                        (z-zo)*(z-zo));
 243
 244             if(dist<fMinDist) {
 245                test=0;
 246                break;
 247             }
 248          }
 249          if (test) break; //found nucleuon outside of mindist
 250       }
 251
 252       sumx += nucleon->GetX();
 253       sumy += nucleon->GetY();
 254       sumz += nucleon->GetZ();
 255    }
 256
 257    if(1) { // set the centre-of-mass to be at zero (+xshift)
 258       sumx = sumx/fN;
 259       sumy = sumy/fN;
 260       sumz = sumz/fN;
 261       for (Int_t i = 0; i<fN; i++) {
 262          AliGlauberNucleon *nucleon=(AliGlauberNucleon*)(fNucleons->At(i));
 263          nucleon->SetXYZ(nucleon->GetX()-sumx-xshift,
 264                          nucleon->GetY()-sumy,
 265                          nucleon->GetZ()-sumz);
 266       }
 267    }
 268 }
 269