HBTAN/AliHBTPositionRandomizer.cxx

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