few extra options for momentum dependence v_n
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowEventSimpleMakerOnTheFly.cxx
1 /*************************************************************************
2 * Copyright(c) 1998-2008, 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 /********************************** 
17  * create an event and perform    *
18  * flow analysis 'on the fly'     * 
19  *                                * 
20  * authors: Raimond Snellings     *
21  *           (snelling@nikhef.nl) * 
22  *          Ante Bilandzic        * 
23  *           (anteb@nikhef.nl)    *
24  *********************************/
25   
26 #include "Riostream.h"
27 #include "TMath.h"
28 #include "TF1.h"
29 #include "TRandom3.h"
30 //#include "TUnuran.h"
31
32 #include "AliFlowEventSimpleMakerOnTheFly.h"
33 #include "AliFlowEventSimple.h"
34 #include "AliFlowTrackSimple.h"
35
36 ClassImp(AliFlowEventSimpleMakerOnTheFly)
37
38
39 //========================================================================
40
41
42 AliFlowEventSimpleMakerOnTheFly::AliFlowEventSimpleMakerOnTheFly(UInt_t iseed):
43   fMultiplicityOfRP(0),
44   fMultiplicitySpreadOfRP(0.),
45   fUseConstantHarmonics(kFALSE),
46   fV1RP(0.), 
47   fV1SpreadRP(0.), 
48   fV2RP(0.), 
49   fV2SpreadRP(0.), 
50   fV4RP(0.), 
51   fV4SpreadRP(0.), 
52   fV2RPMax(0.), 
53   fPtCutOff(0.), 
54   fPtSpectra(NULL),
55   fPhiDistribution(NULL),
56   fMyTRandom3(NULL),
57   fCount(0),
58   fNoOfLoops(1)
59  {
60   // constructor
61   fMyTRandom3 = new TRandom3(iseed);   
62   gRandom->SetSeed(fMyTRandom3->Integer(65539));
63  
64   //fMyUnuran = new TUnuran(); 
65   //fMyUnuran->SetSeed(iseed);  
66  }
67
68
69 //========================================================================
70
71
72 AliFlowEventSimpleMakerOnTheFly::~AliFlowEventSimpleMakerOnTheFly()
73 {
74  // destructor
75  if (fPtSpectra) delete fPtSpectra;
76  if (fPhiDistribution) delete fPhiDistribution;
77  if (fMyTRandom3) delete fMyTRandom3;
78 }
79
80
81 //========================================================================
82
83
84 void AliFlowEventSimpleMakerOnTheFly::Init()
85 {
86  // define the pt spectra and phi distribution
87  
88  // pt spectra:   
89  Double_t dPtMin = 0.; // to be improved (move this to the body of contstructor?)
90  Double_t dPtMax = 10.; // to be improved (move this to the body of contstructor?) 
91   
92  fPtSpectra = new TF1("fPtSpectra","[0]*x*TMath::Exp(-x*x)",dPtMin,dPtMax);  
93  fPtSpectra->SetParName(0,"Multiplicity of RPs");  
94  
95  // phi distribution:
96  Double_t dPhiMin = 0.; // to be improved (move this to the body of contstructor?)
97  Double_t dPhiMax = TMath::TwoPi(); // to be improved (move this to the body of contstructor?)
98   
99  fPhiDistribution = new TF1("fPhiDistribution","1+2.*[0]*TMath::Cos(x-[2])+2.*[1]*TMath::Cos(2*(x-[2]))+2.*[3]*TMath::Cos(4*(x-[2]))",dPhiMin,dPhiMax);
100  fPhiDistribution->SetParName(0,"directed flow");
101  fPhiDistribution->SetParName(1,"elliptic flow"); 
102  fPhiDistribution->SetParName(2,"Reaction Plane");
103  fPhiDistribution->SetParName(3,"harmonic 4"); // to be improved (name)
104 }
105
106 //========================================================================
107
108 AliFlowEventSimple* AliFlowEventSimpleMakerOnTheFly::CreateEventOnTheFly()
109 {
110   // method to create event on the fly
111   
112   AliFlowEventSimple* pEvent = new AliFlowEventSimple(fMultiplicityOfRP);
113    
114   // sampling the multiplicity:
115   Int_t iNewMultiplicityOfRP = fMultiplicityOfRP;
116   if(fMultiplicitySpreadOfRP>0.0) iNewMultiplicityOfRP = (Int_t)fMyTRandom3->Gaus(fMultiplicityOfRP,fMultiplicitySpreadOfRP);
117   fPtSpectra->SetParameter(0,iNewMultiplicityOfRP);
118   
119   // sampling the reaction plane
120   Double_t dMCReactionPlaneAngle = fMyTRandom3->Uniform(0.,TMath::TwoPi());
121   fPhiDistribution->SetParameter(2,dMCReactionPlaneAngle);
122
123   // sampling the V1:
124   Double_t dNewV1RP=fV1RP;
125   if(fV1SpreadRP>0.0) {dNewV1RP = fMyTRandom3->Gaus(fV1RP,fV1SpreadRP);}
126   fPhiDistribution->SetParameter(0,dNewV1RP);
127  
128   // sampling the V2:
129   if(fUseConstantHarmonics)
130   {
131    Double_t dNewV2RP = fV2RP;
132    if(fV2SpreadRP>0.0) dNewV2RP = fMyTRandom3->Gaus(fV2RP,fV2SpreadRP);
133    fPhiDistribution->SetParameter(1,dNewV2RP);
134   }
135   
136   // sampling the V4:
137   Double_t dNewV4RP = fV4RP;
138   if(fV4SpreadRP>0.0) dNewV4RP = fMyTRandom3->Gaus(fV4RP,fV4SpreadRP);
139   fPhiDistribution->SetParameter(3,dNewV4RP);
140    
141   // eta:
142   Double_t dEtaMin = -1.; // to be improved 
143   Double_t dEtaMax = 1.; // to be improved 
144   
145   Int_t iGoodTracks = 0;
146   Int_t iSelParticlesRP = 0;
147   Int_t iSelParticlesPOI = 0;
148   Double_t dTmpPt = 0.;
149   Double_t dTmpEta = 0.;
150   Double_t dTmpV2 = 0.;
151   Double_t dTmpPhi = 0.;
152   for(Int_t i=0;i<iNewMultiplicityOfRP;i++) {
153     dTmpPt = fPtSpectra->GetRandom();
154     dTmpEta = fMyTRandom3->Uniform(dEtaMin,dEtaMax);
155     // to be improved:
156     if(!fUseConstantHarmonics) {
157       if(dTmpPt >= fPtCutOff) { 
158         dTmpV2 = fV2RPMax;
159       } else {
160         dTmpV2 = fV2RPMax*(dTmpPt/fPtCutOff);
161       }  
162       fPhiDistribution->SetParameter(1,dTmpV2);         
163     }
164     dTmpPhi = fPhiDistribution->GetRandom();
165     // add the track to the event
166     for(Int_t d=0;d<fNoOfLoops;d++) {
167       AliFlowTrackSimple* pTrack = new AliFlowTrackSimple();
168       pTrack->SetPt(dTmpPt);
169       pTrack->SetEta(dTmpEta);
170       pTrack->SetPhi(dTmpPhi);
171       pTrack->SetForRPSelection(kTRUE);
172       iSelParticlesRP++;
173       pTrack->SetForPOISelection(kTRUE);
174       iSelParticlesPOI++;
175       pEvent->TrackCollection()->Add(pTrack);
176       iGoodTracks++;
177     }
178   }
179   
180   // update the event quantities
181   pEvent->SetEventNSelTracksRP(iSelParticlesRP);  
182   pEvent->SetNumberOfTracks(iGoodTracks);//tracks used either for RP or for POI selection
183   pEvent->SetMCReactionPlaneAngle(dMCReactionPlaneAngle);
184
185   if (!dMCReactionPlaneAngle == 0) cout<<" MC Reaction Plane Angle = "<<  dMCReactionPlaneAngle << endl;
186   else cout<<" MC Reaction Plane Angle = unknown "<< endl;
187
188   cout<<" iGoodTracks = "<< iGoodTracks << endl;
189   cout<<" # of RP selected tracks = "<<iSelParticlesRP<<endl;
190   cout<<" # of POI selected tracks = "<<iSelParticlesPOI<<endl;  
191   cout << "# " << ++fCount << " events processed" << endl;
192   
193   return pEvent;  
194  
195 } // end of CreateEventOnTheFly()
196
197
198