Line data Source code
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 : /* $Id: AliGenTunedOnPbPb.cxx 51126 2013-08-19 13:37:49Z fnoferin $ */
17 :
18 : // Parameterisation based on 5.5 ATeV PbPb data
19 : // pi, K, p, neutron, K0, lambda, phi, Xi, Omega spectra, v2, v3 (no jets!)
20 : // Author: fnoferin@cern.ch
21 :
22 : #include <TArrayF.h>
23 : #include <TCanvas.h>
24 : #include <TClonesArray.h>
25 : #include <TDatabasePDG.h>
26 : #include <TF1.h>
27 : #include <TH1.h>
28 : #include <TPDGCode.h>
29 : #include <TParticle.h>
30 : #include <TROOT.h>
31 : #include <TVirtualMC.h>
32 : #include <TLorentzVector.h>
33 :
34 : #include "AliConst.h"
35 : #include "AliDecayer.h"
36 : #include "AliGenEventHeaderTunedPbPb.h"
37 : #include "AliLog.h"
38 : #include "AliRun.h"
39 : #include "AliGenTunedOnPbPb.h"
40 :
41 6 : ClassImp(AliGenTunedOnPbPb)
42 :
43 : // set default parameters for 10-20% centrality
44 : Int_t AliGenTunedOnPbPb::fgPdgInput[fgNspecies] = {211,-211,111,321,-321,2212,-2212,310,3122,-3122,333,3312,-3312,3334,-3334,2112,-2112};
45 : Float_t AliGenTunedOnPbPb::fgMult[fgNspecies] = {450,450,450,70,70,21,21,70,20,20,8,2.4,2.4,0.4,0.4,21,21};
46 :
47 : Float_t AliGenTunedOnPbPb::fgV3Overv2 = 6.25000000000000000e-01;
48 : Float_t AliGenTunedOnPbPb::fgEventplane=0;
49 :
50 : TF1 *AliGenTunedOnPbPb::fgV2 = NULL;
51 :
52 : //_____________________________________________________________________________
53 : AliGenTunedOnPbPb::AliGenTunedOnPbPb()
54 0 : :AliGenerator(),
55 0 : fCmin(0.),
56 0 : fCmax(100.),
57 0 : fChangeWithCentrality(kFALSE),
58 0 : fYMaxValue(2.0),
59 0 : fYlimitForFlatness(2.0),
60 0 : fYdecreaseSp(0.2),
61 0 : fYdecreaseV2(0.2)
62 0 : {
63 : //
64 : // Default constructor
65 : //
66 0 : SetCutVertexZ();
67 0 : SetPtRange();
68 :
69 0 : for(Int_t i=0;i < fgNspecies;i++){
70 0 : fgHSpectrum[i] = NULL;
71 0 : fgHv2[i] = NULL;
72 : }
73 0 : }
74 :
75 : //_____________________________________________________________________________
76 0 : AliGenTunedOnPbPb::~AliGenTunedOnPbPb()
77 0 : {
78 : //
79 : // Standard destructor
80 : //
81 0 : }
82 :
83 : //_____________________________________________________________________________
84 : void AliGenTunedOnPbPb::Init()
85 : {
86 : //
87 : // Initialise the generator
88 : //
89 :
90 : // define histos
91 0 : }
92 :
93 :
94 : //_____________________________________________________________________________
95 : void AliGenTunedOnPbPb::Generate()
96 : {
97 : //
98 : // Generate one trigger
99 : //
100 :
101 0 : Float_t avCentr = (fCmin+fCmax)*0.5;
102 :
103 : Float_t centrality = avCentr;
104 :
105 0 : if(fChangeWithCentrality) centrality = fCmin + gRandom->Rndm()*(fCmax-fCmin);
106 :
107 0 : SetParameters(centrality);
108 :
109 0 : if(!fChangeWithCentrality){
110 : Float_t in=0;
111 0 : for(Int_t i=0;i < fgNspecies;i++){
112 : in=0;
113 0 : if(fgHSpectrum[i]){
114 0 : for(Int_t j=1;j<=fgHSpectrum[i]->GetNbinsX();j++){
115 0 : in += fgHSpectrum[i]->GetBinContent(j)*fgHSpectrum[i]->GetBinWidth(j);
116 : }
117 0 : }
118 :
119 : // replace n-particles with the one in input file if centralidy dependece was disable
120 0 : fgMult[i] = in;
121 : }
122 0 : }
123 :
124 :
125 0 : TMCProcess statusPdg[fgNspecies] = {kPPrimary,kPPrimary,kPPrimary,kPPrimary,kPPrimary,kPPrimary,kPPrimary,kPPrimary,kPPrimary,kPPrimary,kPPrimary,kPPrimary,kPPrimary,kPPrimary,kPPrimary,kPPrimary,kPPrimary};
126 :
127 : Float_t parV2scaling[3] = {1,0.202538,-0.00214468};
128 :
129 : Float_t scaleV2 = 1.0;
130 :
131 0 : TDatabasePDG *pdgD = TDatabasePDG::Instance();
132 :
133 0 : if(fChangeWithCentrality){
134 0 : parV2scaling[0] = 1. / (1 + parV2scaling[1]*avCentr + parV2scaling[2]*avCentr*avCentr); // normalize to average centrality
135 :
136 0 : scaleV2 = parV2scaling[0]*(1 + parV2scaling[1]*centrality + parV2scaling[2]*centrality*centrality); // apply the trand of v2 w.r.t. centrality
137 0 : }
138 :
139 0 : fgV2->SetParameter(2,fgV3Overv2);
140 :
141 0 : Float_t psi = gRandom->Rndm()*TMath::Pi();
142 0 : fgEventplane = psi;
143 0 : Float_t psi3 = gRandom->Rndm()*TMath::Pi()*2/3;
144 0 : Float_t psi4 = gRandom->Rndm()*TMath::Pi()*2/4;
145 0 : fgV2->SetParameter(1,psi);
146 0 : fgV2->SetParameter(3,psi3);
147 0 : fgV2->SetParameter(4,psi4);
148 :
149 0 : Int_t npart = 0;
150 :
151 0 : Float_t origin0[3];
152 0 : for (Int_t j=0;j<3;j++) origin0[j]=fOrigin[j];
153 : Float_t time;
154 0 : time = fTimeOrigin;
155 0 : if(fVertexSmear == kPerEvent) {
156 0 : Vertex();
157 0 : for (Int_t j=0; j < 3; j++) origin0[j] = fVertex[j];
158 0 : time = fTime;
159 0 : } // if kPerEvent
160 :
161 0 : printf("Generate event with centrality = %3.1f%c, |y|<%4.1f\n",centrality,'%',fYMaxValue);
162 :
163 0 : for(Int_t isp=0;isp < fgNspecies;isp++){
164 0 : if(! fgHSpectrum[isp]) continue;
165 :
166 0 : Int_t npartSp = Int_t(fgMult[isp]*2*fYMaxValue + gRandom->Rndm());
167 :
168 0 : printf("Total number of %i = %i\n",fgPdgInput[isp],npartSp);
169 :
170 0 : for(Int_t ipart =0; ipart < npartSp; ipart++){
171 0 : Int_t pdg = fgPdgInput[isp];
172 :
173 0 : Double_t y = gRandom->Rndm()*2*fYMaxValue - fYMaxValue;
174 0 : Double_t ytanh = TMath::TanH(y);
175 0 : Double_t pt = fgHSpectrum[isp]->GetRandom();
176 0 : Double_t mass = pdgD->GetParticle(pdg)->Mass();
177 0 : Double_t mt2 = pt*pt + mass*mass;
178 0 : Double_t pz = ytanh*TMath::Sqrt(mt2)/TMath::Sqrt(1-ytanh*ytanh);
179 0 : Double_t etot = TMath::Sqrt(mt2 + pz*pz);
180 0 : TLorentzVector tempVect(pt,0,pz,etot);
181 : // Double_t eta = tempVect.PseudoRapidity();
182 : Double_t scaleEtaV2 = 1; // set the eta dependence
183 0 : if(TMath::Abs(y)> fYlimitForFlatness) scaleEtaV2 = 1 - fYdecreaseV2*(TMath::Abs(y) - fYlimitForFlatness);
184 :
185 0 : if(fgHv2[isp]) fgV2->SetParameter(0,fgHv2[isp]->Interpolate(pt) * scaleV2 * scaleEtaV2);
186 0 : else fgV2->SetParameter(0,0.);
187 0 : Double_t phi = fgV2->GetRandom(-TMath::Pi(),TMath::Pi());
188 0 : Double_t px = pt*TMath::Cos(phi);
189 0 : Double_t py = pt*TMath::Sin(phi);
190 0 : Float_t p[3] = {static_cast<Float_t>(px),static_cast<Float_t>(py),static_cast<Float_t>(pz)};
191 0 : Float_t polar[3] = {0.,0.,0.};
192 :
193 0 : if(TMath::Abs(y)< fYlimitForFlatness || gRandom->Rndm() < 1 - fYdecreaseSp*(TMath::Abs(y) - fYlimitForFlatness)){// check on pseudorapidity distribution
194 : // printf("%f %f\n",eta,phi - psi); // for debugging
195 0 : PushTrack(1, -1, pdg, p, origin0, polar, time, statusPdg[isp], npart, 1., 1);
196 0 : KeepTrack(npart);
197 0 : npart++;
198 0 : }
199 0 : }
200 0 : }
201 :
202 0 : TArrayF eventVertex;
203 0 : eventVertex.Set(3);
204 0 : eventVertex[0] = origin0[0];
205 0 : eventVertex[1] = origin0[1];
206 0 : eventVertex[2] = origin0[2];
207 :
208 : // Header
209 0 : AliGenEventHeaderTunedPbPb* header = new AliGenEventHeaderTunedPbPb("tunedOnPbPb");
210 : // Event Vertex
211 0 : header->SetPrimaryVertex(eventVertex);
212 0 : header->SetInteractionTime(time);
213 0 : header->SetNProduced(npart);
214 0 : header->SetCentrality(centrality);
215 0 : header->SetPsi2(psi);
216 0 : header->SetPsi3(psi3);
217 0 : header->SetPsi4(psi4);
218 0 : gAlice->SetGenEventHeader(header);
219 0 : }
220 :
221 : void AliGenTunedOnPbPb::SetPtRange(Float_t ptmin, Float_t ptmax) {
222 0 : AliGenerator::SetPtRange(ptmin, ptmax);
223 0 : }
224 : //_____________________________________________________________________________
225 : TH1F *AliGenTunedOnPbPb::GetMultVsCentrality(Int_t species){
226 0 : char title[100];
227 0 : snprintf(title,100,"pdg = %i;centrality;dN/dy",fgPdgInput[species]);
228 0 : TH1F *h = new TH1F("multVsCentr",title,100,0,100);
229 :
230 0 : for(Int_t i=1;i<=100;i++){
231 0 : Float_t x = i+0.5;
232 0 : SetParameters(x);
233 0 : h->SetBinContent(i,fgMult[species]);
234 : }
235 :
236 0 : return h;
237 0 : }
238 : //_____________________________________________________________________________
239 : void AliGenTunedOnPbPb::SetParameters(Float_t centrality){
240 :
241 0 : if(!fgV2) fgV2 = new TF1("fv2Par","TMath::Max(0.,(1 + 2*[0]*cos(2*(x-[1])) + 2*[0]*[2]*cos(3*(x-[3])) + [0]*[2]*cos(4*(x-[4]))))",-TMath::Pi(),TMath::Pi()); // v4 is approx. 0.5*v3
242 :
243 0 : Float_t fr[9] = {0.,0.,0.,0.,0.,0.,0.,0.,0.};
244 :
245 0 : if(centrality < 7.5){
246 0 : fr[0] = (7.5 - centrality)/5.;
247 0 : fr[1] = (centrality-2.5)/5.;
248 0 : }
249 0 : else if(centrality < 15){
250 0 : fr[1] = (15-centrality)/7.5;
251 0 : fr[2] = (centrality-7.5)/7.5;
252 0 : }
253 0 : else if(centrality < 25){
254 0 : fr[2] = (25-centrality)/10.;
255 0 : fr[3] = (centrality-15)/10.;
256 0 : }
257 0 : else if(centrality < 35){
258 0 : fr[3] = (35-centrality)/10.;
259 0 : fr[4] = (centrality-25)/10.;
260 0 : }
261 0 : else if(centrality < 45){
262 0 : fr[4] = (45-centrality)/10.;
263 0 : fr[5] = (centrality-35)/10.;
264 0 : }
265 0 : else if(centrality < 55){
266 0 : fr[5] = (55-centrality)/10.;
267 0 : fr[6] = (centrality-45)/10.;
268 0 : }
269 0 : else if(centrality < 65){
270 0 : fr[6] = (65-centrality)/10.;
271 0 : fr[7] = (centrality-55)/10.;
272 0 : }
273 0 : else if(centrality < 75){
274 0 : fr[7] = (75-centrality)/10.;
275 0 : fr[8] = (centrality-65)/10.;
276 0 : }
277 : else{
278 0 : fr[8] = 1.0;
279 : }
280 :
281 : // parameters as a function of centrality
282 0 : Float_t multCent[9*fgNspecies] = {
283 : 733.,733.,733.,109.,109.,34.0,34.0,109.,28.,28.,11.5,3.1 ,3.1 ,0.5 ,0.5 ,34.0,34.0,
284 : 606.,606.,606.,91.0,91.0,28.0,28.0,91. ,24.,24.,9. ,2.7 ,2.7 ,0.45 ,0.45 ,28.0,28.0,
285 : 455.,455.,455.,68.0,68.0,21.0,21.0,68. ,20.,20.,8. ,2.4 ,2.4 ,0.40 ,0.40 ,21.0,21.0,
286 : 307.,307.,307.,46.0,46.0,14.5,14.5,46. ,14.,14.,5.5 ,1.5 ,1.5 ,0.2 ,0.2 ,14.5,14.5,
287 : 201.,201.,201.,30.0,30.0,9.60,9.60,30. ,9. ,9. ,3.5 ,0.9 ,0.9 ,0.08 ,0.08 ,9.60,9.60,
288 : 124.,124.,124.,18.3,18.3,6.10,6.10,18.3,5.1,5.1,2.2 ,0.6 ,0.6 ,0.055,0.055,6.10,6.10,
289 : 71.0,71.0,71.0,10.2,10.2,3.60,3.60,10.2,2.6,2.6,1.4 ,0.36 ,0.36 ,0.035,0.035,3.60,3.60,
290 : 37.0,37.0,37.0,5.10,5.10,2.00,2.00,5.10,1.5,1.5,0.50,0.020,0.020,0.015,0.015,2.00,2.00,
291 : 17.0,17.0,17.0,2.30,2.30,0.90,0.90,2.30,0.6,0.6,0.16,0.006,0.006,0.005,0.005,0.90,0.90
292 : };
293 :
294 0 : Float_t v3Overv2Cent[9] = {1.2,0.82,0.625,0.5,0.45,0.4,0.37,0.3,0.3};
295 :
296 0 : fgV3Overv2 = 0;
297 0 : for(Int_t j=0;j < 9;j++)
298 0 : fgV3Overv2 += fr[j]*v3Overv2Cent[j];
299 :
300 : // set parameters for current centrality
301 0 : for(Int_t i=0;i < fgNspecies;i++){
302 0 : fgMult[i] = 0;
303 :
304 0 : for(Int_t j=0;j < 9;j++){
305 0 : fgMult[i] += fr[j]*multCent[i+j*fgNspecies];
306 : }
307 : }
308 :
309 0 : if(centrality > 80){
310 0 : for(Int_t i=0;i < fgNspecies;i++)
311 0 : fgMult[i] /= TMath::Log(centrality-77.);
312 0 : }
313 0 : }
|
