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)
  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 Bool_t AliGlauberMC::NextEvent(Double_t bgen)
 256 {
 257    if(bgen<0)
 258       bgen = TMath::Sqrt((fBMax*fBMax-fBMin*fBMin)*gRandom->Rndm()+fBMin*fBMin);
 259
 260    return CalcEvent(bgen);
 261 }
 262
 263 //______________________________________________________________________________
 264 void AliGlauberMC::Run(Int_t nevents)
 265 {
 266    TString name(Form("nt_%s_%s",fANucleus.GetName(),fBNucleus.GetName()));
 267    TString title(Form("%s + %s (x-sect = %d mb)",fANucleus.GetName(),fBNucleus.GetName(),(Int_t) fXSect));
 268    if (fnt == 0) {
 269       fnt = new TNtuple(name,title,
 270                         "Npart:Ncoll:B:MeanX:MeanY:MeanX2:MeanY2:MeanXY:VarX:VarY:VarXY:MeanXSystem:MeanYSystem:MeanXA:MeanYA:MeanXB:MeanYB");
 271       fnt->SetDirectory(0);
 272    }
 273
 274    for (int i = 0;i<nevents;i++) {
 275
 276       while(!NextEvent()) {}
 277
 278       Float_t v[17];
 279       v[0]  = GetNpart();
 280       v[1]  = GetNcoll();
 281       v[2]  = fB_MC;
 282       v[3]  = fMeanXParts;
 283       v[4]  = fMeanYParts;
 284       v[5]  = fMeanX2;
 285       v[6]  = fMeanY2;
 286       v[7]  = fMeanXY;
 287       v[8]  = fSx2;
 288       v[9]  = fSy2;
 289       v[10] = fSxy;
 290       v[11] = fMeanXSystem;
 291       v[12] = fMeanYSystem;
 292       v[13] = fMeanX_A;
 293       v[14] = fMeanY_A;
 294       v[16] = fMeanX_B;
 295       v[17] = fMeanY_B;
 296
 297       fnt->Fill(v);
 298
 299       if (!(i%50)) std::cout << "Event # " << i << " x-sect = " << GetTotXSect() << " +- " << GetTotXSectErr() << " b        \r" << flush;
 300    }
 301    std::cout << std::endl << "Done!" << std::endl;
 302 }
 303
 304 //---------------------------------------------------------------------------------
 305 void AliGlauberMC::runAndSaveNtuple( Int_t n,
 306                                      Option_t *sysA,
 307                                      Option_t *sysB,
 308                                      Double_t signn,
 309                                      Double_t mind,
 310                                      const char *fname)
 311 {
 312    AliGlauberMC *mcg=new AliGlauberMC(sysA,sysB,signn);
 313    mcg->SetMinDistance(mind);
 314    mcg->Run(n);
 315    TNtuple  *nt=mcg->GetNtuple();
 316    TFile out(fname,"recreate",fname,9);
 317    if(nt) nt->Write();
 318    out.Close();
 319 }
 320
 321 //---------------------------------------------------------------------------------
 322 void AliGlauberMC::runAndSaveNucleons( Int_t n,
 323                                        Option_t *sysA,
 324                                        Option_t *sysB,
 325                                        Double_t signn,
 326                                        Double_t mind,
 327                                        Bool_t verbose,
 328                                        const char *fname)
 329 {
 330    AliGlauberMC *mcg=new AliGlauberMC(sysA,sysB,signn);
 331    mcg->SetMinDistance(mind);
 332    TFile *out=0;
 333    if(fname)
 334       out=new TFile(fname,"recreate",fname,9);
 335
 336    for(Int_t ievent=0;ievent<n;ievent++){
 337
 338       //get an event with at least one collision
 339       while(!mcg->NextEvent()) {}
 340
 341       //access, save and (if wanted) print out nucleons
 342       TObjArray* nucleons=mcg->GetNucleons();
 343       if(!nucleons) continue;
 344       if(out)
 345          nucleons->Write(Form("nucleonarray%d",ievent),TObject::kSingleKey);
 346
 347       if(verbose) {
 348          cout<<endl<<endl<<"EVENT NO: "<<ievent<<endl;
 349          cout<<"B = "<<mcg->GetB()<<"  Npart = "<<mcg->GetNpart()<<endl<<endl;
 350          printf("Nucleus\t X\t Y\t Z\tNcoll\n");
 351          Int_t nNucls=nucleons->GetEntries();
 352          for(Int_t iNucl=0;iNucl<nNucls;iNucl++) {
 353             AliGlauberNucleon *nucl=(AliGlauberNucleon *)nucleons->At(iNucl);
 354             Char_t nucleus='A';
 355             if(nucl->IsInNucleusB()) nucleus='B';
 356             Double_t x=nucl->GetX();
 357             Double_t y=nucl->GetY();
 358             Double_t z=nucl->GetZ();
 359             Int_t ncoll=nucl->GetNColl();
 360             printf("   %c\t%2.2f\t%2.2f\t%2.2f\t%3d\n",nucleus,x,y,z,ncoll);
 361          }
 362       }
 363    }
 364    if(out) delete out;
 365 }
 366