PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.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 #include <TNamed.h>
  17 #include <TH1D.h>
  18 #include <TString.h>
  19 #include <TFile.h>
  20 #include <TMath.h>
  21 #include <TRandom1.h>
  22 #include <TROOT.h>
  23 #include <TH2F.h>
  24 #include <TF1.h>
  25 #include <TNtuple.h>
  26 #include <TStyle.h>
  27 #include <TGraphErrors.h>
  28 #include <vector>
  29 #include <TMinuit.h>
  30 #include <TCanvas.h>
  31 #include <TRandom.h>
  32 #include <TDatabasePDG.h>
  33 #include <TPDGCode.h>
  34 #include "AliCentralityGlauberFit.h"
  35 #include "AliLog.h"
  36
  37 ClassImp(AliCentralityGlauberFit)
  38
  39 //______________________________________________________________________________
  40 AliCentralityGlauberFit::AliCentralityGlauberFit(const char *filename) :
  41   fNk(0),
  42   fKlow(0),
  43   fKhigh(0),
  44   fNalpha(0),
  45   fAlphalow(0),
  46   fAlphahigh(0),
  47   fNsigma(0),
  48   fSigmalow(0),
  49   fSigmahigh(0),
  50   fNbog(0),
  51   fBoglow(0),
  52   fBoghigh(0),
  53   fNgamma(0),
  54   fgammalow(0),
  55   fgammahigh(0),
  56   fNsigmap(0),
  57   fsigmaplow(0),
  58   fsigmaphigh(0),
  59   fRebinFactor(0),
  60   fScalemin(0),
  61   fMultmin(0),
  62   fMultmax(0),
  63   fGlauntuple(0),
  64   fNpart(0),
  65   fNcoll(0),
  66   fB(0),
  67   fTaa(0),
  68   fTempHist(0),
  69   fGlauberHist(0),
  70   fUseChi2(kTRUE),
  71   fNevents(100000),
  72   fIsZN(kTRUE),
  73   fNTrials(0),
  74   fChi2Min(2000.),
  75   fInrootfilename(0),
  76   fInntuplename(0),
  77   fOutrootfilename(0),
  78   fOutntuplename(0),
  79   fAncfilename("ancestor_hists.root"),
  80   fHistnames()
  81 {
  82   // Standard constructor.
  83   TFile *f = 0;
  84   if (filename) {  // glauber file
  85     f = TFile::Open(filename);
  86     fGlauntuple = (TNtuple*) f->Get("nt_p_Pb");
  87     fGlauntuple->SetBranchAddress("Npart",&fNpart);
  88     fGlauntuple->SetBranchAddress("Ncoll",&fNcoll);
  89     fGlauntuple->SetBranchAddress("B",&fB);
  90     fGlauntuple->SetBranchAddress("tAA",&fTaa);
  91   }
  92   for(int i=0; i<4; i++) fParamnSlow[i] = 0.;
  93
  94 }
  95
  96 //--------------------------------------------------------------------------------------------------
  97 void AliCentralityGlauberFit::SetRangeToFit(Double_t fmultmin, Double_t fmultmax)
  98 {
  99   // Set fit range.
 100
 101   fMultmin=fmultmin;
 102   fMultmax=fmultmax;
 103 }
 104
 105 //--------------------------------------------------------------------------------------------------
 106 void AliCentralityGlauberFit::SetRangeToScale(Double_t fscalemin)
 107 {
 108   // Set range where to scale simulated histo to data
 109
 110   fScalemin=fscalemin;
 111 }
 112
 113 //--------------------------------------------------------------------------------------------------
 114 void AliCentralityGlauberFit::SetGlauberParam(
 115                                               Int_t   Nk,
 116                                               Double_t klow,
 117                                               Double_t khigh,
 118                                               Int_t   Nalpha,
 119                                               Double_t alphalow,
 120                                               Double_t alphahigh,
 121                                               Int_t   Nsigma,
 122                                               Double_t sigmalow,
 123                                               Double_t sigmahigh,
 124                                               Int_t   Nbog,
 125                                               Double_t boglow,
 126                                               Double_t boghigh,
 127                                               Int_t   Ngamma,
 128                                               Double_t gammalow,
 129                                               Double_t gammahigh,
 130                                               Int_t   Nsigmap,
 131                                               Double_t sigmaplow,
 132                                               Double_t sigmaphigh)
 133 {
 134   // Set Glauber parameters.
 135   fNk=Nk;
 136   fKlow=klow;
 137   fKhigh=khigh;
 138   fNalpha=Nalpha;
 139   fAlphalow=alphalow;
 140   fAlphahigh=alphahigh;
 141   fNsigma=Nsigma;
 142   fSigmalow=sigmalow;
 143   fSigmahigh=sigmahigh;
 144   fNbog=Nbog;
 145   fBoglow=boglow;
 146   fBoghigh=boghigh;
 147   fNgamma=Ngamma;
 148   fgammalow=gammalow;
 149   fgammahigh=gammahigh;
 150   fNsigmap=Nsigmap;
 151   fsigmaplow=sigmaplow;
 152   fsigmaphigh=sigmaphigh;
 153 }
 154
 155 //--------------------------------------------------------------------------------------------------
 156 void AliCentralityGlauberFit::SetNParam(Double_t a, Double_t b, Double_t c, Double_t d)
 157 {
 158    // Setting parameters that are adjusted
 159    fParamnSlow[0] = a;
 160    fParamnSlow[1] = b;
 161    fParamnSlow[2] = c;
 162    fParamnSlow[3] = d;
 163    for(int i=0; i<4; i++) printf("  INITIAL PARAMETER VALUES : paramnSlow[%d] = %f\n", i, fParamnSlow[i]);
 164 }
 165
 166 //--------------------------------------------------------------------------------------------------
 167 void AliCentralityGlauberFit::MakeFits()
 168 {
 169   // Make fits.
 170
 171   TH1F *hDATA = 0x0;
 172   TH1F *thistGlau;
 173   TFile *inrootfile;
 174   TFile *outrootfile;
 175
 176   // open inrootfile, outrootfile
 177   std::cout << "input file "  << fInrootfilename  << std::endl;
 178   std::cout << "output file " << fOutrootfilename << std::endl << std::endl;
 179   inrootfile  = new TFile(fInrootfilename,"OPEN");
 180   outrootfile = new TFile(fOutrootfilename,"RECREATE");
 181
 182   // loop over all distribution names
 183   std::vector<TString>::const_iterator hni;
 184   for(hni=fHistnames.begin(); hni!=fHistnames.end(); hni++) {
 185     hDATA  = (TH1F*) (inrootfile->Get(*hni));
 186     std::cout << " ->  Getting histogram " << *hni << std::endl << std::endl;
 187     if (!hDATA) {
 188       printf("   @@@@@@ No DATA histo in input -> EXITING!!!!!\n\n");
 189       return;
 190       //TList *list  = (TList*) (inrootfile->Get("CentralityStat"));
 191       //hDATA  = (TH1F*) (list->FindObject(*hni));
 192     }
 193     //hDATA->Rebin(fRebinFactor);
 194     //TH1F *hGLAU = new TH1F("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
 195
 196     Double_t gamma_min=-1;
 197     Double_t sigmap_min=-1;
 198     Double_t bog_min=-1;
 199     Double_t sigma_min=-1;
 200     Double_t alpha_min=-1;
 201     Double_t k_min=-1, chi2=0.;
 202     Double_t alpha=0., k=0., sigma=0., bog=0., gamma=0., sigmap=0.;
 203
 204     for(Int_t nsigmap=0; nsigmap<fNsigmap; nsigmap++) {
 205       sigmap = fsigmaplow   + ((Double_t) nsigmap ) * (fsigmaphigh  - fsigmaplow ) / fNsigmap;
 206
 207      for(Int_t ngamma=0; ngamma<fNgamma; ngamma++) {
 208       gamma = fgammalow   + ((Double_t) ngamma ) * (fgammahigh  - fgammalow ) / fNgamma;
 209
 210       for(Int_t nbog=0; nbog<fNbog; nbog++) {
 211         bog = fBoglow   + ((Double_t) nbog ) * (fBoghigh  - fBoglow ) / fNbog;
 212
 213         for(Int_t nsigma=0; nsigma<fNsigma; nsigma++) {
 214           sigma = fSigmalow   + ((Double_t) nsigma ) * (fSigmahigh  - fSigmalow ) / fNsigma;
 215
 216           for(Int_t nalpha=0; nalpha<fNalpha; nalpha++) {
 217             alpha = fAlphalow   + ((Double_t) nalpha ) * (fAlphahigh  - fAlphalow ) / fNalpha;
 218
 219             for(Int_t nk=0; nk<fNk; nk++) {
 220               k = fKlow  + ((Double_t) nk ) * (fKhigh  - fKlow ) / fNk;
 221
 222
 223                for(int i=0; i<fNTrials; i++){
 224                  printf("  **** calling GlauberHisto with acc = %1.3f R = %1.4f k = %1.4f bog = %1.3f  gamma = %1.3f sigmap = %1.3f\n",
 225                         alpha, sigma, k, bog, gamma, sigmap);
 226                  if(fIsZN){
 227                    // Updating parameters for slow n
 228                    fParamnSlow[0] = fParamnSlow[0]; // + 0.5*i;
 229                    fParamnSlow[1] = fParamnSlow[1]; // + 10.*i;
 230                    fParamnSlow[2] = fParamnSlow[2] + 0.05*i;
 231                    fParamnSlow[3] = fParamnSlow[3];// + 0.01*i;
 232                    //
 233                    printf("  \t\t trial #%d  n slow param.: %1.2f %1.2f %1.2f %1.2f\n",
 234                         i,fParamnSlow[0], fParamnSlow[1], fParamnSlow[2], fParamnSlow[3]);
 235                  }
 236
 237                  //thistGlau = GlauberHisto(k,alpha,sigma,bog,gamma,sigmap,hDATA,kFALSE);
 238                  thistGlau = GlauberHisto(k,alpha,sigma,bog,gamma,sigmap,hDATA,kTRUE);
 239                  //
 240                  chi2 = CalculateChi2(hDATA, thistGlau);
 241                  printf("   ->  Chi2 %f \n", chi2);
 242                  if(chi2 < fChi2Min){
 243                    fChi2Min = chi2;
 244                    alpha_min = alpha;
 245                    sigma_min = sigma;
 246                    bog_min = bog;
 247                    gamma_min = gamma;
 248                    sigmap_min = sigmap;
 249                    k_min = k;
 250                  }
 251                }
 252             }
 253           }
 254         }
 255       }
 256      }
 257     }
 258     if(fNsigmap>1 || fNgamma>1 || fNbog>1 || fNsigma>1 || fNalpha>1 || fNk>1 || fNTrials>1){
 259       thistGlau->Reset();
 260       printf("  **** calling GlauberHisto with acc = %1.3f R = %1.4f bog = %1.3f  gamma = %1.3f sigmap = %1.3f\n",
 261                 alpha, sigma, bog, gamma, sigmap);
 262       thistGlau = GlauberHisto(k_min,alpha_min,sigma_min,bog_min,gamma_min, sigmap_min, hDATA,kTRUE);
 263     }
 264
 265     TH1F * hGLAU = 0x0;
 266     hGLAU = (TH1F *) thistGlau->Clone("hGLAU");
 267
 268     hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU")));
 269     if(fIsZN) hGLAU->SetTitle(((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f_%.2f-%.2f-%.3f-%.2f",
 270                                       bog,sigmap,gamma,fParamnSlow[0], fParamnSlow[1], fParamnSlow[2], fParamnSlow[3])));
 271     else hGLAU->SetTitle(((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f",
 272                                       alpha,bog,sigmap)));
 273
 274     Int_t firstBin=1;
 275     if(fIsZN) firstBin = 8;//ZN
 276     Int_t lastBin=0;
 277     if(fIsZN) lastBin = 800;//ZN
 278     else lastBin = hDATA->GetNbinsX(); // ZP
 279     //
 280     Double_t mcintegral = hGLAU->Integral(firstBin, lastBin,"width");
 281     Double_t dataintegral = hDATA->Integral(firstBin, lastBin,"width");
 282     Double_t scale = 1.;
 283     if(mcintegral>0.) (scale = dataintegral/mcintegral);
 284     printf("\n scale %f \n", scale);
 285     //
 286     std::cout << "hGLAU -> Integral in bin range:" << firstBin << "-" << lastBin << " = " <<
 287         mcintegral << std::endl;
 288     std::cout << "hDATA -> Integral in bin range:" << firstBin << "-" << lastBin << " = " <<
 289         dataintegral << std::endl;
 290
 291     hGLAU->Scale(scale);
 292     //
 293     std::cout << "hGLAU -> Integral (whole range): " << hGLAU->Integral(1, hGLAU->GetNbinsX(),"width") << std::endl;
 294     std::cout << "hDATA -> Integral (whole range): " << hDATA->Integral(1, hDATA->GetNbinsX(),"width") << std::endl;
 295
 296     TCanvas *c2 = new TCanvas("c2","checkGLAU",100,0,700,700);
 297     c2->cd();
 298     hGLAU->SetLineColor(kPink);
 299     hGLAU->Draw("");
 300     hDATA->Draw("same");
 301
 302     SaveHisto(hDATA, hGLAU, outrootfile);
 303
 304     printf("\n \t k = %1.4f  alpha = %1.3f  sigma = %1.2f  bog = %1.4f  gamma = %1.4f sigmap = %.3f \n", k,alpha, sigma, bog, gamma, sigmap);
 305     if(fIsZN) printf(" \t a = %1.2f  b = %1.3f  c = %1.2f  d = %1.2f\n\n", fParamnSlow[0], fParamnSlow[1], fParamnSlow[2], fParamnSlow[3]);
 306     std::cout << "chi2 min is " << fChi2Min << std::endl << std::endl;
 307     std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin),
 308                                               hGLAU->FindBin(fMultmax))/hGLAU->Integral()
 309               << " of the total cross section" << std::endl << std::endl;
 310     fTempHist = hDATA;
 311     fGlauberHist = hGLAU;
 312   }
 313
 314   // close inrootfile, outrootfile
 315   inrootfile->Close();
 316   outrootfile->Close();
 317 }
 318
 319
 320 //--------------------------------------------------------------------------------------------------
 321 TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t k, Double_t alpha, Double_t sigma, Double_t bog,
 322                 Double_t gamma, Double_t sigmap, TH1F *hDATA, Bool_t save)
 323 {
 324   // Get Glauber histogram.
 325   static TH1F *h1 = (TH1F*) hDATA->Clone("h1");
 326   h1->Reset();
 327   h1->SetName(Form("fit_%.3f_%.3f",k,alpha));
 328
 329   TFile *outFile = NULL;
 330   TNtuple *ntuple = NULL;
 331
 332   Double_t ntot=0, nblack=0, ngray=0, Etot=0, lcp=0, nslowp=0, nslown=0;
 333
 334   if(save){
 335     outFile = new TFile(fOutntuplename,"RECREATE");
 336     ntuple = new TNtuple("gnt", "Glauber ntuple", "fNpart:fNcoll:fB:fTaa:ntot:nblack:ngray:Etot:lcp:nslowp:nslown");
 337   }
 338
 339   Int_t nents=0, evtot=0, evn=0, evzeron=0, evzerop=0, evenzero=0;
 340   if(fGlauntuple) nents = fGlauntuple->GetEntries();
 341
 342   for (Int_t i=0;i<fNevents;++i) {
 343     if(fGlauntuple) fGlauntuple->GetEntry(i % nents);
 344     else{
 345       printf(" NOT GETTING ENTRY FROM TNTUPLE for ev. %d -> setting Ncoll=1 \n\n",i);
 346       fNpart = 2;
 347       fNcoll = 1;
 348     }
 349     //printf(" Getting entry %d from ntuple  Ncoll = %f\n",i%nents, fNcoll);
 350
 351     // Slow Nucleon Model from Chiara
 352     Double_t n=0., nbn=0., ngn=0., nbp=0., ngp=0.;
 353     //
 354     //MakeSlowNucleons(fNcoll,nbn,ngn,nbp,ngp); // pure Sikler
 355     //MakeSlowNucleons2(fNcoll,bog,gamma, nbn,ngn,nbp,ngp,lcp); // smear in Ncoll
 356     MakeSlowNucleons2s(fNcoll,bog,gamma,sigmap, nbn,ngn,nbp,ngp,lcp);  // smear in Nslow
 357
 358     evtot++;
 359     if(fNcoll==1) evn++;
 360     //
 361     if(fIsZN){
 362       n = alpha*(ngn+nbn);
 363       nblack = nbn;
 364       ngray = ngn;
 365     }
 366     else{
 367       n = alpha*(ngp+nbp); // ZPA
 368       nblack = nbp;
 369       ngray = ngp;
 370     }
 371     nslown = nbn+ngn;
 372     nslowp = nbp+ngp;
 373
 374     if((alpha*(ngn+nbn)<0.5)) evzeron++;
 375     if((alpha*(ngp+nbp)<0.5)) evzerop++;
 376
 377     //------ experimental resolution -> Gaussian smearing ------------------------------------
 378     Double_t resexp = 0.;
 379     //if (n>0) resexp = sigma*gRandom->Gaus(0.0,1.0)/TMath::Sqrt(n);
 380     //else resexp=0;
 381     //ntot = n*(1+resexp);
 382
 383     if(n>=0){
 384       resexp = sigma*TMath::Sqrt(n);
 385       ntot = (gRandom->Gaus(n, resexp));
 386     }
 387     //else n=0.;
 388
 389     // non-lineary effect -------------------------------
 390     //ntot = ntot + k*ntot*ntot;
 391
 392     Double_t nFact = 4.*82./208.;
 393     Etot = nFact*n;
 394
 395     // USE ONLY IF ZDC TDC IS REQUESTED TO CONDITION THE SPECTRUM!!!!!!!!!
 396     float resexpe=0.;
 397     if(n>0.){
 398       resexpe = 1./TMath::Sqrt(n)*sigma*gRandom->Gaus(0.0, 1.0);
 399       //resexpe = 1./TMath::Sqrt(n)*sigma;
 400       //
 401       Etot = Etot*(1+resexpe);
 402     }
 403     // non-lineary effect -------------------------------
 404     //Etot = Etot + k*Etot*Etot;
 405
 406     if(Etot<0.5) evenzero++;
 407
 408     h1->Fill(Etot);
 409     if(save) ntuple->Fill(fNpart, fNcoll, fB, fTaa, ntot, nblack, ngray, Etot, lcp, nslowp, nslown);
 410   }
 411
 412   if(save) {
 413     //printf("\n ******** Fraction of events with Ncoll=1 in the ntuple %f\n",(float) (1.*evn/evtot));
 414     if(fIsZN){
 415       printf(" ******** Fraction of events with no neutrons %f\n",(float) (1.*evzeron/evtot));
 416       printf(" ******** Fraction of events with no neutron energy %f\n",(float) (1.*evenzero/evtot));
 417       printf(" DATA: fraction of events with no ZN signal = %f\n\n", 1-alpha);
 418     }
 419     else printf("\n ******** Fraction of events with no protons %f \n",(float) (1.*evzerop/evtot));
 420
 421     ntuple->Write();
 422     outFile->Close();
 423   }
 424
 425   return h1;
 426 }
 427
 428 //--------------------------------------------------------------------------------------------------
 429 Double_t AliCentralityGlauberFit::CalculateChi2(TH1F *hDATA, TH1F *thistGlau)
 430 {
 431   // Note that for different values of neff the mc distribution is identical, just re-scaled below
 432   // normalize the two histogram, scale MC up but leave data alone - that preserves error bars !!!!
 433
 434   Int_t lowchibin  =  8;
 435   Int_t highchibin =  0;
 436   if(fIsZN) highchibin = 800;
 437   else highchibin = 25;
 438
 439   Double_t mcintegral = thistGlau->Integral(1, thistGlau->GetNbinsX());
 440   Double_t scale = (hDATA->Integral(1, hDATA->GetNbinsX())/mcintegral);
 441   thistGlau->Scale(scale);
 442
 443   // calculate the chi2 between MC and real data over some range ????
 444   if (fUseChi2) {
 445     Double_t chi2 = 0.0;
 446     for (Int_t i=lowchibin; i<=highchibin; i++) {
 447       if (hDATA->GetBinContent(i) < 1.0) continue;
 448       Double_t diff = TMath::Power((thistGlau->GetBinContent(i) - hDATA->GetBinContent(i)),2);
 449       diff = diff / (TMath::Power(hDATA->GetBinError(i),2) + TMath::Power(thistGlau->GetBinError(i),2));
 450       chi2 += diff;
 451     }
 452     chi2 = chi2 / (highchibin - lowchibin + 1);
 453     return chi2;
 454   }
 455   // "-2 log likelihood ratio(mu;n) = 2[(mu - n) + n ln(n/mu)]"
 456   else {
 457     std::cout << "LL" << std::endl;
 458
 459     Double_t ll = 0.0;
 460     for (Int_t i=lowchibin; i<=highchibin; i++) {
 461       Double_t data = hDATA    ->GetBinContent(i);
 462       Double_t mc   = thistGlau->GetBinContent(i);
 463       Int_t idata = TMath::Nint(data);
 464       if (mc < 1.e-9) mc = 1.e-9;
 465       Double_t fsub = - mc + idata * TMath::Log(mc);
 466       Double_t fobs = 0;
 467       if (idata > 0) {
 468         for(Int_t istep = 0; istep < idata; istep++){
 469           if (istep > 1)
 470             fobs += TMath::Log(istep);
 471         }
 472       }
 473       fsub -= fobs;
 474       ll -=  fsub ;
 475     }
 476     return 2*ll;
 477   }
 478 }
 479
 480 //--------------------------------------------------------------------------------------------------
 481 void AliCentralityGlauberFit::SaveHisto(TH1F *hist1, TH1F *hist2, TFile *outrootfile)
 482 {
 483   // Save histograms.
 484
 485   outrootfile->cd();
 486   hist1->Write();
 487   hist2->Write();
 488 }
 489
 490 //--------------------------------------------------------------------------------------------------
 491 void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t &nbn, Double_t &ngn,Double_t &nbp, Double_t &ngp)
 492 {
 493 // from AliGenSlowNucleonModelExp (Chiara Oppedisano)
 494 // Return the number of black and gray nucleons
 495 //
 496 // Number of collisions
 497
 498   Float_t fP = 82;
 499   Float_t fN = 126;
 500   Float_t nu = (Float_t) ncoll;
 501   //
 502 // Mean number of gray nucleons
 503     Float_t nGray         = 2. * nu;
 504     Float_t nGrayNeutrons = nGray * fN / (fN + fP);
 505     Float_t nGrayProtons  = nGray - nGrayNeutrons;
 506
 507 // Mean number of black nucleons
 508     Float_t nBlack  = 0.;
 509     if(nGray>=15.) nBlack = 28.;
 510     Float_t nBlackNeutrons = nBlack * 0.84;
 511     Float_t nBlackProtons  = nBlack - nBlackNeutrons;
 512
 513 // Actual number (including fluctuations) from binomial distribution
 514 //  gray neutrons
 515     ngn = (double) gRandom->Binomial((Int_t) fN, nGrayNeutrons/fN);
 516
 517 //  gray protons
 518     ngp = (double) gRandom->Binomial((Int_t) fP, nGrayProtons/fP);
 519
 520 //  black neutrons
 521     nbn = (double) gRandom->Binomial((Int_t) fN, nBlackNeutrons/fN);
 522
 523 //  black protons
 524     nbp = (double) gRandom->Binomial((Int_t) fP, nBlackProtons/fP);
 525
 526    // printf(" Ncoll %d  ngp %f nbp %f  ngn %f nbn %f \n",ncoll,ngp, nbp, ngn, nbn);
 527
 528 }
 529
 530
 531 //--------------------------------------------------------------------------------------------------
 532 void AliCentralityGlauberFit::MakeSlowNucleons2(Int_t ncoll, Double_t bog, Double_t gamma,
 533         Double_t &nbn, Double_t &ngn, Double_t &nbp, Double_t &ngp, Double_t &lcp)
 534 {
 535 // from AliGenSlowNucleonModelExp (Chiara Oppedisano)
 536 //
 537 // Return the number of black and gray nucleons
 538 //
 539
 540  // PROTONS ----------------------------------------------------------------
 541
 542   Int_t fP = 82;
 543   Int_t fN = 126;
 544
 545   Float_t nu = (Float_t) ncoll;
 546   //
 547   nu = gRandom->Gaus(nu,0.5);
 548   if(nu<1.) nu = 1.;
 549
 550   // PROTONS ----------------------------------------------------------------
 551   //  gray protons
 552   Float_t  poverpd = 0.843;
 553   Float_t  zAu2zPb = 82./79.;
 554   Float_t  nGrayp = (-0.27 + 0.63 * nu - 0.0008 *nu *nu)*poverpd*zAu2zPb;
 555   if(nGrayp<0.) nGrayp=0;
 556
 557   Double_t p;
 558   p =  nGrayp/fP;
 559   ngp = (double) gRandom->Binomial(fP, p);
 560
 561   //  black protons
 562   //Float_t blackovergray = 3./7.;// from spallation
 563   //Float_t blackovergray = 0.65; // from COSY
 564   Float_t blackovergray = bog;
 565   Float_t nBlackp  = blackovergray*nGrayp;
 566   if(nBlackp<0.) nBlackp=0;
 567
 568   p =  nBlackp/(fP-ngp);
 569   nbp = (double) gRandom->Binomial((Int_t) (fP-ngp), p);
 570   // SATURATION in no. of black protons
 571   //if(ngp>=9.) nbp = 15;
 572
 573  // NEUTRONS ----------------------------------------------------------------
 574
 575   /*if(nu<3.){
 576     nGrayp = -0.836 + 0.9112 *nu - 0.05381 *nu *nu;
 577     nBlackp  = blackovergray*nGrayp;
 578   } */
 579
 580   //  gray neutrons
 581   Float_t nGrayNeutrons = 0.;
 582   Float_t nBlackNeutrons = 0.;
 583   lcp = (nGrayp+nBlackp)/gamma;
 584
 585   Float_t nSlow  = 0.;
 586   if(lcp>0.){
 587     nSlow  = fParamnSlow[0]-fParamnSlow[1]/(fParamnSlow[2]+lcp)+fParamnSlow[3]*lcp;
 588     //
 589 //changed
 590 //    float xconj = fParamnSlow[1]/fParamnSlow[0]-fParamnSlow[2];
 591     float xconj = 1.;
 592     float yconj = fParamnSlow[0]-fParamnSlow[1]/(fParamnSlow[2]+xconj)+fParamnSlow[3]*xconj;
 593     if(lcp<xconj) nSlow = 0.+(yconj-0.)*(lcp-0.)/(xconj-0.);
 594     //
 595     // Factor to scale COSY to obtain ALICE n mltiplicities
 596     //    nSlow = 1.68*nSlow;  // now integrated in fParamnSlow[0], fParamnSlow[1]
 597     //
 598     nBlackNeutrons = nSlow * 0.9;
 599     //nBlackNeutrons = nSlow * 0.5;
 600     nGrayNeutrons  = nSlow - nBlackNeutrons;
 601   }
 602   else{
 603     // slightly adjusted Sikler corrected for neutron yield...
 604     nBlackNeutrons = 126./208. * 4. * nu;
 605     //nGrayNeutrons  = 126./208. * 2. * nu; // a la Sikler (<Ngn>~0.50*<Nbn>)
 606     //nGrayNeutrons  = 0.47 * 2. * nu; // a la DPMJET (<Ngn>~0.39*<Nbn>)
 607     nGrayNeutrons = nBlackNeutrons/9.;   // a la spallation (<Ngn>~0.11*<Nbn>)
 608   }
 609   //
 610   if(nGrayNeutrons<0.) nGrayNeutrons=0.;
 611   if(nBlackNeutrons<0.) nBlackNeutrons=0.;
 612
 613   p =  nGrayNeutrons/fN;
 614   ngn = gRandom->Binomial(fN, p);
 615 //changed
 616   if(nGrayNeutrons>=10) ngn = gRandom->Gaus(nGrayNeutrons+0.5, TMath::Sqrt(fN*p*(1-p)));
 617   //else ngn = gRandom->PoissonD(nGrayNeutrons);
 618   //else ngn = gRandom->Binomial((Int_t) fN, p);
 619   if(ngn<0.) ngn=0.;
 620
 621   //  black neutrons
 622   p =  nBlackNeutrons/(fN-ngn);
 623   nbn = gRandom->Binomial((Int_t) (fN-ngn), p);
 624 //changed
 625   if(nBlackNeutrons>=10) nbn = gRandom->Gaus(nBlackNeutrons+0.5, TMath::Sqrt(fN*p*(1-p)));
 626   else nbn = gRandom->PoissonD(nBlackNeutrons);
 627   //else nbn = gRandom->Binomial((Int_t) fN, p);
 628   if(nbn<0.) nbn=0.;
 629
 630 }
 631
 632 //--------------------------------------------------------------------------------------------------
 633  void AliCentralityGlauberFit::MakeSlowNucleons2s(Float_t ncoll, Double_t bog, Double_t gamma, Double_t sigmap,
 634         Double_t &nbn, Double_t &ngn, Double_t &nbp, Double_t &ngp, Double_t &lcp)
 635 {
 636 // from AliGenSlowNucleonModelExp (Chiara Oppedisano)
 637 //
 638 // Return the number of black and gray nucleons
 639 //
 640
 641   float fP = 82;
 642   float fN = 126;
 643
 644   Float_t nu = ncoll;
 645   //Float_t nu = ncoll*(gRandom->Rndm()+0.5);
 646
 647  // PROTONS ----------------------------------------------------------------
 648   //  gray protons
 649   Float_t  poverpd = 0.843;
 650   Float_t  zAu2zPb = 82./79.;
 651   Float_t  grayp = (-0.27 + 0.63 * nu - 0.0008 *nu *nu)*poverpd*zAu2zPb;
 652   //if(nu<3.) grayp = -0.836 + 0.9112 *nu - 0.05381 *nu *nu;
 653   //
 654 //changed
 655   //Float_t  nGrayp = grayp;
 656   Float_t  nGrayp = gRandom->Gaus(grayp, sigmap);
 657 //changed
 658 //  Float_t  nGrayp = gRandom->Gaus(grayp, (0.8*sigmap-sigmap)*grayp/90+sigmap);
 659   if(nGrayp<0.) nGrayp=0;
 660
 661   Double_t p=0.;
 662   p =  nGrayp/fP;
 663   ngp = (double) gRandom->Binomial((int) fP, p);
 664
 665   //  black protons
 666   //Float_t blackovergray = 3./7.;// =0.43 from spallation
 667   //Float_t blackovergray = 0.65; // from COSY
 668   Float_t blackovergray = bog;
 669   //
 670   Float_t nBlackp = blackovergray*nGrayp;
 671   //if(nBlackp<0.) nBlackp=0;
 672
 673   p =  nBlackp/(fP-ngp);
 674   nbp = (double) gRandom->Binomial((int) (fP-ngp), p);
 675   // SATURATION in no. of black protons
 676   //if(ngp>=9.) nbp = 15;
 677
 678  // NEUTRONS ----------------------------------------------------------------
 679   // Relevant ONLY for n
 680   // NB-> DOESN'T WORK FOR SMEARING IN Nslow!!!!!!!!!!
 681   /*if(nu<=3.){
 682     grayp = -0.836 + 0.9112 *nu - 0.05381 *nu *nu;
 683     // smearing
 684     //nGrayp = gRandom->Gaus(grayp, sigmap);
 685     nGrayp = grayp;
 686     nBlackp  = blackovergray*nGrayp;
 687   }*/
 688
 689   //  gray neutrons
 690   lcp = (nGrayp+nBlackp)/gamma;
 691
 692   Float_t nGrayNeutrons = 0.;
 693   Float_t nBlackNeutrons = 0.;
 694   Float_t nSlow  = 0.;
 695   if(lcp>0.){
 696     nSlow  = fParamnSlow[0]-fParamnSlow[1]/(fParamnSlow[2]+lcp)+fParamnSlow[3]*lcp;
 697     //
 698 //changed
 699 //    float xconj = fParamnSlow[1]/fParamnSlow[0]-fParamnSlow[2];
 700     /*float xconj = 1.;
 701     float yconj = fParamnSlow[0]-fParamnSlow[1]/(fParamnSlow[2]+xconj)+fParamnSlow[3]*xconj;
 702     if(lcp<xconj) nSlow = 0.+(yconj-0.)*(lcp-0.)/(xconj-0.);*/
 703     //
 704     // Factor to scale COSY to obtain ALICE n mltiplicities
 705     //    nSlow = 1.68*nSlow;  // now integrated in fParamnSlow[0], fParamnSlow[1]
 706     //
 707     nBlackNeutrons = nSlow * 0.9;
 708     //nBlackNeutrons = nSlow * 0.5;
 709     nGrayNeutrons  = nSlow - nBlackNeutrons;
 710   }
 711   else{
 712    // taking into account the prob4bility p to have 0 neutrons when we have 0 protons
 713    // p = 1.-(0.82*0.96) ~ 0.22 (assuming that 100% of ev. without n are also without p)
 714    //Float_t r = gRandom->Rndm();
 715    //if(r>0.1){
 716      // Sikler corrected for neutron yield...
 717      nBlackNeutrons = 126./208. * 4. * nu;
 718      //nGrayNeutrons  = 126./208. * 2. * nu; // a la Sikler (<Ngn>~0.50*<Nbn>)
 719      //nGrayNeutrons  = 0.47 * 2. * nu; // a la DPMJET (<Ngn>~0.39*<Nbn>)
 720      //nGrayNeutrons = nBlackNeutrons/9.;   // a la spallation (<Ngn>~0.11*<Nbn>)
 721      // smearing!
 722 //changed
 723      nBlackNeutrons = gRandom->Gaus(nBlackNeutrons, sigmap);
 724      //nGrayNeutrons  = nBlackNeutrons/2.; // a la Sikler (<Ngn>~0.50*<Nbn>)
 725      nGrayNeutrons = nBlackNeutrons/9.;   // a la spallation (<Ngn>~0.11*<Nbn>)
 726      //printf("  Sikler -> Ncoll %f <Ngrayn> %f  <Nblackn> %f \n",nu,nGrayNeutrons, nBlackNeutrons);
 727    //}
 728   }
 729   //
 730   if(nGrayNeutrons<0.) nGrayNeutrons=0.;
 731   if(nBlackNeutrons<0.) nBlackNeutrons=0.;
 732
 733   //  black neutrons
 734   p =  nBlackNeutrons/fN;
 735   nbn = gRandom->Binomial((Int_t) fN, p);
 736   TRandom1 r;
 737   //nbn = (double) r.Binomial((int) fN, p);
 738 //changed
 739 //  if(nBlackNeutrons>=10) nbn = r.Gaus(nBlackNeutrons+0.5, TMath::Sqrt(fN*p*(1-p)));
 740   if(nbn<0.) nbn=0.;
 741
 742   //  gray neutrons
 743   p =  nGrayNeutrons/(fN-nbn);
 744   ngn = gRandom->Binomial((Int_t) (fN-nbn), p);
 745   //ngn = (double) (r.Binomial((Int_t) (fN-nbn), p));
 746 //changed
 747   /*if(nGrayNeutrons>=10) ngn = gRandom->Gaus(nGrayNeutrons+0.5, TMath::Sqrt(fN*p*(1-p)));
 748   else ngn = gRandom->Binomial((Int_t) fN, p);*/
 749   if(ngn<0.) ngn=0.;
 750
 751   //printf(" Ncoll %1.2f   Ngp %1.2f Nbp %1.2f   Nbn %1.2f Ngn %1.2f\n",nu, ngp,nbp,nbn,ngn);
 752
 753 }
 754
 755
 756 Double_t AliCentralityGlauberFit::ConvertToEnergy(Double_t T)
 757 {
 758    TDatabasePDG * pdg = TDatabasePDG::Instance();
 759     const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
 760     Double_t fPmax =10;
 761
 762     Double_t m, p, pmax, f;
 763     m = kMassNeutron;
 764     pmax = TMath::Sqrt(2*T*(T+TMath::Sqrt(T*T+m*m)));
 765
 766     do
 767       {
 768         p = gRandom->Rndm() * fPmax;
 769         f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
 770       }
 771     while(f < gRandom->Rndm());
 772
 773     return 13.5*TMath::Sqrt(p*p+m*m);
 774 }
 775
 776 Double_t AliCentralityGlauberFit::Maxwell(Double_t m, Double_t p, Double_t T)
 777 {
 778 /* Relativistic Maxwell-distribution */
 779     Double_t ekin;
 780     ekin = TMath::Sqrt(p*p+m*m)-m;
 781     return (p*p * exp(-ekin/T));
 782 }
 783