PWG2/FLOW/Tools/glauberMC/AliGlauberMC.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 //  Glauber MC implementation
  19 //
  20 //  origin: PHOBOS experiment
  21 //  alification: Mikolaj Krzewicki, Nikhef, mikolaj.krzewicki@cern.ch
  22 //
  23 ////////////////////////////////////////////////////////////////////////////////
  24
  25 #include <Riostream.h>
  26 #include <TMath.h>
  27 #include <TEllipse.h>
  28 #include <TRandom.h>
  29 #include <TNamed.h>
  30 #include <TObjArray.h>
  31 #include <TNtuple.h>
  32 #include <TFile.h>
  33
  34 #include "AliGlauberNucleon.h"
  35 #include "AliGlauberNucleus.h"
  36 #include "AliGlauberMC.h"
  37
  38 ClassImp(AliGlauberMC)
  39
  40 //______________________________________________________________________________
  41 AliGlauberMC::AliGlauberMC(Option_t* NA, Option_t* NB, Double_t xsect) :
  42    fANucleus(NA),fBNucleus(NB),
  43    fXSect(0),fNucleonsA(0),fNucleonsB(0),fnt(0),
  44    fMeanX2(0),fMeanY2(0),fMeanXY(0),fMeanXParts(0),
  45    fMeanYParts(0),fMeanXSystem(0),fMeanYSystem(0),
  46    fMeanX_A(0),fMeanY_A(0),fMeanX_B(0),fMeanY_B(0),fB_MC(0),
  47    fEvents(0),fTotalEvents(0),fBMin(0),fBMax(0),fMaxNpartFound(0),
  48    fNpart(0),fNcoll(0),fSx2(0),fSy2(0),fSxy(0), fX(0.13), fNpp(8.)
  49 {
  50    fBMin = 0;
  51    fBMax = 20;
  52    fXSect = xsect;
  53
  54    TString name(Form("Glauber_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
  55    TString title(Form("Glauber %s+%s Version",fANucleus.GetName(),fBNucleus.GetName()));
  56    SetName(name);
  57    SetTitle(title);
  58 }
  59
  60 //______________________________________________________________________________
  61 Bool_t AliGlauberMC::CalcEvent(Double_t bgen)
  62 {
  63    // prepare event
  64    fANucleus.ThrowNucleons(-bgen/2.);
  65    fNucleonsA = fANucleus.GetNucleons();
  66    fAN = fANucleus.GetN();
  67    for (Int_t i = 0; i<fAN; i++) {
  68       AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->At(i));
  69       nucleonA->SetInNucleusA();
  70    }
  71    fBNucleus.ThrowNucleons(bgen/2.);
  72    fNucleonsB = fBNucleus.GetNucleons();
  73    fBN = fBNucleus.GetN();
  74    for (Int_t i = 0; i<fBN; i++) {
  75       AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->At(i));
  76       nucleonB->SetInNucleusB();
  77    }
  78
  79    // "ball" diameter = distance at which two balls interact
  80    Double_t d2 = (Double_t)fXSect/(TMath::Pi()*10); // in fm^2
  81
  82    // for each of the A nucleons in nucleus B
  83    for (Int_t i = 0; i<fBN; i++) {
  84       AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->At(i));
  85       for (Int_t j = 0 ; j < fAN ;j++) {
  86          AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->At(j));
  87          Double_t dx = nucleonB->GetX()-nucleonA->GetX();
  88          Double_t dy = nucleonB->GetY()-nucleonA->GetY();
  89          Double_t dij = dx*dx+dy*dy;
  90          if (dij < d2) {
  91             nucleonB->Collide();
  92             nucleonA->Collide();
  93          }
  94       }
  95   }
  96
  97   return CalcResults(bgen);
  98 }
  99
 100 //______________________________________________________________________________
 101 Bool_t AliGlauberMC::CalcResults(Double_t bgen)
 102 {
 103    // calc results for the given event
 104    fNpart=0;
 105    fNcoll=0;
 106    fMeanX2=0;
 107    fMeanY2=0;
 108    fMeanXY=0;
 109    fMeanXParts=0;
 110    fMeanYParts=0;
 111    fMeanXSystem=0;
 112    fMeanYSystem=0;
 113    fMeanX_A=0;
 114    fMeanY_A=0;
 115    fMeanX_B=0;
 116    fMeanY_B=0;
 117
 118    for (Int_t i = 0; i<fAN; i++) {
 119       AliGlauberNucleon *nucleonA=(AliGlauberNucleon*)(fNucleonsA->At(i));
 120       Double_t xA=nucleonA->GetX();
 121       Double_t yA=nucleonA->GetY();
 122       fMeanXSystem  += xA;
 123       fMeanYSystem  += yA;
 124       fMeanX_A  += xA;
 125       fMeanY_A  += yA;
 126
 127       if(nucleonA->IsWounded()) {
 128          fNpart++;
 129          fMeanXParts  += xA;
 130          fMeanYParts  += yA;
 131          fMeanX2 += xA * xA;
 132          fMeanY2 += yA * yA;
 133          fMeanXY += xA * yA;
 134       }
 135    }
 136
 137    for (Int_t i = 0; i<fBN; i++) {
 138       AliGlauberNucleon *nucleonB=(AliGlauberNucleon*)(fNucleonsB->At(i));
 139       Double_t xB=nucleonB->GetX();
 140       Double_t yB=nucleonB->GetY();
 141       fMeanXSystem  += xB;
 142       fMeanYSystem  += yB;
 143       fMeanX_B  += xB;
 144       fMeanY_B  += yB;
 145
 146       if(nucleonB->IsWounded()) {
 147          fNpart++;
 148          fMeanXParts  += xB;
 149          fMeanYParts  += yB;
 150          fMeanX2 += xB * xB;
 151          fMeanY2 += yB * yB;
 152          fMeanXY += xB * yB;
 153          fNcoll += nucleonB->GetNColl();
 154       }
 155    }
 156
 157    if (fNpart>0) {
 158       fMeanXParts /= fNpart;
 159       fMeanYParts /= fNpart;
 160       fMeanX2 /= fNpart;
 161       fMeanY2 /= fNpart;
 162       fMeanXY /= fNpart;
 163    } else {
 164       fMeanXParts = 0;
 165       fMeanYParts = 0;
 166       fMeanX2 = 0;
 167       fMeanY2 = 0;
 168       fMeanXY = 0;
 169    }
 170
 171    if(fAN+fBN>0) {
 172       fMeanXSystem /= (fAN + fBN);
 173       fMeanYSystem /= (fAN + fBN);
 174    } else {
 175       fMeanXSystem = 0;
 176       fMeanYSystem = 0;
 177    }
 178    if(fAN>0) {
 179       fMeanX_A /= fAN;
 180       fMeanY_A /= fAN;
 181    } else {
 182       fMeanX_A = 0;
 183       fMeanY_A = 0;
 184    }
 185
 186    if(fBN>0) {
 187       fMeanX_B /= fBN;
 188       fMeanY_B /= fBN;
 189    } else {
 190       fMeanX_B = 0;
 191       fMeanY_B = 0;
 192    }
 193
 194    fSx2=fMeanX2-(fMeanXParts*fMeanXParts);
 195    fSy2=fMeanY2-(fMeanYParts*fMeanYParts);
 196    fSxy=fMeanXY-fMeanXParts*fMeanYParts;
 197
 198    fB_MC = bgen;
 199
 200    fTotalEvents++;
 201    if (fNpart>0) fEvents++;
 202
 203    if (fNpart==0) return kFALSE;
 204    if (fNpart > fMaxNpartFound) fMaxNpartFound = fNpart;
 205
 206    return kTRUE;
 207 }
 208
 209 //______________________________________________________________________________
 210 void AliGlauberMC::Draw(Option_t* /*option*/)
 211 {
 212    fANucleus.Draw(fXSect, 2);
 213    fBNucleus.Draw(fXSect, 4);
 214
 215    TEllipse e;
 216    e.SetFillColor(0);
 217    e.SetLineColor(1);
 218    e.SetLineStyle(2);
 219    e.SetLineWidth(1);
 220    e.DrawEllipse(GetB()/2,0,fBNucleus.GetR(),fBNucleus.GetR(),0,360,0);
 221    e.DrawEllipse(-GetB()/2,0,fANucleus.GetR(),fANucleus.GetR(),0,360,0);
 222 }
 223
 224 //______________________________________________________________________________
 225 Double_t AliGlauberMC::GetTotXSect() const
 226 {
 227    return (1.*fEvents/fTotalEvents)*TMath::Pi()*fBMax*fBMax/100;
 228 }
 229
 230 //______________________________________________________________________________
 231 Double_t AliGlauberMC::GetTotXSectErr() const
 232 {
 233    return GetTotXSect()/TMath::Sqrt((Double_t)fEvents) *
 234       TMath::Sqrt(Double_t(1.-fEvents/fTotalEvents));
 235 }
 236
 237 //______________________________________________________________________________
 238 TObjArray *AliGlauberMC::GetNucleons()
 239 {
 240    if(!fNucleonsA || !fNucleonsB) return 0;
 241    fNucleonsA->SetOwner(0);
 242    fNucleonsB->SetOwner(0);
 243    TObjArray *allnucleons=new TObjArray(fAN+fBN);
 244    allnucleons->SetOwner();
 245    for (Int_t i = 0; i<fAN; i++) {
 246       allnucleons->Add(fNucleonsA->At(i));
 247    }
 248    for (Int_t i = 0; i<fBN; i++) {
 249       allnucleons->Add(fNucleonsB->At(i));
 250    }
 251    return allnucleons;
 252 }
 253
 254 //______________________________________________________________________________
 255 Double_t AliGlauberMC::GetdNdEta( const Double_t npp, const Double_t x) const
 256 {
 257   //Get particle density per unit of rapidity
 258   //using the two component model
 259   return (npp*((1.-x)*fNpart/2.+x*fNcoll));
 260 }
 261
 262 //______________________________________________________________________________
 263 Double_t AliGlauberMC::GetdNdEtaGBW( const Double_t delta,
 264                                      const Double_t lambda,
 265                                      const Double_t snn ) const
 266 {
 267   //Get particle density per unit of rapidity
 268   //using the GBW model
 269   return (fNpart*0.47*TMath::Sqrt(TMath::Power(snn,lambda))
 270           * TMath::Power(fNpart,(1.-delta)/3./delta));
 271 }
 272
 273 //______________________________________________________________________________
 274 Double_t AliGlauberMC::GetEccentricityPart() const
 275 {
 276   //get participant eccentricity
 277   return (TMath::Sqrt((fSy2-fSx2)*(fSy2-fSx2)+4*fSxy*fSxy)/(fSy2+fSx2));
 278 }
 279
 280 //______________________________________________________________________________
 281 Double_t AliGlauberMC::GetEccentricity() const
 282 {
 283   //get standard eccentricity
 284   return ((fSy2-fSx2)/(fSy2+fSx2));
 285 }
 286 //______________________________________________________________________________
 287 Bool_t AliGlauberMC::NextEvent(Double_t bgen)
 288 {
 289   //Make a new event
 290   if(bgen<0)
 291      bgen = TMath::Sqrt((fBMax*fBMax-fBMin*fBMin)*gRandom->Rndm()+fBMin*fBMin);
 292
 293   return CalcEvent(bgen);
 294 }
 295
 296 //______________________________________________________________________________
 297 void AliGlauberMC::Run(Int_t nevents)
 298 {
 299    TString name(Form("nt_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
 300    TString title(Form("%s + %s (x-sect = %d mb)",fANucleus.GetName(),fBNucleus.GetName(),(Int_t) fXSect));
 301    if (fnt == 0) {
 302       fnt = new TNtuple(name,title,
 303                         "Npart:Ncoll:B:MeanX:MeanY:MeanX2:MeanY2:MeanXY:VarX:VarY:VarXY:MeanXSystem:MeanYSystem:MeanXA:MeanYA:MeanXB:MeanYB:VarE:VarEpart:Mult:MultGBW");
 304       fnt->SetDirectory(0);
 305    }
 306
 307    for (int i = 0;i<nevents;i++) {
 308
 309       while(!NextEvent()) {}
 310
 311       Float_t v[21];
 312       v[0]  = GetNpart();
 313       v[1]  = GetNcoll();
 314       v[2]  = fB_MC;
 315       v[3]  = fMeanXParts;
 316       v[4]  = fMeanYParts;
 317       v[5]  = fMeanX2;
 318       v[6]  = fMeanY2;
 319       v[7]  = fMeanXY;
 320       v[8]  = fSx2;
 321       v[9]  = fSy2;
 322       v[10] = fSxy;
 323       v[11] = fMeanXSystem;
 324       v[12] = fMeanYSystem;
 325       v[13] = fMeanX_A;
 326       v[14] = fMeanY_A;
 327       v[15] = fMeanX_B;
 328       v[16] = fMeanY_B;
 329       v[17] = GetEccentricity();
 330       v[18] = GetEccentricityPart();
 331       v[19] = GetdNdEta();
 332       v[20] = GetdNdEtaGBW();
 333
 334       fnt->Fill(v);
 335
 336       if (!(i%50)) std::cout << "Event # " << i << " x-sect = " << GetTotXSect() << " +- " << GetTotXSectErr() << " b        \r" << flush;
 337    }
 338    std::cout << std::endl << "Done!" << std::endl;
 339 }
 340
 341 //---------------------------------------------------------------------------------
 342 void AliGlauberMC::runAndSaveNtuple( Int_t n,
 343                                      Option_t *sysA,
 344                                      Option_t *sysB,
 345                                      Double_t signn,
 346                                      Double_t mind,
 347                                      const char *fname)
 348 {
 349    AliGlauberMC *mcg=new AliGlauberMC(sysA,sysB,signn);
 350    mcg->SetMinDistance(mind);
 351    mcg->Run(n);
 352    TNtuple  *nt=mcg->GetNtuple();
 353    TFile out(fname,"recreate",fname,9);
 354    if(nt) nt->Write();
 355    out.Close();
 356 }
 357
 358 //---------------------------------------------------------------------------------
 359 void AliGlauberMC::runAndSaveNucleons( Int_t n,
 360                                        Option_t *sysA,
 361                                        Option_t *sysB,
 362                                        Double_t signn,
 363                                        Double_t mind,
 364                                        Bool_t verbose,
 365                                        const char *fname)
 366 {
 367    AliGlauberMC *mcg=new AliGlauberMC(sysA,sysB,signn);
 368    mcg->SetMinDistance(mind);
 369    TFile *out=0;
 370    if(fname)
 371       out=new TFile(fname,"recreate",fname,9);
 372
 373    for(Int_t ievent=0;ievent<n;ievent++){
 374
 375       //get an event with at least one collision
 376       while(!mcg->NextEvent()) {}
 377
 378       //access, save and (if wanted) print out nucleons
 379       TObjArray* nucleons=mcg->GetNucleons();
 380       if(!nucleons) continue;
 381       if(out)
 382          nucleons->Write(Form("nucleonarray%d",ievent),TObject::kSingleKey);
 383
 384       if(verbose) {
 385          cout<<endl<<endl<<"EVENT NO: "<<ievent<<endl;
 386          cout<<"B = "<<mcg->GetB()<<"  Npart = "<<mcg->GetNpart()<<endl<<endl;
 387          printf("Nucleus\t X\t Y\t Z\tNcoll\n");
 388          Int_t nNucls=nucleons->GetEntries();
 389          for(Int_t iNucl=0;iNucl<nNucls;iNucl++) {
 390             AliGlauberNucleon *nucl=(AliGlauberNucleon *)nucleons->At(iNucl);
 391             Char_t nucleus='A';
 392             if(nucl->IsInNucleusB()) nucleus='B';
 393             Double_t x=nucl->GetX();
 394             Double_t y=nucl->GetY();
 395             Double_t z=nucl->GetZ();
 396             Int_t ncoll=nucl->GetNColl();
 397             printf("   %c\t%2.2f\t%2.2f\t%2.2f\t%3d\n",nucleus,x,y,z,ncoll);
 398          }
 399       }
 400    }
 401    if(out) delete out;
 402 }
 403