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