Updated PaintContour() method
[u/mrichter/AliRoot.git] / HBTAN / AliHBTPositionRandomizer.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 // 
21 // class AliHBTPositionRandomizer
22 //
23 // These class randomizes particle vertex positions
24 // Piotr.Skowronski@cern.ch
25 //
26 ////////////////////////////////////////////////////////////////////////////////
27
28 #include <TMath.h>
29 #include <TRandom.h>
30
31 #include "AliAOD.h"
32 #include "AliHBTPositionRandomizer.h"
33 #include "AliLog.h"
34 #include "AliVAODParticle.h"
35
36
37 ClassImp(AliHBTPositionRandomizer)
38
39 const Int_t AliHBTPositionRandomizer::fgkNumberOfPids = 10;
40
41 /*********************************************************************/
42
43 AliHBTPositionRandomizer::AliHBTPositionRandomizer():
44  fReader(0x0),
45  fDefaultRandomizer(0x0),
46  fRandomizers(0x0),
47  fNPid(0),
48  fPids(0x0),
49  fAddToExistingPos(kFALSE),
50  fOnlyParticlesFromVertex(kFALSE),
51  fRandomizeTracks(kFALSE),
52  fVX(0.0),
53  fVY(0.0),
54  fVZ(0.0)
55 {
56 //constructor
57 }
58 /*********************************************************************/
59
60 AliHBTPositionRandomizer::AliHBTPositionRandomizer(AliReader* reader):
61  fReader(reader),
62  fRandomizers(new TObjArray(fgkNumberOfPids)),
63  fNPid(1),
64  fPids(new Int_t[fgkNumberOfPids]),
65  fAddToExistingPos(kFALSE),
66  fOnlyParticlesFromVertex(kFALSE),
67  fRandomizeTracks(kFALSE),
68  fVX(0.0),
69  fVY(0.0),
70  fVZ(0.0)
71 {
72 //constructor
73  fRandomizers->AddAt(new AliHBTRndmGaussBall(8.0,0.0,0.0),0);
74
75 /*********************************************************************/
76
77 AliHBTPositionRandomizer::AliHBTPositionRandomizer(const AliHBTPositionRandomizer& in):
78  AliReader(in),
79  fReader(),
80  fRandomizers(0x0),
81  fNPid(0),
82  fPids(0x0),
83  fAddToExistingPos(kFALSE),
84  fOnlyParticlesFromVertex(kFALSE),
85  fRandomizeTracks(kFALSE),
86  fVX(0.0),
87  fVY(0.0),
88  fVZ(0.0)
89 {
90   //cpy constructor
91   in.Copy(*this);
92 }
93 /*********************************************************************/
94 AliHBTPositionRandomizer::~AliHBTPositionRandomizer()
95 {
96   //dtor
97   delete fReader;
98   delete fRandomizers;
99   delete [] fPids;
100 }
101 /*********************************************************************/
102 AliHBTPositionRandomizer& AliHBTPositionRandomizer::operator=(const AliHBTPositionRandomizer& in)
103 {
104   //assigment operator
105   in.Copy(*this);
106   return *this;
107 }
108 /*********************************************************************/
109
110 AliAOD* AliHBTPositionRandomizer::GetEventSim() const
111 {
112  // gets from fReader and randomizes current particle event
113  if (fReader == 0x0) 
114   {
115     Error("GetEventSim","Reader is null");
116     return 0x0;
117   } 
118  AliAOD *e =  fReader->GetEventSim();
119  if (e) 
120    if (e->IsRandomized() == kFALSE) 
121      Randomize(e);
122  return e;
123 }
124 /*********************************************************************/
125
126 AliAOD* AliHBTPositionRandomizer::GetEventRec() const
127 {
128  // gets from fReader and randomizes current track event
129  if (fReader == 0x0) 
130   {
131     Error("GetEventRec","Reader is null");
132     return 0x0;
133   }  
134  AliAOD *e =  fReader->GetEventRec();
135  if (fRandomizeTracks && e) if (e->IsRandomized() == kFALSE) Randomize(e);
136  return e;
137 }
138 /*********************************************************************/
139
140 AliAOD* AliHBTPositionRandomizer::GetEventSim(Int_t n)
141 {
142 //returns event n
143  if (fReader == 0x0) return 0x0;
144  AliAOD *e =  fReader->GetEventSim(n);
145  if (e->IsRandomized() == kFALSE) Randomize(e);
146  return e;
147 }
148
149 /*********************************************************************/
150 void AliHBTPositionRandomizer::Randomize(AliAOD* event) const
151 {
152 // randomizes postions of all particles in the event
153   static const Double_t kfmtocm = 1.e-13;
154   AliDebug(5," ");
155   if (event == 0x0) return;
156
157   for (Int_t i = 0; i < event->GetNumberOfParticles(); i++)
158    {
159      AliVAODParticle* p = event->GetParticle(i);
160      Double_t x,y,z,t=0.0;
161      AliHBTRndm* r = GetRandomizer(p->GetPdgCode());
162      r->Randomize(x,y,z,t,p);
163      
164      Double_t nx = x*kfmtocm;
165      Double_t ny = y*kfmtocm;
166      Double_t nz = z*kfmtocm;
167      Double_t nt = t*kfmtocm;
168      
169      if (fAddToExistingPos)
170       {
171        nx += p->Vx();
172        ny += p->Vy();
173        nz += p->Vz();
174        nt += p->T();
175       }
176      p->SetProductionVertex(nx,ny,nz,nt); 
177    }
178   event->SetRandomized();
179 }
180 /*********************************************************************/
181
182 AliHBTRndm* AliHBTPositionRandomizer::GetRandomizer(Int_t pdg) const
183 {
184   //returns randomizer for a given pdg 
185   Int_t idx = GetRandomizerIndex(pdg);//in most of cases 
186   if (idx < 0) idx = 0;//if not found return a default one
187   return (AliHBTRndm*)fRandomizers->At(idx);
188 }
189 /*********************************************************************/
190 Int_t AliHBTPositionRandomizer::GetRandomizerIndex(Int_t pdg) const
191 {
192   //returns randomizer index for a given pdg 
193
194   if (pdg == 0) return 0;
195   
196   for (Int_t i=1; i < fNPid; i++)
197    {
198      if (fPids[i] == pdg) 
199       return i;
200    }
201    
202   return -1;
203 }
204 /*********************************************************************/
205
206 void AliHBTPositionRandomizer::SetRandomizer(Int_t pid, AliHBTRndm* rndm)
207 {
208  //sets the randomizer for a given particle type
209   if (rndm == 0x0)
210    {
211      Error("SetRandomizer","Randomizer is null");
212      return;
213    }
214    
215   Int_t idx = GetRandomizerIndex(pid);
216   if (idx >= 0) 
217    {
218      delete fRandomizers->At(idx);
219      fRandomizers->AddAt(rndm,idx);
220    }  
221   
222   if (fNPid == fgkNumberOfPids)
223    {
224      Error("SetRandomizer","There is no more space in the array");
225      return;
226    }
227
228   fPids[fNPid] = pid;
229   fRandomizers->AddAt(rndm,fNPid);
230   fNPid++;
231 }
232 /*********************************************************************/
233
234 void AliHBTPositionRandomizer::SetGaussianBall(Int_t pid, Double_t r, Double_t meantime, Double_t sigmatime)
235 {
236  //Sets Gaussian Ball Model
237   SetGaussianBall(pid,r,r,r,meantime,sigmatime);
238 }
239 /*********************************************************************/
240
241 void AliHBTPositionRandomizer::SetGaussianBall(Int_t pid, Double_t rx, Double_t ry, Double_t rz, Double_t meantime, Double_t sigmatime)
242 {
243  //Sets Gaussian Ball Model
244   AliHBTRndm* rndm = new AliHBTRndmGaussBall(rx,ry,rz,meantime,sigmatime);
245   SetRandomizer(pid,rndm);
246 }
247 /*********************************************************************/
248
249 void AliHBTPositionRandomizer::SetCyllinderSurface(Int_t pid, Double_t r, Double_t l)
250 {
251  //Sets Cylinder Surface Model
252   AliHBTRndm* rndm = new  AliHBTRndmCyllSurf(r,l);
253   SetRandomizer(pid,rndm);
254 }
255 /*********************************************************************/
256
257 void AliHBTPositionRandomizer::SetEventVertex(Double_t x, Double_t y,Double_t z)
258 {
259 //sets event vertex position
260   fVX = x;
261   fVY = y;
262   fVZ = z;
263 }
264
265
266 void AliHBTPositionRandomizer::SetEllipse(Int_t pid, Double_t rx, Double_t ryz)
267 {
268 //sets the ellipse randomization for the given pid
269   AliHBTRndm* rndm = new AliHBTRndmEllipse(rx,ryz);
270   SetRandomizer(pid,rndm);   
271 }
272
273 /*********************************************************************/
274 //_____________________________________________________________________
275 ///////////////////////////////////////////////////////////////////////
276 //                                                                   //
277 //  class AliHBTRndmGaussBall                                        //
278 //                                                                   //
279 ///////////////////////////////////////////////////////////////////////
280
281 AliHBTRndmGaussBall::AliHBTRndmGaussBall():
282  fRx(0.0),
283  fRy(0.0),
284  fRz(0.0),
285  fTmean(0.0),
286  fTsigma(0.0)
287 {
288   //constructor
289 }
290 /*********************************************************************/
291
292 AliHBTRndmGaussBall::AliHBTRndmGaussBall(Float_t r, Double_t meantime, Double_t sigmatime):
293  fRx(r),
294  fRy(r),
295  fRz(r),
296  fTmean(meantime),
297  fTsigma(sigmatime)
298 {
299   //constructor
300 }
301 /*********************************************************************/
302
303 AliHBTRndmGaussBall::AliHBTRndmGaussBall(Float_t rx, Float_t ry, Float_t rz, Double_t meantime, Double_t sigmatime):
304  fRx(rx),
305  fRy(ry),
306  fRz(rz),
307  fTmean(meantime),
308  fTsigma(sigmatime)
309 {
310   //constructor
311 }
312 /*********************************************************************/
313
314
315 AliHBTRndmEllipse::AliHBTRndmEllipse(Float_t rmin, Float_t rmax):
316  fRmin(rmin),
317  fRmax(rmax)
318 {
319      //constructor
320 }
321
322 /*********************************************************************/
323
324 void AliHBTRndmGaussBall::Randomize(Double_t& x,Double_t& y,Double_t&z,Double_t&t, AliVAODParticle*/*particle*/) const
325 {
326 //randomizez gauss for each coordinate separately
327   x = gRandom->Gaus(0.0,fRx);
328   y = gRandom->Gaus(0.0,fRy);
329   z = gRandom->Gaus(0.0,fRz);
330   
331   if (fTsigma == 0.0)
332    {
333      t = fTmean;
334      return;
335    }
336   
337   t = gRandom->Gaus(fTmean,fTsigma);
338     
339 }
340 /*********************************************************************/
341 //_____________________________________________________________________
342 ///////////////////////////////////////////////////////////////////////
343 //                                                                   //
344 //  class AliHBTRndmGaussBall                                        //
345 //                                                                   //
346 ///////////////////////////////////////////////////////////////////////
347
348 void AliHBTRndmCyllSurf::Randomize(Double_t& x,Double_t& y,Double_t&z,Double_t&/*t*/, AliVAODParticle* particle) const
349 {
350 //Randomizes x,y,z
351    Double_t r = fR + gRandom->Gaus(0.0, 1.0);
352    Double_t sf = r/particle->Pt();//scaling factor for position transformation ->
353                              //we move direction of string momentum but legth defined by r
354    x = sf*particle->Px();
355    y = sf*particle->Py();
356    z = gRandom->Uniform(-fL,fL);
357 }
358
359 /*********************************************************************/
360 /*********************************************************************/
361
362 void AliHBTRndmEllipse::Randomize(Double_t& x, Double_t& y, Double_t& z,Double_t&/*t*/, AliVAODParticle*p) const
363 {
364     // p=0; //workaround - fix this damn little thingy
365    double R;
366      double phi=p->Phi();
367      
368      R=fRmin+(fRmax-fRmin)*TMath::Sin(phi);
369      x=R*TMath::Sin(phi);
370      y=R*TMath::Cos(phi);
371      z=z;
372 }