X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=FASTSIM%2FAliFastGlauber.cxx;h=88089cb944339039a9727aa2b8b3708cf1d2a1d0;hb=a8bdf7d2fd4d224a53ebe00eb8d20eb2a8e6857f;hp=201207b0b698d3badab84ed131ee190c783729b1;hpb=710a8d903cb0630964c4ef2263574ab201afaa5a;p=u%2Fmrichter%2FAliRoot.git diff --git a/FASTSIM/AliFastGlauber.cxx b/FASTSIM/AliFastGlauber.cxx index 201207b0b69..88089cb9443 100644 --- a/FASTSIM/AliFastGlauber.cxx +++ b/FASTSIM/AliFastGlauber.cxx @@ -16,7 +16,7 @@ /* $Id$ */ // // Utility class to make simple Glauber type calculations -// for SYMMTRIC collision geometries (AA): +// for SYMMETRIC collision geometries (AA): // Impact parameter, production points, reaction plane dependence // // The SimulateTrigger method can be used for simple MB and hard-process @@ -50,17 +50,16 @@ // from AliRoot #include "AliFastGlauber.h" // from root -#include -#include -#include +#include #include #include +#include +#include +#include #include -#include +#include #include -#include -#include -#include +#include ClassImp(AliFastGlauber) @@ -82,28 +81,68 @@ TF1* AliFastGlauber::fgWIntRadius = NULL; TF2* AliFastGlauber::fgWKParticipants = NULL; TF1* AliFastGlauber::fgWParticipants = NULL; TF2* AliFastGlauber::fgWAlmondCurrent = NULL; -TF2 AliFastGlauber::fgWAlmondFixedB[40]; +TF2* AliFastGlauber::fgWAlmondFixedB[40]; const Int_t AliFastGlauber::fgkMCInts = 100000; -Int_t AliFastGlauber::fgCounter = 0; - -AliFastGlauber::AliFastGlauber() : fName() -{ - // Default Constructor - // - fgCounter++; - if(fgCounter>1) - Error("AliFastGlauber","More than more instance (%d) is not supported, check your code!",fgCounter); - +AliFastGlauber* AliFastGlauber::fgGlauber = NULL; + + +AliFastGlauber::AliFastGlauber(): + fWSr0(0.), + fWSd(0.), + fWSw(0.), + fWSn(0.), + fSigmaHard(0.), + fSigmaNN(0.), + fA(0), + fBmin(0.), + fBmax(0.), + fEllDef(0), + fName() +{ + // Default Constructor // Defaults for Pb SetMaxImpact(); SetLengthDefinition(); SetPbPbLHC(); + fXY[0] = fXY[1] = 0; + fI0I1[0] = fI0I1[1] = 0; +} + +AliFastGlauber::AliFastGlauber(const AliFastGlauber & gl) + :TObject(gl), + fWSr0(0.), + fWSd(0.), + fWSw(0.), + fWSn(0.), + fSigmaHard(0.), + fSigmaNN(0.), + fA(0), + fBmin(0.), + fBmax(0.), + fEllDef(0), + fName() +{ +// Copy constructor + gl.Copy(*this); + fXY[0] = fXY[1] = 0; + fI0I1[0] = fI0I1[1] = 0; +} + +AliFastGlauber* AliFastGlauber::Instance() +{ +// Set random number generator + if (fgGlauber) { + return fgGlauber; + } else { + fgGlauber = new AliFastGlauber(); + return fgGlauber; + } } AliFastGlauber::~AliFastGlauber() { - fgCounter--; - //if(fgCounter==0) Reset(); +// Destructor + for(Int_t k=0; k<40; k++) delete fgWAlmondFixedB[k]; } void AliFastGlauber::SetAuAuRhic() @@ -148,7 +187,7 @@ void AliFastGlauber::Init(Int_t mode) fgWSb->SetParameter(2, fWSw); fgWSb->SetParameter(3, fWSn); - fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 4); + fgWSbz = new TF2("WSbz", WSbz, 0, fgBMax, 0, fgBMax, 4); fgWSbz->SetParameter(0, fWSr0); fgWSbz->SetParameter(1, fWSd); fgWSbz->SetParameter(2, fWSw); @@ -247,7 +286,7 @@ void AliFastGlauber::Init(Int_t mode) for(Int_t k=0; k<40; k++) { sprintf(almondName,"WAlmondFixedB%d",k); fgWAlmondCurrent = (TF2*)ff->Get(almondName); - new(&fgWAlmondFixedB[k]) TF2(*fgWAlmondCurrent); + fgWAlmondFixedB[k] = fgWAlmondCurrent; } delete ff; } @@ -268,7 +307,7 @@ void AliFastGlauber::Init(Int_t mode) fgWPathLength->SetParameter(2, 0); //Pathlength definition } -void AliFastGlauber::Reset() +void AliFastGlauber::Reset() const { // // Reset dynamic allocated formulas @@ -331,7 +370,7 @@ void AliFastGlauber::DrawOverlap() const TCanvas *c2 = new TCanvas("c2","Overlap",400,10,600,700); c2->cd(); Double_t max=fgWStaa->GetMaximum(0,fgBMax)*1.01; - TH2F *h2f=new TH2F("h2ftaa","Overlap function: T_{AB} [mbarn^{-1}]",2,0,fgBMax,2,0,max); + TH2F *h2f=new TH2F("h2ftaa","Overlap function: T_{AB} [mbarn^{-1}]",2,0,fgBMax,2,0, max); h2f->SetStats(0); h2f->GetXaxis()->SetTitle("b [fm]"); h2f->GetYaxis()->SetTitle("T_{AB} [mbarn^{-1}]"); @@ -471,7 +510,7 @@ void AliFastGlauber::DrawKernel(Double_t b) const l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; - sprintf(label,"b = %.f fm",b); + sprintf(label,"b = %.1f fm",b); l1a->AddEntry(fgWStarfi,label,""); l1a->Draw(); c8->Update(); @@ -495,7 +534,7 @@ void AliFastGlauber::DrawAlmond(Double_t b) const l1a->SetFillStyle(0); l1a->SetBorderSize(0); Char_t label[100]; - sprintf(label,"b = %.f fm",b); + sprintf(label,"b = %.1f fm",b); l1a->AddEntry(fgWAlmond,label,""); l1a->Draw(); c9->Update(); @@ -577,7 +616,7 @@ void AliFastGlauber::DrawIntRadius(Double_t b) const fgWIntRadius->Draw("same"); } -Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par) +Double_t AliFastGlauber::WSb(const Double_t* x, const Double_t* par) { // // Woods-Saxon Parameterisation @@ -592,7 +631,7 @@ Double_t AliFastGlauber::WSb(Double_t* x, Double_t* par) return y; //fm^-3 } -Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par) +Double_t AliFastGlauber::WSbz(const Double_t* x, const Double_t* par) { // // Wood Saxon Parameterisation @@ -609,7 +648,7 @@ Double_t AliFastGlauber::WSbz(Double_t* x, Double_t* par) return y; //fm^-3 } -Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par) +Double_t AliFastGlauber::WSz(const Double_t* x, const Double_t* par) { // // Wood Saxon Parameterisation @@ -626,7 +665,7 @@ Double_t AliFastGlauber::WSz(Double_t* x, Double_t* par) return y; //fm^-3 } -Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/) +Double_t AliFastGlauber::WSta(const Double_t* x, const Double_t* /*par*/) { // // Thickness function T_A @@ -638,7 +677,7 @@ Double_t AliFastGlauber::WSta(Double_t* x, Double_t* /*par*/) return y; //fm^-2 } -Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par) +Double_t AliFastGlauber::WStarfi(const Double_t* x, const Double_t* par) { // // Kernel for overlap function: T_A(s)*T_A(s-b) @@ -651,7 +690,7 @@ Double_t AliFastGlauber::WStarfi(Double_t* x, Double_t* par) return y; //fm^-3 } -Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par) +Double_t AliFastGlauber::WStaa(const Double_t* x, const Double_t* par) { // // Overlap function @@ -681,8 +720,11 @@ Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par) // MC Integration // Double_t y = 0; + + for (Int_t i = 0; i < fgkMCInts; i++) { + const Double_t kphi = TMath::Pi() * gRandom->Rndm(); const Double_t kb1 = fgBMax * gRandom->Rndm(); y += fgWStarfi->Eval(kb1, kphi); @@ -692,7 +734,7 @@ Double_t AliFastGlauber::WStaa(Double_t* x, Double_t* par) return y; } -Double_t AliFastGlauber::WKParticipants(Double_t* x, Double_t* par) +Double_t AliFastGlauber::WKParticipants(const Double_t* x, const Double_t* par) { // // Kernel for number of participants @@ -721,7 +763,7 @@ Double_t AliFastGlauber::WKParticipants(Double_t* x, Double_t* par) return y; //fm^-1 } -Double_t AliFastGlauber::WParticipants(Double_t* x, Double_t* par) +Double_t AliFastGlauber::WParticipants(const Double_t* x, const Double_t* par) { // // Number of Participants as @@ -748,7 +790,7 @@ Double_t AliFastGlauber::WParticipants(Double_t* x, Double_t* par) return y; //no units } -Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* par) +Double_t AliFastGlauber::WSgeo(const Double_t* x, const Double_t* par) { // // Geometrical Cross-Section @@ -761,7 +803,7 @@ Double_t AliFastGlauber::WSgeo(Double_t* x, Double_t* par) return y; //fm } -Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par) +Double_t AliFastGlauber::WSbinary(const Double_t* x, const Double_t* par) { // // Number of binary hard collisions @@ -774,7 +816,7 @@ Double_t AliFastGlauber::WSbinary(Double_t* x, Double_t* par) return y; //fm } -Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/) +Double_t AliFastGlauber::WSN(const Double_t* x, const Double_t* /*par*/) { // // Number of hard processes per event @@ -784,7 +826,7 @@ Double_t AliFastGlauber::WSN(Double_t* x, Double_t* /*par*/) return y; //no units } -Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par) +Double_t AliFastGlauber::WEnergyDensity(const Double_t* x, const Double_t* par) { // // Initial energy density @@ -803,7 +845,7 @@ Double_t AliFastGlauber::WEnergyDensity(Double_t* x, Double_t* par) return y; //fm^-4 } -Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par) +Double_t AliFastGlauber::WAlmond(const Double_t* x, const Double_t* par) { // // Almond shaped interaction region @@ -822,7 +864,7 @@ Double_t AliFastGlauber::WAlmond(Double_t* x, Double_t* par) return y; //fm^-4 } -Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par) +Double_t AliFastGlauber::WIntRadius(const Double_t* x, const Double_t* par) { // // Average interaction density over radius @@ -847,7 +889,7 @@ Double_t AliFastGlauber::WIntRadius(Double_t* x, Double_t* par) return y; //fm^-3 } -Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par) +Double_t AliFastGlauber::WPathLength0(const Double_t* x, const Double_t* par) { // // Path Length as a function of phi @@ -884,16 +926,16 @@ Double_t AliFastGlauber::WPathLength0(Double_t* x, Double_t* par) } // radial steps Double_t y=0.; - if (!kiopt) // My length definition (is exact for hard disk) - if(w) y= 2. * rw / w; - else { - const Double_t knorm=fgWSta->Eval(1e-4); - if(knorm) y = TMath::Sqrt(2. * rw * kDr / knorm / knorm); + if (!kiopt) { // My length definition (is exact for hard disk) + if(w) y= 2. * rw / w; + } else { + const Double_t knorm=fgWSta->Eval(1e-4); + if(knorm) y = TMath::Sqrt(2. * rw * kDr / knorm / knorm); } return y; //fm } -Double_t AliFastGlauber::WPathLength(Double_t* x, Double_t* par) +Double_t AliFastGlauber::WPathLength(const Double_t* x, const Double_t* par) { // // Path Length as a function of phi @@ -972,26 +1014,77 @@ Double_t AliFastGlauber::HardCrossSection(Double_t b1, Double_t b2) const Double_t AliFastGlauber::FractionOfHardCrossSection(Double_t b1, Double_t b2) const { // - // Return raction of hard cross-section integrated from b1 to b2 + // Return fraction of hard cross-section integrated from b1 to b2 // return fgWSbinary->Integral(b1, b2)/fgWSbinary->Integral(0., 100.); } +Double_t AliFastGlauber::NHard(const Double_t b1, const Double_t b2) const +{ + // + // Number of binary hard collisions + // as a function of b (nucl/ex/0302016 eq. 19) + // + const Double_t kshard=HardCrossSection(b1,b2); + const Double_t ksgeo=CrossSection(b1,b2); + if(ksgeo>0) + return kshard/ksgeo; + else return -1; +} + Double_t AliFastGlauber::Binaries(Double_t b) const { // // Return number of binary hard collisions normalized to 1 at b=0 // - if(b==0) b=1e-4; + if(b < 1.e-4) b = 1e-4; return fgWSN->Eval(b)/fgWSN->Eval(1e-4); } +Double_t AliFastGlauber::MeanOverlap(Double_t b1, Double_t b2) +{ +// +// Calculate the mean overlap for impact parameter range b1 .. b2 +// + Double_t sum = 0.; + Double_t sumc = 0.; + Double_t b = b1; + + while (b < b2-0.005) { + Double_t nc = GetNumberOfCollisions(b); + sum += 10. * fgWStaa->Eval(b) * fgWSgeo->Eval(b) * 0.01 / (1. - TMath::Exp(-nc)); + sumc += 10. * fgWSgeo->Eval(b) * 0.01; + b += 0.01; + } + return (sum / CrossSection(b1, b2)); +} + + +Double_t AliFastGlauber::MeanNumberOfCollisionsPerEvent(Double_t b1, Double_t b2) +{ +// +// Calculate the mean number of collisions per event for impact parameter range b1 .. b2 +// + Double_t sum = 0.; + Double_t sumc = 0.; + Double_t b = b1; + + while (b < b2-0.005) { + Double_t nc = GetNumberOfCollisions(b); + sum += nc / (1. - TMath::Exp(-nc)) * 10. * fgWSgeo->Eval(b) * 0.01; + sumc += 10. * fgWSgeo->Eval(b) * 0.01; + b += 0.01; + } + return (sum / CrossSection(b1, b2)); +} + + Double_t AliFastGlauber::GetNumberOfBinaries(Double_t b) const { // // Return number of binary hard collisions at b // - if(b==0) b=1e-4; + if(b<1.e-4) b=1e-4; return fgWSN->Eval(b); } @@ -1000,7 +1093,7 @@ Double_t AliFastGlauber::Participants(Double_t b) const // // Return the number of participants normalized to 1 at b=0 // - if(b==0) b=1e-4; + if(b<1.e-4) b=1e-4; return (fgWParticipants->Eval(b)/fgWParticipants->Eval(1e-4)); } @@ -1009,7 +1102,7 @@ Double_t AliFastGlauber::GetNumberOfParticipants(Double_t b) const // // Return the number of participants for impact parameter b // - if(b==0) b=1e-4; + if(b<1.e-4) b=1e-4; return (fgWParticipants->Eval(b)); } @@ -1018,10 +1111,24 @@ Double_t AliFastGlauber::GetNumberOfCollisions(Double_t b) const // // Return the number of collisions for impact parameter b // - if(b==0) b=1e-4; + if(b<1.e-4) b=1e-4; return (fgWStaa->Eval(b)*fSigmaNN); } +Double_t AliFastGlauber::GetNumberOfCollisionsPerEvent(Double_t b) const +{ + // + // Return the number of collisions per event (at least one collision) + // for impact parameter b + // + Double_t n = GetNumberOfCollisions(b); + if (n > 0.) { + return (n / (1. - TMath::Exp(- n))); + } else { + return (0.); + } +} + void AliFastGlauber::SimulateTrigger(Int_t n) { // @@ -1190,7 +1297,7 @@ void AliFastGlauber::GetRandomBHard(Double_t& b) b = fgWSbinary->GetRandom(); Int_t bin = 2*(Int_t)b; if( (b-(Int_t)b) > 0.5) bin++; - fgWAlmondCurrent = &fgWAlmondFixedB[bin]; + fgWAlmondCurrent = fgWAlmondFixedB[bin]; return; } @@ -1295,34 +1402,59 @@ Double_t AliFastGlauber::CalculateLength(Double_t b,Double_t x0,Double_t y0,Doub } } -void AliFastGlauber::GetLength(Double_t& ell,Double_t b) +void AliFastGlauber::GetLengthAndPhi(Double_t& ell,Double_t& phi,Double_t b) { // // Return length from random b, x0, y0, phi0 + // Return also phi0 // Double_t x0,y0,phi0; if(b<0.) GetRandomBHard(b); GetRandomXY(x0,y0); GetRandomPhi(phi0); + phi = phi0; ell = CalculateLength(b,x0,y0,phi0); return; } -void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2,Double_t b) +void AliFastGlauber::GetLength(Double_t& ell,Double_t b) { // - // Return 2 lengths back to back from random b, x0, y0, phi0 + // Return length from random b, x0, y0, phi0 + // + Double_t phi; + GetLengthAndPhi(ell,phi,b); + return; +} + +void AliFastGlauber::GetLengthsBackToBackAndPhi(Double_t& ell1,Double_t& ell2,Double_t &phi,Double_t b) +{ // + // Return 2 lengths back to back from random b, x0, y0, phi0 + // Return also phi0 + // Double_t x0,y0,phi0; if(b<0.) GetRandomBHard(b); GetRandomXY(x0,y0); GetRandomPhi(phi0); const Double_t kphi0plusPi = phi0+TMath::Pi(); + phi = phi0; ell1 = CalculateLength(b,x0,y0,phi0); ell2 = CalculateLength(b,x0,y0,kphi0plusPi); return; } +void AliFastGlauber::GetLengthsBackToBack(Double_t& ell1,Double_t& ell2, + Double_t b) +{ + // + // Return 2 lengths back to back from random b, x0, y0, phi0 + // + Double_t phi; + GetLengthsBackToBackAndPhi(ell1,ell2,phi,b); + return; +} + void AliFastGlauber::GetLengthsForPythia(Int_t n,Double_t* phi,Double_t* ell, Double_t b) { // @@ -1356,7 +1488,7 @@ void AliFastGlauber::PlotBDistr(Int_t n) return; } -void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,Char_t *fname) +void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,const char *fname) { // // Plot length distribution @@ -1383,7 +1515,7 @@ void AliFastGlauber::PlotLengthDistr(Int_t n,Bool_t save,Char_t *fname) return; } -void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,Char_t *fname) +void AliFastGlauber::PlotLengthB2BDistr(Int_t n,Bool_t save,const char *fname) { // // Plot lengths back-to-back distributions @@ -1418,13 +1550,13 @@ void AliFastGlauber::PlotAlmonds() const gStyle->SetPalette(1,0); c->Divide(2,2); c->cd(1); - fgWAlmondFixedB[0].Draw("cont1"); + fgWAlmondFixedB[0]->Draw("cont1"); c->cd(2); - fgWAlmondFixedB[10].Draw("cont1"); + fgWAlmondFixedB[10]->Draw("cont1"); c->cd(3); - fgWAlmondFixedB[20].Draw("cont1"); + fgWAlmondFixedB[20]->Draw("cont1"); c->cd(4); - fgWAlmondFixedB[30].Draw("cont1"); + fgWAlmondFixedB[30]->Draw("cont1"); return; } @@ -1436,8 +1568,8 @@ void AliFastGlauber::CalculateI0I1(Double_t& integral0,Double_t& integral1, { // // Calculate integrals: - // integral0 = \int_0^ellCut dl*(T_A*T_B)(x0+l*ux,y0+l*uy) - // integral1 = \int_0^ellCut dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy) + // integral0 = \int_0^ellCut dl*(T_A*T_B)(x0+l*ux,y0+l*uy) + // integral1 = \int_0^ellCut dl*l*(T_A*T_B)(x0+l*ux,y0+l*uy) // // for a parton with production point (x0,y0) // and propagation direction (ux=cos(phi0),uy=sin(phi0)) @@ -1474,37 +1606,87 @@ void AliFastGlauber::CalculateI0I1(Double_t& integral0,Double_t& integral1, return; } +void AliFastGlauber::GetI0I1AndPhi(Double_t& integral0,Double_t& integral1, + Double_t& phi, + Double_t ellCut,Double_t b) +{ + // + // Return I0 and I1 from random b, x0, y0, phi0 + // Return also phi + // + Double_t x0,y0,phi0; + if(b<0.) GetRandomBHard(b); + GetRandomXY(x0,y0); + GetRandomPhi(phi0); + phi = phi0; + CalculateI0I1(integral0,integral1,b,x0,y0,phi0,ellCut); + return; +} + void AliFastGlauber::GetI0I1(Double_t& integral0,Double_t& integral1, Double_t ellCut,Double_t b) { // // Return I0 and I1 from random b, x0, y0, phi0 // + Double_t phi; + GetI0I1AndPhi(integral0,integral1,phi,ellCut,b); + return; +} + +void AliFastGlauber::GetI0I1BackToBackAndPhi(Double_t& integral01,Double_t& integral11, + Double_t& integral02,Double_t& integral12, + Double_t& phi, + Double_t ellCut,Double_t b) +{ + // + // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0 + // Return also phi0 + // Double_t x0,y0,phi0; if(b<0.) GetRandomBHard(b); GetRandomXY(x0,y0); GetRandomPhi(phi0); - CalculateI0I1(integral0,integral1,b,x0,y0,phi0,ellCut); + phi = phi0; + const Double_t kphi0plusPi = phi0+TMath::Pi(); + CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut); + CalculateI0I1(integral02,integral12,b,x0,y0,kphi0plusPi,ellCut); return; } -void AliFastGlauber::GetI0I1BackToBack(Double_t& integral01,Double_t& integral11, - Double_t& integral02,Double_t& integral12, - Double_t ellCut,Double_t b) +void AliFastGlauber::GetI0I1BackToBackAndPhiAndXY(Double_t& integral01,Double_t& integral11, + Double_t& integral02,Double_t& integral12, + Double_t& phi,Double_t &x,Double_t &y, + Double_t ellCut,Double_t b) { // // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0 + // Return also phi0 // Double_t x0,y0,phi0; if(b<0.) GetRandomBHard(b); GetRandomXY(x0,y0); GetRandomPhi(phi0); + phi = phi0; x=x0; y=y0; const Double_t kphi0plusPi = phi0+TMath::Pi(); CalculateI0I1(integral01,integral11,b,x0,y0,phi0,ellCut); CalculateI0I1(integral02,integral12,b,x0,y0,kphi0plusPi,ellCut); return; } +void AliFastGlauber::GetI0I1BackToBack(Double_t& integral01,Double_t& integral11, + Double_t& integral02,Double_t& integral12, + Double_t ellCut,Double_t b) +{ + // + // Return 2 pairs of I0 and I1 back to back from random b, x0, y0, phi0 + // + Double_t phi; + GetI0I1BackToBackAndPhi(integral01,integral11,integral02,integral12, + phi,ellCut,b); + return; +} + void AliFastGlauber::GetI0I1ForPythia(Int_t n,Double_t* phi, Double_t* integral0,Double_t* integral1, Double_t ellCut,Double_t b) @@ -1521,20 +1703,43 @@ void AliFastGlauber::GetI0I1ForPythia(Int_t n,Double_t* phi, return; } +void AliFastGlauber::GetI0I1ForPythiaAndXY(Int_t n,Double_t* phi, + Double_t* integral0,Double_t* integral1, + Double_t &x,Double_t& y, + Double_t ellCut,Double_t b) +{ + // + // Returns I0 and I1 pairs for n partons with azimuthal angles phi[n] + // from random b, x0, y0 and return x0,y0 + // + Double_t x0,y0; + if(b<0.) GetRandomBHard(b); + GetRandomXY(x0,y0); + for(Int_t i=0; iSetXTitle("I_{0} [fm^{-3}]"); + hI0I1s->SetYTitle("I_{1} [fm^{-2}]"); + TH1F *hI0 = new TH1F("hI0","I_{0} = #hat{q}L / k", - 100,0,0.001); + 1000,0,0.001); hI0->SetXTitle("I_{0} [fm^{-3}]"); hI0->SetYTitle("Probability"); hI0->SetFillColor(3); TH1F *hI1 = new TH1F("hI1","I_{1} = #omega_{c} / k", - 100,0,0.01); + 1000,0,0.01); hI1->SetXTitle("I_{1} [fm^{-2}]"); hI1->SetYTitle("Probability"); hI1->SetFillColor(4); @@ -1556,6 +1761,7 @@ void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut, for(Int_t i=0; iFill(i0,i1); hI0->Fill(i0); hI1->Fill(i1); h2->Fill(2.*i1*i1/i0); @@ -1564,6 +1770,10 @@ void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut, } hI0->Scale(1/(Double_t)n); hI1->Scale(1/(Double_t)n); + h2->Scale(1/(Double_t)n); + h3->Scale(1/(Double_t)n); + h4->Scale(1/(Double_t)n); + hI0I1s->Scale(1/(Double_t)n); TCanvas *cI0I1 = new TCanvas("cI0I1","I0 and I1",0,0,900,700); cI0I1->Divide(3,2); @@ -1577,9 +1787,13 @@ void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut, h3->Draw(); cI0I1->cd(5); h4->Draw(); + cI0I1->cd(6); + gStyle->SetPalette(1,0); + hI0I1s->Draw("col,Z"); if(save) { TFile *f = new TFile(fname,"recreate"); + hI0I1s->Write(); hI0->Write(); hI1->Write(); h2->Write(); @@ -1591,10 +1805,10 @@ void AliFastGlauber::PlotI0I1Distr(Int_t n,Double_t ellCut, } void AliFastGlauber::PlotI0I1B2BDistr(Int_t n,Double_t ellCut, - Bool_t save,Char_t *fname) + Bool_t save,const char *fname) { // - // Plot lengths back-to-back distributions + // Plot I0-I1 back-to-back distributions // Double_t i01,i11,i02,i12; TH2F *hI0s = new TH2F("hI0s","I_{0}'s back-to-back",100,0,100,100,0,100); @@ -1629,3 +1843,18 @@ void AliFastGlauber::PlotI0I1B2BDistr(Int_t n,Double_t ellCut, return; } +AliFastGlauber& AliFastGlauber::operator=(const AliFastGlauber& rhs) +{ +// Assignment operator + rhs.Copy(*this); + return *this; +} + +void AliFastGlauber::Copy(TObject&) const +{ + // + // Copy + // + Fatal("Copy","Not implemented!\n"); +} +