Updates
[u/mrichter/AliRoot.git] / EVGEN / AliGenAfterBurnerFlow.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 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 // 
20 // AliGenAfterBurnerFlow is a After Burner event generator applying flow.
21 // The generator changes Phi coordinate of the particle momentum.
22 // Flow (directed and elliptical) can be defined on particle type level
23 //
24 // Author:
25 // Sylwester Radomski, 2002
26 // Martin Poghosyan, 2008
27 // Constantin Loizides, 2010
28 //////////////////////////////////////////////////////////////////////////////
29
30 #include <Riostream.h>
31 #include <TParticle.h>
32 #include <TLorentzVector.h>
33 #include <TList.h>
34 #include <TRandom.h>
35 #include "AliStack.h"
36 #include "AliGenAfterBurnerFlow.h"
37 #include "AliGenCocktailAfterBurner.h"
38 #include "AliMC.h"
39 #include "AliRun.h"
40 #include "AliCollisionGeometry.h"
41 #include "AliGenCocktailEntry.h"
42
43
44 using std::cout;
45 using std::endl;
46 ClassImp(AliGenAfterBurnerFlow)
47
48
49 AliGenAfterBurnerFlow::AliGenAfterBurnerFlow():AliGenerator(),
50   fReactionPlane(0),
51   fHow(0),
52   fCounter(0),
53   fStack(0)
54 {
55   //
56   // Default Construction
57   InitPrimaries();
58   SetNpParams();
59 }
60
61 AliGenAfterBurnerFlow::AliGenAfterBurnerFlow(Float_t reactionPlane):AliGenerator(),
62   fReactionPlane(TMath::Pi()*reactionPlane/180.),
63   fHow(1),
64   fCounter(0),
65   fStack(0)
66 {
67   // reactionPlane - Reaction Plane Angle given in Deg [0-360]
68   // but stored and applied in radiants (standard for TParticle & AliCollisionGeometry) 
69
70   InitPrimaries();
71   SetNpParams();
72 }
73
74 ////////////////////////////////////////////////////////////////////////////////////////////////////
75
76 AliGenAfterBurnerFlow::~AliGenAfterBurnerFlow() 
77 {
78   // def. dest.
79 }
80
81 void AliGenAfterBurnerFlow::SetDirectedSimple(Int_t pdg, Float_t v1) 
82 {
83   //
84   // Set Directed Flow 
85   // The same directed flow is applied to all specified particles 
86   // independently on transverse momentum or rapidity
87   //
88   // PDG - particle type to apply directed flow
89   //       if (PDG == 0) use as default  
90   //
91
92   SetFlowParameters(pdg, 1, 0, v1, 0, 0, 0);
93 }
94
95 ////////////////////////////////////////////////////////////////////////////////////////////////////
96
97 void AliGenAfterBurnerFlow::SetDirectedParam(Int_t pdg, Float_t v11, Float_t v12, 
98                                              Float_t v13, Float_t v14) 
99 {
100   //
101   // Set Directed Flow 
102   // Directed flow is parameterised as follows
103   //
104   // V1(Pt,Y) = (V11 + V12*Pt) * sign(Y) * (V13 + V14 * abs(Y)^3)
105   //
106   // where sign = 1 for Y > 0 and -1 for Y < 0
107   // 
108   // Defaults values
109   // v12 = v14 = 0
110   // v13 = 1
111   //  
112   // PDG - particle type to apply directed flow
113   //       if (PDG == 0) use as default  
114   //
115   
116   SetFlowParameters(pdg, 1, 1, v11, v12, v13, v14);
117 }
118
119 ////////////////////////////////////////////////////////////////////////////////////////////////////
120
121 void AliGenAfterBurnerFlow::SetEllipticSimple(Int_t pdg, Float_t v2) 
122 {
123   //
124   // Set Elliptic Flow
125   // The same Elliptic flow is applied to all specified particles
126   // independently on transverse momentum or rapidity
127   //
128   // PDG - particle type to apply directed flow
129   //       if (PDG == 0) use as default  
130   //
131   // V2 - flow coefficient
132   //      
133   // NOTE: for starting playing with FLOW 
134   //       start with this function and values 0.05 - 0.1
135   //
136
137   SetFlowParameters(pdg, 2, 0, v2, 0, 0, 0);
138 }
139
140 ////////////////////////////////////////////////////////////////////////////////////////////////////
141
142 void AliGenAfterBurnerFlow::SetEllipticParam(Int_t pdg, 
143                                              Float_t v00, Float_t v10, Float_t v11,
144                                              Float_t v22) 
145 {
146   //
147   // Set Elliptic Flow
148   //
149   // Elliptic flow is parametrised to reproduce 
150   // V2 of Pions at RHIC energies and is given by:
151   // 
152   // V2 = (v00 + v10*pt + v11*pt^2) * exp (-v22 * y^2) and zero if V2<0.
153   //
154
155   SetFlowParameters(pdg, 2, 3, v00, v10, v11, v22);
156 }
157
158 void AliGenAfterBurnerFlow::SetEllipticParamPion(Int_t pdg, Float_t v21, 
159                                                  Float_t pTmax, Float_t v22) 
160 {
161   //
162   // Set Elliptic Flow
163   //
164   // Elliptic flow is parametrised to reproduce 
165   // V2 of Pions at RHIC energies and is given by:
166   // 
167   // V2 = v21 * (pT/pTMax ) * exp (-v22 * y^2)    where pT <= pTmax  
168   //      v21 * exp (-v22 * y^2)                  where pT > pTmax  
169   //
170   // v21   - value at saturation
171   // pTmax - saturation transverse momentum
172   // v22   - rapidity decreasing
173   //
174
175   SetFlowParameters(pdg, 2, 1, v21, pTmax, v22, 0);
176 }
177
178 ////////////////////////////////////////////////////////////////////////////////////////////////////
179
180 void AliGenAfterBurnerFlow::SetEllipticParamOld(Int_t pdg, Float_t v21, Float_t v22, Float_t v23) 
181 {
182   //
183   // Set Elliptic Flow
184   //
185   // Elliptic flow is parameterised using 
186   // old MevSim parameterisation
187   // 
188   // V2 = (V21 + V22 pT^2) * exp (-v22 * y^2)
189   //
190
191   SetFlowParameters(pdg, 2, 2, v21, v22, v23, 0);
192 }
193
194 ////////////////////////////////////////////////////////////////////////////////////////////////////
195
196 void AliGenAfterBurnerFlow::SetNpParams(Int_t order, Float_t p0, Float_t p1, Float_t p2, Float_t p3)
197 {
198   //
199   // Set npart parameterization.
200   //
201
202   fNpParams[0] = order;
203   fNpParams[1] = p0;
204   fNpParams[2] = p1;
205   fNpParams[3] = p2;
206   fNpParams[4] = p3;
207 }
208
209 ////////////////////////////////////////////////////////////////////////////////////////////////////
210
211 void AliGenAfterBurnerFlow::SetFlowParameters(Int_t pdg, Int_t order, Int_t type, 
212                                               Float_t v1, Float_t v2,Float_t v3,Float_t v4) 
213 {
214   // 
215   // private function
216   // 
217
218   if(TMath::Abs(pdg)>=fgkPDG){
219     Error("AliAfterBurnerFlow","Overflow");
220     return;
221   }
222   fIsPrim[TMath::Abs(pdg)]=kTRUE;
223   
224   Int_t index = 0;
225   Bool_t newEntry = kTRUE;
226
227   // Defaults
228   if (pdg == 0) {
229     index = fgkN - order;
230     newEntry = kFALSE;
231   }
232
233   // try to find existing entry
234   for (Int_t i=0; i<fCounter; i++) {
235     if (pdg == (Int_t)fParams[i][0] && 
236         order == (Int_t)fParams[i][1]) {
237       
238       index = i;
239       newEntry = kFALSE;
240     }
241   }
242   
243   // check fCounter
244
245   if (newEntry && (fCounter > fgkN-3)) {
246     Error("AliAfterBurnerFlow","Overflow");
247     return;
248   }
249   
250   if (newEntry) {
251     index = fCounter;
252     fCounter++;
253   }
254   
255   // Set new particle type
256
257   fParams[index][0] = pdg;
258   fParams[index][1] = order;
259   fParams[index][2] = type;
260   fParams[index][3] = v1;
261   fParams[index][4] = v2;
262   fParams[index][5] = v3;
263   fParams[index][6] = v4;  
264 }
265
266 ////////////////////////////////////////////////////////////////////////////////////////////////////
267
268 void AliGenAfterBurnerFlow::Init() 
269 {
270   // 
271   // Standard AliGenerator Initializer
272   //
273
274   if(fHow == 0)     { Info("AliGenAfterBurnerFlow", "Using the Hijing R.P. Angle event by event "); }
275   else if(fHow == 1){ Info("AliGenAfterBurnerFlow", "Using a fixed R.P. Angle for every event ") ; }
276   else { Info("AliGenAfterBurnerFlow", 
277               "Using a random R.P. Angle event by event ( ! not the same used by Hijing ! ) "); }
278 }
279
280 ////////////////////////////////////////////////////////////////////////////////////////////////////
281
282 Float_t AliGenAfterBurnerFlow::GetCoefficient(Int_t pdg, Int_t n, Float_t Pt, Float_t Y) const
283 {
284   //
285   // private function
286   // Return Flow Coefficient for a given particle type flow order
287   // and particle momentum (Pt, Y)
288   //
289
290   Int_t index = fgkN - n;  // default index (for all pdg)
291   Float_t v = 0;
292
293   // try to find specific parametrs
294
295   for (Int_t i=0; i<fCounter; i++) {
296     if ((Int_t)fParams[i][0] == pdg &&
297         (Int_t)fParams[i][1] == n) {
298       
299       index = i;
300       break;
301     }
302   } 
303   
304   // calculate v
305   
306   Int_t type = (Int_t)fParams[index][2];
307
308   if ((Int_t)fParams[index][1] == 1) { // Directed
309
310     if (type == 0 )
311       v = fParams[index][3];
312     else 
313       v = (fParams[index][3] + fParams[index][4] * Pt) * TMath::Sign((Float_t)1.,Y) *
314         (fParams[index][5] + fParams[index][6] * TMath::Abs(Y*Y*Y) );
315
316   } else {  // Elliptic
317
318     if (type == 0) v = fParams[index][3];
319
320     // Pion parameterisation 
321     if (type == 1) { 
322       if (Pt < fParams[index][4]) 
323         v = fParams[index][3] * (Pt / fParams[index][4]) ;
324       else 
325         v = fParams[index][3];
326       
327       v *= TMath::Exp( - fParams[index][5] * Y * Y);
328     }
329
330     // Old parameterisation
331     if (type == 2) 
332       v = (fParams[index][3] + fParams[index][4] * Pt * Pt) *
333         TMath::Exp( - fParams[index][5] * Y * Y);
334
335     // New v2 parameterisation
336     if (type == 3) {
337       v = (fParams[index][3] + fParams[index][4] *Pt + fParams[index][5] *Pt*Pt) *
338         TMath::Exp( - fParams[index][6] * Y*Y);
339       if (v<0)
340         v = 0;
341     }
342   }
343   
344   return v;
345 }
346
347 ////////////////////////////////////////////////////////////////////////////////////////////////////
348
349 Float_t AliGenAfterBurnerFlow::GetNpNorm(Int_t npart) const
350 {
351   //
352   // Calculate npart norm.
353   //
354
355   if (npart<0)
356     return 1;
357
358   Int_t order = (Int_t)fNpParams[0];
359   if (order<0)
360     return 1;
361
362   Float_t ret = 0;
363   Int_t npp = 1;
364   for (Int_t i=0; i<=order; i++) {
365     ret += npp*fNpParams[i+1];
366     npp *= npart;
367   } 
368   return ret;
369 }
370
371 ////////////////////////////////////////////////////////////////////////////////////////////////////
372
373 Bool_t AliGenAfterBurnerFlow::IsPrimary(Int_t pdg) const
374 {
375   if(pdg>=fgkPDG) return kFALSE;
376   return fIsPrim[pdg];
377 }
378
379 ////////////////////////////////////////////////////////////////////////////////////////////////////
380
381 Double_t CalcAngle(Double_t phi, Double_t phi0, Double_t phiRP, Double_t v2, Double_t v1=0.)
382 {
383     // Calculate relative angle
384   Double_t phi1 = phi-(phi+2*v1*TMath::Sin(phi-phiRP)+v2*TMath::Sin(2*(phi-phiRP))-phi0)/
385     (1.+2*v1*TMath::Cos(phi-phiRP)+ 2*v2*TMath::Cos(2*(phi-phiRP)));
386   if(TMath::Abs(phi/phi1-1.)<0.00001) return phi1;
387   return CalcAngle(phi1, phi0, phiRP, v2, v1);
388 }
389
390 ////////////////////////////////////////////////////////////////////////////////////////////////////
391
392 void AliGenAfterBurnerFlow::InitPrimaries()
393 {
394   // Init the primary particle list
395   for(Int_t i=0; i<fgkPDG; i++) fIsPrim[i]=kFALSE;
396
397   //mesons
398   fIsPrim[211]=kTRUE;
399   fIsPrim[311]=kTRUE;
400   fIsPrim[321]=kTRUE;
401   fIsPrim[411]=kTRUE;
402   fIsPrim[421]=kTRUE;
403   fIsPrim[431]=kTRUE;
404   fIsPrim[511]=kTRUE;
405   fIsPrim[521]=kTRUE;
406   fIsPrim[531]=kTRUE;
407   fIsPrim[541]=kTRUE;
408   fIsPrim[111]=kTRUE;
409   fIsPrim[221]=kTRUE;
410   fIsPrim[331]=kTRUE;
411   fIsPrim[441]=kTRUE;
412   fIsPrim[551]=kTRUE;
413   fIsPrim[130]=kTRUE;
414   fIsPrim[310]=kTRUE;
415   fIsPrim[213]=kTRUE;
416   fIsPrim[313]=kTRUE;
417   fIsPrim[323]=kTRUE;
418   fIsPrim[413]=kTRUE;
419   fIsPrim[423]=kTRUE;
420   fIsPrim[433]=kTRUE;
421   fIsPrim[513]=kTRUE;
422   fIsPrim[523]=kTRUE;
423   fIsPrim[533]=kTRUE;
424   fIsPrim[543]=kTRUE;
425   fIsPrim[113]=kTRUE;
426   fIsPrim[223]=kTRUE;
427   fIsPrim[333]=kTRUE;
428   fIsPrim[443]=kTRUE;
429   fIsPrim[553]=kTRUE;
430
431   //baryons
432   fIsPrim[2112]=kTRUE;
433   fIsPrim[2212]=kTRUE;
434   fIsPrim[3112]=kTRUE;
435   fIsPrim[3122]=kTRUE;
436   fIsPrim[3212]=kTRUE;
437   fIsPrim[3222]=kTRUE;
438   fIsPrim[3312]=kTRUE;
439   fIsPrim[3322]=kTRUE;
440   fIsPrim[4112]=kTRUE;
441   fIsPrim[4122]=kTRUE;
442   fIsPrim[4212]=kTRUE;
443   fIsPrim[4222]=kTRUE;
444   fIsPrim[4132]=kTRUE;
445   fIsPrim[4312]=kTRUE;
446   fIsPrim[4232]=kTRUE;
447   fIsPrim[4322]=kTRUE;
448   fIsPrim[4332]=kTRUE;
449   fIsPrim[5112]=kTRUE;
450   fIsPrim[5122]=kTRUE;
451   fIsPrim[5212]=kTRUE;
452   fIsPrim[5222]=kTRUE;
453   fIsPrim[1114]=kTRUE;
454   fIsPrim[2114]=kTRUE;
455   fIsPrim[2214]=kTRUE;
456   fIsPrim[2224]=kTRUE;
457   fIsPrim[3114]=kTRUE;
458   fIsPrim[3214]=kTRUE;
459   fIsPrim[3224]=kTRUE;
460   fIsPrim[3314]=kTRUE;
461   fIsPrim[3324]=kTRUE;
462   fIsPrim[3334]=kTRUE;
463   fIsPrim[4114]=kTRUE;
464   fIsPrim[4214]=kTRUE;
465   fIsPrim[4224]=kTRUE;
466   fIsPrim[4314]=kTRUE;
467   fIsPrim[4324]=kTRUE;
468   fIsPrim[4334]=kTRUE;
469   fIsPrim[5114]=kTRUE;
470   fIsPrim[5214]=kTRUE;
471   fIsPrim[5224]=kTRUE;
472 }
473
474 ////////////////////////////////////////////////////////////////////////////////////////////////////
475
476 void AliGenAfterBurnerFlow::Generate() 
477 {
478   // 
479   // AliGenerator generate function doing actual job.
480   // Algorythm:
481   //
482   // 1. loop over particles on the stack and choose primaries
483   // 2. calculate delta phi
484   // 3. change phi of primary particle and if it is non-stable 
485   //    then its daughters' phi and vertex also
486   // 
487   // For more details see :
488   // M.G. Poghosyan 
489   // PWG2 meeting on 06.05.2008 and 03.06.2008
490
491
492   if (0) 
493     for(Int_t ii=0; ii<fCounter;ii++)
494     {
495       printf("%d  %f  %f  %f  %f\n",ii,fParams[ii][0],fParams[ii][1],fParams[ii][2],fParams[ii][3]);
496     } 
497
498   AliGenCocktailAfterBurner *gen;
499   
500   TParticle *particle;
501   TParticle *particleM;
502   TLorentzVector momentum;
503   TLorentzVector vertex;
504
505   Int_t pdg;
506   Float_t phi;
507   Float_t pt, y;
508
509   // Get Stack of the first Generator
510   //  gen = (AliGenCocktailAfterBurner *)gAlice->Generator();
511   gen = (AliGenCocktailAfterBurner *)gAlice->GetMCApp()->Generator();
512
513
514   AliGenerator* genHijing = 0 ;
515   AliCollisionGeometry* geom = 0 ;
516   AliGenCocktailEntry* entry = 0 ;
517   TList* fEntries = 0 ;
518
519   TRandom* rand = new TRandom(0) ;
520   for(Int_t ns=0;ns<gen->GetNumberOfEvents();ns++) 
521   {
522     gen->SetActiveEventNumber(ns) ;
523    
524     fStack = gen->GetStack(ns);
525     fEntries = gen->Entries() ;
526
527     TIter next(fEntries) ;
528     Int_t npart = -1;
529
530     if(fHow == 0) // hijing R.P.
531     {
532       while((entry = (AliGenCocktailEntry*)next())) 
533       {
534         Info("Generate (e)","Using R.P. from HIJING ... ");       
535         genHijing = entry->Generator() ;        
536         if(genHijing->ProvidesCollisionGeometry()) 
537         { 
538           geom = gen->GetCollisionGeometry(ns) ;
539           fReactionPlane = geom->ReactionPlaneAngle() ; 
540           npart =  geom->ProjectileParticipants() + geom->TargetParticipants();
541           break;
542         }
543         else 
544         {
545           Error("Generate (e)", "NO CollisionGeometry !!!  -  using fixed R.P. angle = 0. ") ; 
546           fReactionPlane = 0. ; 
547         }
548       }
549     }
550     else if(fHow ==1 ) //  fixed R.P. 
551     {
552       Info("Generate (e)","Using fixed R.P. ...");       
553     }
554     else
555     {
556       Info("Generate (e)","Using random R.P.s ... ");       
557       fReactionPlane = 2 * TMath::Pi() * rand->Rndm() ;
558     }    
559     
560     cout << " * Reaction Plane Angle (event " << ns << ") = " << fReactionPlane << 
561       " rad. ( = " << (360*fReactionPlane/(2*TMath::Pi())) << " deg.) Npart = " << npart << "* " << endl ;
562
563     Int_t nParticles = fStack->GetNprimary();
564     for (Int_t i=0; i<nParticles; i++) 
565     {
566       particle = fStack->Particle(i);
567  
568       Int_t iM=particle->GetMother(0);
569       pdg = particle->GetPdgCode();
570
571       //exclude incoming protons in PYTHIA     
572       if(particle->GetPdgCode()==21) continue;
573
574       if(TMath::Abs(pdg)>=fgkPDG) continue; 
575       // is particle primary?
576       if(!fIsPrim[TMath::Abs(pdg)]) continue; 
577
578       if(iM>0)
579       {
580         particleM = fStack->Particle(iM);
581         Int_t pdgM = TMath::Abs(particleM->GetPdgCode());
582         // is mother primary?
583         if((TMath::Abs(pdgM)<fgkPDG)&&fIsPrim[TMath::Abs(pdgM)]) continue;
584       }
585    
586       particle->Momentum(momentum);
587       phi = particle->Phi();
588
589       // get Pt, y    
590       pt = momentum.Pt() ; 
591       y = 10000.;
592
593       if(TMath::Abs(momentum.Z()) != TMath::Abs(momentum.T())) 
594         y = momentum.Rapidity() ; 
595       
596       Double_t v1 = GetCoefficient(pdg, 1, pt, y);
597       Double_t v2 = GetCoefficient(pdg, 2, pt, y);
598       Double_t npartnorm = GetNpNorm(npart);
599       v2 *= npartnorm;
600
601       //printf("ntup %d %f  %f %f %f %f\n ",npart, v1, v2, pt, y, npartnorm);
602      
603       Double_t phi1 = CalcAngle(phi, phi, fReactionPlane,v2,v1);
604      
605       Rotate(i, phi1-phi);
606     }
607   }
608
609   Info("Generate","Flow After Burner: DONE");  
610 }
611
612 ////////////////////////////////////////////////////////////////////////////////////////////////////
613
614 void AliGenAfterBurnerFlow::Rotate(Int_t i, Double_t phi, Bool_t IsPrim)
615 {
616   // Rotation
617   TParticle*  particle = fStack->Particle(i);
618   
619   TLorentzVector momentum;
620   particle->Momentum(momentum);
621   momentum.RotateZ(phi);
622   particle->SetMomentum(momentum);
623  
624   if(!IsPrim)
625   {
626     TLorentzVector vertex;
627     particle->ProductionVertex(vertex);
628     vertex.RotateZ(phi);
629     particle->SetProductionVertex(vertex);
630   }
631   
632   if(particle->GetFirstDaughter()<0) return;     
633   for(Int_t iD=particle->GetFirstDaughter(); iD<=particle->GetLastDaughter(); iD++) Rotate(iD, phi, kFALSE);
634
635   return;
636 }
637
638