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 /*   Origin: Alberica Toia, CERN, Alberica.Toia@cern.ch                   */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  class to determine centrality percentiles from 1D distributions          //
  21 //                                                                           //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24 #include <TNamed.h>
  25 #include <TH1D.h>
  26 #include <TString.h>
  27 #include <TFile.h>
  28 #include <TMath.h>
  29 #include <TROOT.h>
  30 #include <TH2F.h>
  31 #include <TF1.h>
  32 #include <TNtuple.h>
  33 #include <TStyle.h>
  34 #include <TGraphErrors.h>
  35 #include <vector>
  36 #include <TMinuit.h>
  37 #include <TStopwatch.h>
  38 #include <TRandom.h>
  39 #include <TDatabasePDG.h>
  40 #include <TPDGCode.h>
  41 #include "AliCentralityGlauberFit.h"
  42 #include "AliLog.h"
  43
  44 ClassImp(AliCentralityGlauberFit)
  45
  46 //______________________________________________________________________________
  47 AliCentralityGlauberFit::AliCentralityGlauberFit(const char *filename) :
  48   fNk(0),
  49   fKlow(0),
  50   fKhigh(0),
  51   fNalpha(0),
  52   fAlphalow(0),
  53   fAlphahigh(0),
  54   fNsigma(0),
  55   fSigmalow(0),
  56   fSigmahigh(0),
  57   fNbog(0),
  58   fBoglow(0),
  59   fBoghigh(0),
  60   fNCP(0),
  61   fCPlow(0),
  62   fCPhigh(0),
  63   fRebinFactor(0),
  64   fScalemin(0),
  65   fMultmin(0),
  66   fMultmax(0),
  67   fGlauntuple(0),
  68   fNpart(0),
  69   fNcoll(0),
  70   fB(0),
  71   fTaa(0),
  72   fTempHist(0),
  73   fGlauberHist(0),
  74   fUseChi2(kTRUE),
  75   fNevents(100000),
  76   fInrootfilename(0),
  77   fInntuplename(0),
  78   fOutrootfilename(0),
  79   fOutntuplename(0),
  80   fAncfilename("ancestor_hists.root"),
  81   fHistnames()
  82 {
  83   // Standard constructor.
  84   TFile *f = 0;
  85   if (filename) {  // glauber file
  86     f = TFile::Open(filename);
  87     fGlauntuple = (TNtuple*) f->Get("nt_p_Pb");
  88     fGlauntuple->SetBranchAddress("Npart",&fNpart);
  89     fGlauntuple->SetBranchAddress("Ncoll",&fNcoll);
  90     fGlauntuple->SetBranchAddress("B",&fB);
  91     fGlauntuple->SetBranchAddress("tAA",&fTaa);
  92   }
  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   NCP,
 128                                               Double_t CPlow,
 129                                               Double_t CPhigh)
 130 {
 131   // Set Glauber parameters.
 132   fNk=Nk;
 133   fKlow=klow;
 134   fKhigh=khigh;
 135   fNalpha=Nalpha;
 136   fAlphalow=alphalow;
 137   fAlphahigh=alphahigh;
 138   fNsigma=Nsigma;
 139   fSigmalow=sigmalow;
 140   fSigmahigh=sigmahigh;
 141   fNbog=Nbog;
 142   fBoglow=boglow;
 143   fBoghigh=boghigh;
 144   fNCP=NCP;
 145   fCPlow=CPlow;
 146   fCPhigh=CPhigh;
 147 }
 148
 149 //--------------------------------------------------------------------------------------------------
 150 void AliCentralityGlauberFit::MakeFits()
 151 {
 152   // Make fits.
 153
 154   TH1F *hDATA;
 155   TH1F *thistGlau;
 156   TFile *inrootfile;
 157   TFile *outrootfile;
 158   //FILE* fTxt = fopen ("parameters.txt","w");
 159
 160   // open inrootfile, outrootfile
 161   std::cout << "input file "  << fInrootfilename  << std::endl;
 162   std::cout << "output file " << fOutrootfilename << std::endl;
 163   inrootfile  = new TFile(fInrootfilename,"OPEN");
 164   outrootfile = new TFile(fOutrootfilename,"RECREATE");
 165
 166   // loop over all distribution names
 167   std::vector<TString>::const_iterator hni;
 168   for(hni=fHistnames.begin(); hni!=fHistnames.end(); hni++) {
 169     hDATA  = (TH1F*) (inrootfile->Get(*hni));
 170     if (!hDATA) {
 171       TList *list  = (TList*) (inrootfile->Get("CentralityStat"));
 172       //TList *list  = (TList*) (inrootfile->Get("VZEROEquaFactorStat"));
 173       hDATA  = (TH1F*) (list->FindObject(*hni));
 174     }
 175     hDATA->Rebin(fRebinFactor);
 176     //TH1F *hGLAU = new TH1F("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
 177
 178     Double_t chi2min = 9999999.0;
 179     Double_t CP_min=-1;
 180     Double_t bog_min=-1;
 181     Double_t sigma_min=-1;
 182     Double_t alpha_min=-1;
 183     Double_t k_min=-1;
 184     Double_t alpha, k, sigma, bog, CP, chi2;
 185
 186     for (Int_t nCP=0;nCP<fNCP;nCP++) {
 187       CP = fCPlow   + ((Double_t) nCP ) * (fCPhigh  - fCPlow ) / fNCP;
 188
 189       for (Int_t nbog=0;nbog<fNbog;nbog++) {
 190         bog = fBoglow   + ((Double_t) nbog ) * (fBoghigh  - fBoglow ) / fNbog;
 191
 192         for (Int_t nsigma=0;nsigma<fNsigma;nsigma++) {
 193           sigma = fSigmalow   + ((Double_t) nsigma ) * (fSigmahigh  - fSigmalow ) / fNsigma;
 194
 195           for (Int_t nalpha=0;nalpha<fNalpha;nalpha++) {
 196             alpha = fAlphalow   + ((Double_t) nalpha ) * (fAlphahigh  - fAlphalow ) / fNalpha;
 197
 198             for (Int_t nk=0; nk<fNk; nk++) {
 199               k = fKlow  + ((Double_t) nk ) * (fKhigh  - fKlow ) / fNk;
 200
 201               thistGlau = GlauberHisto(k,alpha,sigma,bog,CP,hDATA,kFALSE);
 202               chi2 = CalculateChi2(hDATA,thistGlau);
 203               //cout << "chi2 " << chi2<< endl;
 204               if ( chi2 < chi2min ) {
 205                 chi2min = chi2;
 206                 alpha_min=alpha;
 207                 sigma_min=sigma;
 208                 bog_min=bog;
 209                 CP_min=CP;
 210                 k_min=k;
 211               }
 212             }
 213           }
 214         }
 215       }
 216     }
 217     thistGlau = GlauberHisto(k_min,alpha_min,sigma_min,bog_min,CP_min, hDATA,kTRUE);
 218
 219     TH1F * hGLAU = 0x0;
 220     hGLAU = (TH1F *) thistGlau->Clone("hGLAU");
 221     hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU")));
 222     hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f_%.3f_%.3f",
 223                                                              k_min,alpha_min,sigma_min,bog_min,CP_min)));
 224
 225     Double_t mcintegral = hGLAU->Integral(hGLAU->FindBin(fScalemin),hGLAU->GetNbinsX());
 226     Double_t scale = (hDATA->Integral(hDATA->FindBin(fScalemin),hDATA->GetNbinsX())/mcintegral);
 227     std::cout << hGLAU->FindBin(fScalemin) << " " << hGLAU->Integral(hGLAU->FindBin(fScalemin),hGLAU->GetNbinsX()) << std::endl;
 228     std::cout << hDATA->FindBin(fScalemin) << " " << hDATA->Integral(hDATA->FindBin(fScalemin),hDATA->GetNbinsX()) << std::endl;
 229     hGLAU->Scale(scale);
 230
 231     SaveHisto(hDATA,hGLAU,outrootfile);
 232     //fclose (fTxt);
 233
 234     std::cout << "chi2 min is " << chi2min << std::endl;
 235     std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin),
 236                                               hGLAU->FindBin(fMultmax))/hGLAU->Integral()
 237               << " of the total cross section" << std::endl;
 238     fTempHist=hDATA;
 239     fGlauberHist=hGLAU;
 240   }
 241
 242   // close inrootfile, outrootfile
 243   inrootfile->Close();
 244   outrootfile->Close();
 245 }
 246
 247
 248 //--------------------------------------------------------------------------------------------------
 249 TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t k, Double_t alpha, Double_t sigma, Double_t bog, Double_t CP,
 250                                             TH1F *hDATA, Bool_t save)
 251 {
 252   // Get Glauber histogram.
 253   static TH1F *h1 = (TH1F*)hDATA->Clone();
 254   h1->Reset();
 255   h1->SetName(Form("fit_%.3f_%.3f",k,alpha));
 256
 257   TFile *outFile = NULL;
 258   TNtuple *ntuple = NULL;
 259
 260   if (save) {
 261     outFile = new TFile(fOutntuplename,"RECREATE");
 262     ntuple = new TNtuple("gnt", "Glauber ntuple", "Npart:Ncoll:B:tAA:ntot:nbn:ngn");
 263   }
 264
 265   Int_t nents = 0;
 266   if (fGlauntuple)
 267     nents = fGlauntuple->GetEntries();
 268
 269   for (Int_t i=0;i<fNevents;++i) {
 270     if (fGlauntuple)
 271       fGlauntuple->GetEntry(i % nents);
 272     else {
 273       fNpart = 2;
 274       fNcoll = 1;
 275     }
 276
 277     // Slow Nucleon Model from Chiara
 278     Double_t n;
 279     Double_t ntot, nbn, ngn, nbp, ngp;
 280     MakeSlowNucleons2(fNcoll,alpha,k,bog,CP, nbn,ngn,nbp,ngp);
 281
 282
 283     //conversion to energy *13.5; // conversion to energy------------------------
 284     //Double_t  Engn, Enbn;
 285     //Engn = ConvertToEnergy(0.05);
 286     //Enbn = ConvertToEnergy(0.005);
 287
 288     //Engn = 1;//gRandom->Gaus(13.5, TMath::Sqrt(13.5));
 289     //Enbn = 1;//gRandom->Gaus(13.5, TMath::Sqrt(13.5));
 290
 291     //Enbn =(int)nbn*Enbn;
 292     //Engn =(int)ngn*Engn;
 293
 294     //Enbn=0;
 295     //for (int i=0; i<(int)nbn; i++){
 296     // Enbn += ConvertToEnergy(0.005);
 297     //}
 298     //for (int i=0; i<(int)ngn; i++){
 299     //  Engn += ConvertToEnergy(0.05);
 300     // }
 301
 302     //cout << "GRAY " << Engn << " BLACK " << Enbn << endl;
 303
 304     // acceptance
 305     n=alpha*ngn+alpha*nbn;
 306     // n=alpha*ngp+alpha*nbp;
 307     //----------------------------------------
 308
 309     //------ experimental resolution -> Gaussian smearing ------------------------------------
 310     Double_t resexp=0;
 311     //if (n>0) resexp = 1./TMath::Sqrt(n)*sigma*gRandom->Gaus(0.0,1.0);
 312     //else resexp=0;
 313     //ntot=n*(1+resexp);
 314
 315     if (n>0)resexp = sigma*TMath::Sqrt(n)/2;
 316     else resexp=0;
 317     ntot = (Int_t)(gRandom->Gaus(n,resexp));
 318     //----------------------------------------
 319
 320     // non-lineary effect -------------------------------
 321     //ntot = k*ntot;
 322     ntot = ntot + k*ntot*ntot;
 323
 324     //cout << ngn << " " << nbn << " "  << ntot << endl;
 325
 326     if (ntot>0) {
 327       h1->Fill(ntot);
 328       if (save)
 329         ntuple->Fill(fNpart,fNcoll,fB,fTaa,ntot,nbn,ngn);
 330     }
 331   }
 332
 333   if (save) {
 334     ntuple->Write();
 335     outFile->Close();
 336   }
 337   return h1;
 338 }
 339
 340 //--------------------------------------------------------------------------------------------------
 341 Double_t AliCentralityGlauberFit::CalculateChi2(TH1F *hDATA, TH1F *thistGlau)
 342 {
 343   // Note that for different values of neff the mc distribution is identical, just re-scaled below
 344   // normalize the two histogram, scale MC up but leave data alone - that preserves error bars !!!!
 345
 346   Int_t lowchibin =   hDATA->FindBin(fMultmin);
 347   Int_t highchibin =  hDATA->FindBin(fMultmax);
 348
 349   Double_t mcintegral = thistGlau->Integral(1,thistGlau->GetNbinsX());
 350   Double_t scale = (hDATA->Integral(1,hDATA->GetNbinsX())/mcintegral);
 351   thistGlau->Scale(scale);
 352
 353   // calculate the chi2 between MC and real data over some range ????
 354   if (fUseChi2) {
 355     Double_t chi2 = 0.0;
 356     for (Int_t i=lowchibin; i<=highchibin; i++) {
 357       if (hDATA->GetBinContent(i) < 1.0) continue;
 358       Double_t diff = TMath::Power((thistGlau->GetBinContent(i) - hDATA->GetBinContent(i)),2);
 359       diff = diff / (TMath::Power(hDATA->GetBinError(i),2) + TMath::Power(thistGlau->GetBinError(i),2));
 360       chi2 += diff;
 361     }
 362     chi2 = chi2 / (highchibin - lowchibin + 1);
 363     return chi2;
 364   }
 365   // "-2 log likelihood ratio(mu;n) = 2[(mu - n) + n ln(n/mu)]"
 366   else {
 367     std::cout << "LL" << std::endl;
 368
 369     Double_t ll = 0.0;
 370     for (Int_t i=lowchibin; i<=highchibin; i++) {
 371       Double_t data = hDATA    ->GetBinContent(i);
 372       Double_t mc   = thistGlau->GetBinContent(i);
 373       Int_t idata = TMath::Nint(data);
 374       if (mc < 1.e-9) mc = 1.e-9;
 375       Double_t fsub = - mc + idata * TMath::Log(mc);
 376       Double_t fobs = 0;
 377       if (idata > 0) {
 378         for(Int_t istep = 0; istep < idata; istep++){
 379           if (istep > 1)
 380             fobs += TMath::Log(istep);
 381         }
 382       }
 383       fsub -= fobs;
 384       ll -=  fsub ;
 385     }
 386     return 2*ll;
 387   }
 388 }
 389
 390 //--------------------------------------------------------------------------------------------------
 391 void AliCentralityGlauberFit::SaveHisto(TH1F *hist1, TH1F *hist2, TFile *outrootfile)
 392 {
 393   // Save histograms.
 394
 395   outrootfile->cd();
 396   hist1->Write();
 397   hist2->Write();
 398 }
 399
 400 //--------------------------------------------------------------------------------------------------
 401 void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t alpha, Double_t k, Double_t &nbn, Double_t &ngn)
 402 {
 403 // from AliGenSlowNucleonModelExp (Chiara Oppedisano)
 404 // Return the number of black and gray nucleons
 405 //
 406 // Number of collisions
 407
 408   Int_t fP = 82;
 409   Int_t fN = 126;
 410   Float_t nu = (Float_t) ncoll;
 411   //
 412 // Mean number of gray nucleons
 413     Float_t nGray         = alpha * nu;
 414     Float_t nGrayNeutrons = nGray * fN / (fN + fP);
 415     Float_t nGrayProtons  = nGray - nGrayNeutrons;
 416
 417 // Mean number of black nucleons
 418     Float_t nBlack  = 0.;
 419     if(!fApplySaturation || (fApplySaturation && nGray<fnGraySaturation)) nBlack = k * nu;
 420     else if(fApplySaturation && nGray>=fnGraySaturation) nBlack = fnBlackSaturation;
 421     Float_t nBlackNeutrons = nBlack * 0.84;
 422     Float_t nBlackProtons  = nBlack - nBlackNeutrons;
 423
 424 // Actual number (including fluctuations) from binomial distribution
 425     Double_t p, ngp, nbp;
 426
 427 //  gray neutrons
 428     p =  nGrayNeutrons/fN;
 429     ngn = gRandom->Binomial((Int_t) fN, p);
 430
 431 //  gray protons
 432     p =  nGrayProtons/fP;
 433     ngp = gRandom->Binomial((Int_t) fP, p);
 434
 435 //  black neutrons
 436     p =  nBlackNeutrons/fN;
 437     nbn = gRandom->Binomial((Int_t) fN, p);
 438
 439 //  black protons
 440     p =  nBlackProtons/fP;
 441     nbp = gRandom->Binomial((Int_t) fP, p);
 442
 443 }
 444
 445
 446
 447 //--------------------------------------------------------------------------------------------------
 448  void AliCentralityGlauberFit::MakeSlowNucleons2(Int_t ncoll, Double_t alpha, Double_t k, Double_t bog, Double_t CP, Double_t &nbn, Double_t &ngn, Double_t &nbp, Double_t &ngp)
 449 //void AliCentralityGlauberFit::MakeSlowNucleons2(Int_t ncoll, Double_t alpha, Double_t k, Double_t bog, Double_t CP, Double_t &nbn, Double_t &ngn)
 450 {
 451 // from AliGenSlowNucleonModelExp (Chiara Oppedisano)
 452 //
 453 // Return the number of black and gray nucleons
 454 //
 455 // Number of collisions
 456
 457    // based on E910 model ================================================================
 458
 459   Int_t fP = 82;
 460   Int_t fN = 126;
 461
 462   Float_t nu = (Float_t) ncoll;
 463   //
 464   nu = nu+1.*gRandom->Rndm();
 465   //
 466
 467   //  gray protons
 468   Float_t  poverpd = 0.843;
 469   Float_t  zAu2zPb = 82./79.;
 470   Float_t  nGrayp = (-0.27 + 0.63 * nu - 0.0008 *nu *nu)*poverpd*zAu2zPb;
 471
 472   Double_t p;
 473   p =  nGrayp/fP;
 474   ngp = gRandom->Binomial((Int_t) fP, p);
 475   //ngp = gRandom->Gaus(nGrayp, TMath::Sqrt(fP*p*(1-p)));
 476   if(nGrayp<0.) ngp=0;
 477
 478   //  black protons
 479   //Float_t blackovergray = 3./7.;// from spallation
 480   //Float_t blackovergray = 0.65; // from COSY
 481   Float_t blackovergray = bog;
 482   Float_t nBlackp  = blackovergray*nGrayp;
 483
 484   p =  nBlackp/fP;
 485   nbp = gRandom->Binomial((Int_t) fP, p);
 486   //nbp = gRandom->Gaus(nBlackp, TMath::Sqrt(fP*p*(1-p)));
 487   if(nBlackp<0.) nbp=0;
 488
 489   if(nu<3.){
 490     nGrayp = -0.836 + 0.9112 *nu - 0.05381 *nu *nu;
 491     nBlackp  = blackovergray*nGrayp;
 492   }
 493
 494
 495   //  gray neutrons
 496   Float_t nGrayNeutrons = 0.;
 497   Float_t nBlackNeutrons = 0.;
 498   //Float_t cp = (nGrayp+nBlackp)/0.24;
 499   Float_t cp = (nGrayp+nBlackp)/CP;
 500
 501   // if(cp>0.){
 502   //   //Float_t nSlow  = 51.5+469.2/(-8.762-cp);
 503   //   Float_t nSlow  = 60.0+469.2/(-8.762-cp);
 504   //   //if(cp<2.5) nSlow = 1+(9.9-1)/(2.5-0)*(cp-0);
 505   //   if(cp<3.) nSlow = 0.+(11.6-0.)/(3.-0.)*(cp-0.);
 506
 507   //   nGrayNeutrons = nSlow * 0.1;
 508   //   nBlackNeutrons = nSlow - nGrayNeutrons;
 509   // }
 510   if(cp>0.){
 511     Float_t paramnSlow[3] = {60., 469.2, 8.762};
 512     Float_t nSlow  = paramnSlow[0]+paramnSlow[1]/(-paramnSlow[2]-cp);
 513     float paramRetta = paramnSlow[0]+paramnSlow[1]/(-paramnSlow[2]-3);
 514     if(cp<3.) nSlow = 0.+(paramRetta-0.)/(3.-0.)*(cp-0.);
 515
 516     nGrayNeutrons = nSlow * 0.1;
 517     nBlackNeutrons = nSlow - nGrayNeutrons;
 518   }
 519   else{
 520     // Sikler "pasturato"
 521     nGrayNeutrons = 0.47 * 2.3 *  nu; // fAlphaGray=2.3
 522     nBlackNeutrons = 0.88 * 3.6 * nu; // fAlphaBlack=3.6
 523     printf("nslowp=0 -> ncoll = %1.0f -> ngrayn = %1.0f  nblackn = %1.0f \n", nu, nGrayNeutrons, nBlackNeutrons);
 524   }
 525
 526   p =  nGrayNeutrons/fN;
 527   //ngn = gRandom->Binomial((Int_t) fN, p);
 528   ngn = gRandom->Gaus(nGrayNeutrons, TMath::Sqrt(fN*p*(1-p)));
 529
 530   //  black neutrons
 531   p =  nBlackNeutrons/fN;
 532   //nbn = gRandom->Binomial((Int_t) fN, p);
 533   nbn = gRandom->Gaus(nBlackNeutrons, TMath::Sqrt(fN*p*(1-p)));
 534 }
 535
 536 Double_t AliCentralityGlauberFit::ConvertToEnergy(Double_t T)
 537 {
 538    TDatabasePDG * pdg = TDatabasePDG::Instance();
 539     const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
 540     Double_t fPmax =10;
 541
 542     Double_t m, p, pmax, ekin, f;
 543     m = kMassNeutron;
 544     pmax = TMath::Sqrt(2*T*(T+TMath::Sqrt(T*T+m*m)));
 545
 546     do
 547       {
 548         p = gRandom->Rndm() * fPmax;
 549         f = Maxwell(m, p, T) / Maxwell(m , pmax, T);
 550       }
 551     while(f < gRandom->Rndm());
 552
 553     return 13.5*TMath::Sqrt(p*p+m*m);
 554 }
 555
 556 Double_t AliCentralityGlauberFit::Maxwell(Double_t m, Double_t p, Double_t T)
 557 {
 558 /* Relativistic Maxwell-distribution */
 559     Double_t ekin;
 560     ekin = TMath::Sqrt(p*p+m*m)-m;
 561     return (p*p * exp(-ekin/T));
 562 }
 563