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
  31 #include "AliFlowEventSimpleMakerOnTheFly.h"
  32 #include "AliFlowEventSimple.h"
  33 #include "AliFlowTrackSimple.h"
  34
  35 ClassImp(AliFlowEventSimpleMakerOnTheFly)
  36
  37
  38 //========================================================================
  39
  40
  41 AliFlowEventSimpleMakerOnTheFly::AliFlowEventSimpleMakerOnTheFly(UInt_t iseed):
  42   fUseGlauberModel(kFALSE),
  43   fMultDistrOfRPsIsGauss(kFALSE),
  44   fMultiplicityOfRP(0),
  45   fMultiplicitySpreadOfRP(0.),
  46   fMinMultOfRP(0),
  47   fMaxMultOfRP(0),
  48   fTemperatureOfRP(0.),
  49   fPtDependentHarmonicV1(kFALSE),
  50   fEtaDependentHarmonicV1(kFALSE),
  51   fPtDependentHarmonicV2(kFALSE),
  52   fEtaDependentHarmonicV2(kFALSE),
  53   fPtDependentHarmonicV4(kFALSE),
  54   fEtaDependentHarmonicV4(kFALSE),
  55   fV1RP(0.),
  56   fV1SpreadRP(0.),
  57   fConstantV2IsSampledFromGauss(kFALSE),
  58   fV2RP(0.),
  59   fV2SpreadRP(0.),
  60   fMinV2RP(0.),
  61   fMaxV2RP(0.),
  62   fV3RP(0.),
  63   fV3SpreadRP(0.),
  64   fV4RP(0.),
  65   fV4SpreadRP(0.),
  66   fV1vsPtEtaMax(0.),
  67   fV1PtCutOff(0.),
  68   fV2vsPtEtaMax(0.),
  69   fV2PtCutOff(0.),
  70   fV2vsEtaSpread(0.),
  71   fPhiMin1(0.),
  72   fPhiMax1(0.),
  73   fProbability1(0.),
  74   fPhiMin2(0.),
  75   fPhiMax2(0.),
  76   fProbability2(0.),
  77   fPtSpectra(NULL),
  78   fPhiDistribution(NULL),
  79   fMyTRandom3(NULL),
  80   fCount(0),
  81   fNoOfLoops(1),
  82   fPhiRange(0.),
  83   fPtRange(0.),
  84   fEtaRange(0.),
  85   fNonflowSectorMin(0.),
  86   fNonflowSectorMax(TMath::TwoPi()),
  87   fEtaMinA(-1.0),
  88   fEtaMaxA(-0.01),
  89   fEtaMinB(0.01),
  90   fEtaMaxB(1.0),
  91   fCreateJets(kFALSE),
  92   fJetProbability(0.25),
  93   fJetTracksFraction(0.1),
  94   fJetCone(10.)
  95 {
  96   // constructor
  97   fMyTRandom3 = new TRandom3(iseed);
  98   gRandom->SetSeed(fMyTRandom3->Integer(65539));
  99 }
 100
 101
 102 //========================================================================
 103
 104
 105 AliFlowEventSimpleMakerOnTheFly::~AliFlowEventSimpleMakerOnTheFly()
 106 {
 107  // destructor
 108  if (fPtSpectra) delete fPtSpectra;
 109  if (fPhiDistribution) delete fPhiDistribution;
 110  if (fMyTRandom3) delete fMyTRandom3;
 111 }
 112
 113
 114 //========================================================================
 115
 116
 117 void AliFlowEventSimpleMakerOnTheFly::Init()
 118 {
 119  // define the pt spectra and phi distribution
 120
 121  // pt spectra of pions (Boltzman):
 122  Double_t dPtMin = 0.; // to be improved (move this to the body of contstructor?)
 123  Double_t dPtMax = 10.; // to be improved (move this to the body of contstructor?)
 124
 125  fPtSpectra = new TF1("fPtSpectra","[0]*x*TMath::Exp(-pow(0.13957*0.13957+x*x,0.5)/[1])",dPtMin,dPtMax);
 126  fPtSpectra->SetParName(0,"Multiplicity of RPs");
 127  fPtSpectra->SetParName(1,"Temperature of RPs");
 128
 129  // phi distribution:
 130  Double_t dPhiMin = 0.; // to be improved (move this to the body of contstructor?)
 131  Double_t dPhiMax = TMath::TwoPi(); // to be improved (move this to the body of contstructor?)
 132
 133  fPhiDistribution = new TF1("fPhiDistribution","1+2.*[0]*TMath::Cos(x-[2])+2.*[1]*TMath::Cos(2*(x-[2]))+2.*[4]*TMath::Cos(3*(x-[2]))+2.*[3]*TMath::Cos(4*(x-[2]))",dPhiMin,dPhiMax);
 134  fPhiDistribution->SetParName(0,"directed flow");
 135  fPhiDistribution->SetParName(1,"elliptic flow");
 136  fPhiDistribution->SetParName(2,"Reaction Plane");
 137  fPhiDistribution->SetParName(4,"triangular flow");
 138  fPhiDistribution->SetParName(3,"harmonic 4"); // to be improved (name)
 139
 140 }
 141
 142 //========================================================================
 143
 144 Int_t AliFlowEventSimpleMakerOnTheFly::DetermineMultiplicity()
 145 {
 146  // Determine multiplicity for current event.
 147
 148  Int_t iNewMultiplicityOfRP = fMultiplicityOfRP;
 149
 150  if(fMultDistrOfRPsIsGauss) {
 151     if (fMultiplicitySpreadOfRP>0.0) iNewMultiplicityOfRP = (Int_t)fMyTRandom3->Gaus(fMultiplicityOfRP,fMultiplicitySpreadOfRP);
 152     fPtSpectra->SetParameter(0,iNewMultiplicityOfRP);
 153   } else {
 154     if (fMinMultOfRP != fMaxMultOfRP) {
 155       iNewMultiplicityOfRP = (Int_t)fMyTRandom3->Uniform(fMinMultOfRP,fMaxMultOfRP);
 156       fPtSpectra->SetParameter(0,iNewMultiplicityOfRP);
 157     } else {
 158       fPtSpectra->SetParameter(0,fMinMultOfRP);
 159     }
 160   }
 161
 162  return iNewMultiplicityOfRP;
 163
 164 } // end of AliFlowEventSimpleMakerOnTheFly::DetermineMultiplicity()
 165
 166 //========================================================================
 167
 168 void AliFlowEventSimpleMakerOnTheFly::DetermineV1()
 169 {
 170  // Determine flow harmonics v1 for current event (if v1 is not pt or eta dependent).
 171
 172  Double_t dNewV1RP=fV1RP;
 173  if(fV1SpreadRP>0.0) {dNewV1RP = fMyTRandom3->Gaus(fV1RP,fV1SpreadRP);}
 174  fPhiDistribution->SetParameter(0,dNewV1RP);
 175
 176 } // end of void AliFlowEventSimpleMakerOnTheFly::DetermineV1()
 177
 178 //========================================================================
 179
 180 void AliFlowEventSimpleMakerOnTheFly::DetermineV3()
 181 {
 182  // Determine flow harmonics V3 for current event (if V3 is not pt or eta dependent).
 183
 184  Double_t dNewV3RP = fV3RP;
 185  if(fV3SpreadRP>0.0) dNewV3RP = fMyTRandom3->Gaus(fV3RP,fV3SpreadRP);
 186  fPhiDistribution->SetParameter(4,dNewV3RP);
 187
 188 } // end of void AliFlowEventSimpleMakerOnTheFly::DetermineV3()
 189
 190 //========================================================================
 191
 192 void AliFlowEventSimpleMakerOnTheFly::DetermineV4()
 193 {
 194  // Determine flow harmonics v4 for current event (if v4 is not pt or eta dependent).
 195
 196  Double_t dNewV4RP = fV4RP;
 197  if(fV4SpreadRP>0.0) dNewV4RP = fMyTRandom3->Gaus(fV4RP,fV4SpreadRP);
 198  fPhiDistribution->SetParameter(3,dNewV4RP);
 199
 200 } // end of void AliFlowEventSimpleMakerOnTheFly::DetermineV4()
 201
 202 //========================================================================
 203
 204 void AliFlowEventSimpleMakerOnTheFly::DetermineV2()
 205 {
 206  // Determine flow harmonics v2 for current event (if v2 is not pt or eta dependent).
 207
 208  Double_t dNewV2RP = fV2RP;
 209
 210  if(fConstantV2IsSampledFromGauss)
 211  {
 212   if(fV2SpreadRP>0.0)
 213   {
 214    dNewV2RP = fMyTRandom3->Gaus(fV2RP,fV2SpreadRP);}
 215   fPhiDistribution->SetParameter(1,dNewV2RP);
 216  } else if(fMinV2RP < fMaxV2RP)
 217    {
 218     dNewV2RP = fMyTRandom3->Uniform(fMinV2RP,fMaxV2RP);
 219     fPhiDistribution->SetParameter(1,dNewV2RP);
 220    } else if(fMinV2RP == fMaxV2RP)
 221      {
 222       dNewV2RP = fMinV2RP;
 223       fPhiDistribution->SetParameter(1,dNewV2RP);
 224      }
 225
 226 } // end of void AliFlowEventSimpleMakerOnTheFly::DetermineV2()
 227
 228 //========================================================================
 229
 230 Int_t AliFlowEventSimpleMakerOnTheFly::GlauberModel()
 231 {
 232  // Determine multiplicity and flow harmonics for current event from Glauber moder
 233
 234  Int_t multiplicity = 0;
 235  Double_t v1 = 0.;
 236  Double_t v2 = 0.;
 237  Double_t v3 = 0.;
 238  Double_t v4 = 0.;
 239
 240  // Determine multiplicity, v1, v2, v3 and v4 from Glauber model:
 241
 242  // multiplicity = ...
 243  // v1 = ...
 244  // v2 = ...
 245  // v3 = ...
 246  // v4 = ...
 247
 248  // Set obtained values as parameters in relevant distributions:
 249  fPtSpectra->SetParameter(0,multiplicity);
 250  fPhiDistribution->SetParameter(0,v1);
 251  fPhiDistribution->SetParameter(1,v2);
 252  fPhiDistribution->SetParameter(4,v3);
 253  fPhiDistribution->SetParameter(3,v4);
 254
 255  return multiplicity;
 256
 257 } // end of Int_t AliFlowEventSimpleMakerOnTheFly::GlauberModel()
 258
 259 //========================================================================
 260
 261 AliFlowEventSimple* AliFlowEventSimpleMakerOnTheFly::CreateEventOnTheFly(AliFlowTrackSimpleCuts *cutsRP, AliFlowTrackSimpleCuts *cutsPOI)
 262 {
 263   // Method to create event on the fly.
 264
 265   // Determine multiplicity and flow harmonics (if not pt or eta dependent) for current event:
 266   Int_t multiplicityRP = 0;
 267   if(!fUseGlauberModel)
 268   {
 269    multiplicityRP = DetermineMultiplicity();
 270    if(!(fPtDependentHarmonicV1||fEtaDependentHarmonicV1)) {DetermineV1();}
 271    if(!(fPtDependentHarmonicV2||fEtaDependentHarmonicV2)) {DetermineV2();}
 272    DetermineV3(); // to be improved - add also pt and eta dependence for v3
 273    if(!(fPtDependentHarmonicV4||fEtaDependentHarmonicV4)) {DetermineV4();}
 274   } else
 275     {
 276      // Determine multipliciy and flow harmonics from Glauber model:
 277      multiplicityRP = GlauberModel();
 278     }
 279
 280   AliFlowEventSimple* pEvent = new AliFlowEventSimple(multiplicityRP);
 281
 282   // decide if you will create jets or not in this event:
 283   if(fCreateJets)
 284   {
 285    if(fMyTRandom3->Uniform(0,1) > 1 - fJetProbability)
 286    {
 287     fCreateJets = kTRUE;
 288    } else
 289      {
 290       fCreateJets = kFALSE;
 291      }
 292   }
 293
 294   // set the 'temperature' of RPs
 295   fPtSpectra->SetParameter(1,fTemperatureOfRP);
 296
 297   // sampling the reaction plane
 298   Double_t dMCReactionPlaneAngle = fMyTRandom3->Uniform(0.,TMath::TwoPi());
 299   fPhiDistribution->SetParameter(2,dMCReactionPlaneAngle);
 300
 301   // eta:
 302   Double_t dEtaMin = -1.; // to be improved
 303   Double_t dEtaMax = 1.; // to be improved
 304
 305   Int_t iGoodTracks = 0;
 306   Int_t iSelParticlesRP = 0;
 307   Int_t iSelParticlesPOI = 0;
 308   // parameters of original tracks:
 309   Double_t dPhiOriginalTrack = 0.;
 310   Double_t dPtOriginalTrack = 0.;
 311   Double_t dEtaOriginalTrack = 0.;
 312   // parameters of splitted tracks:
 313   Double_t dPhiSplittedTrack = 0.;
 314   Double_t dPtSplittedTrack = 0.;
 315   Double_t dEtaSplittedTrack = 0.;
 316
 317   Double_t dTmpV1 = 0.;
 318   Double_t dTmpV2 = 0.;
 319   //Double_t dTmpV3 = 0.;
 320   Double_t dTmpV4 = 0.;
 321   Bool_t bUniformAcceptance = kTRUE;
 322   Double_t Pi = TMath::Pi();
 323
 324   if(!((fPhiMin1==0.) && (fPhiMax1==0.) && (fPhiMin2==0.) && (fPhiMax2==0.))) {
 325     bUniformAcceptance = kFALSE;
 326   }
 327   // loop over original tracks:
 328   for(Int_t i=0;i<multiplicityRP;i++)
 329   {
 330    // sample the pt and eta for original track:
 331    dPtOriginalTrack = fPtSpectra->GetRandom();
 332    dEtaOriginalTrack = fMyTRandom3->Uniform(dEtaMin,dEtaMax);
 333    // generate flow harmonics which will determine the azimuthal distribution (to be improved - optimized):
 334    // V2:
 335    if(fPtDependentHarmonicV2 || fEtaDependentHarmonicV2)
 336    {
 337     if(fEtaDependentHarmonicV2)
 338     {
 339      if(fV2vsEtaSpread>0.)
 340      {
 341       dTmpV2 = TMath::Exp(-pow(dEtaOriginalTrack/fV2vsEtaSpread,2.));
 342      }
 343      if(!fPtDependentHarmonicV2)
 344      {
 345       dTmpV2*=fV2vsPtEtaMax;
 346      }
 347     } // end of if(fEtaDependentHarmonicV2)
 348     if(fPtDependentHarmonicV2)
 349     {
 350      if(!fEtaDependentHarmonicV2)
 351      {
 352       if(dPtOriginalTrack >= fV2PtCutOff) {dTmpV2 = fV2vsPtEtaMax;}
 353       else {dTmpV2 = fV2vsPtEtaMax*(dPtOriginalTrack/fV2PtCutOff);}
 354      } else
 355        {
 356         if(dPtOriginalTrack >= fV2PtCutOff) {dTmpV2 *= fV2vsPtEtaMax;}
 357         else {dTmpV2 *= fV2vsPtEtaMax*(dPtOriginalTrack/fV2PtCutOff);}
 358        }
 359     } // end of if(fPtDependentHarmonicV2)
 360     // flow harmonic is determined and plugged in as a parameter in the predefined azimuthal distribution:
 361     fPhiDistribution->SetParameter(1,dTmpV2);
 362    }
 363    // V1:
 364    if(fPtDependentHarmonicV1 || fEtaDependentHarmonicV1)
 365    {
 366     if(fEtaDependentHarmonicV1)
 367     {
 368      dTmpV1 = -1.*dEtaOriginalTrack;
 369      if(!fPtDependentHarmonicV1)
 370      {
 371       dTmpV1*=fV1vsPtEtaMax;
 372      }
 373     } // end of if(fEtaDependentHarmonicV1)
 374     if(fPtDependentHarmonicV1)
 375     {
 376      if(!fEtaDependentHarmonicV1)
 377      {
 378       if(dPtOriginalTrack >= fV1PtCutOff) {dTmpV1 = fV1vsPtEtaMax;}
 379       else {dTmpV1 = fV1vsPtEtaMax*(dPtOriginalTrack/fV1PtCutOff);}
 380      } else
 381        {
 382         if(dPtOriginalTrack >= fV1PtCutOff) {dTmpV1 *= fV1vsPtEtaMax;}
 383         else {dTmpV1 *= fV1vsPtEtaMax*(dPtOriginalTrack/fV1PtCutOff);}
 384        }
 385     } // end of if(fPtDependentHarmonicV1)
 386     // flow harmonic is determined and plugged in as a parameter in the predefined azimuthal distribution:
 387     fPhiDistribution->SetParameter(0,dTmpV1);
 388    } // end of if(fPtDependentHarmonicV1 || fEtaDependentHarmonicV1)
 389    // V4:
 390    if(fPtDependentHarmonicV4 || fEtaDependentHarmonicV4)
 391    {
 392     dTmpV4 = pow(dTmpV2,2.);
 393     fPhiDistribution->SetParameter(3,dTmpV4);
 394    }
 395    // sample the phi angle for original track:
 396    dPhiOriginalTrack = fPhiDistribution->GetRandom();
 397    // from each original track make fNoOfLoops splitted tracks if the particle is ongoing in
 398    // detector's sector ranging from fNonflowSectorMin to fNonflowSectorMax
 399    // (simulating nonflow correlations between fNoOfLoops tracks in certain detector's sector):
 400    for(Int_t d=0;d<fNoOfLoops;d++)
 401    {
 402     if(d>0 && (dPhiOriginalTrack>=fNonflowSectorMin && dPhiOriginalTrack<fNonflowSectorMax))
 403     {
 404      dPhiSplittedTrack = dPhiOriginalTrack;
 405      dPtSplittedTrack = dPtOriginalTrack;
 406      dEtaSplittedTrack = dEtaOriginalTrack;
 407      // phi:
 408      if(fPhiRange>0.)
 409      {
 410       dPhiSplittedTrack = fMyTRandom3->Uniform(dPhiOriginalTrack-fPhiRange,dPhiOriginalTrack+fPhiRange);
 411       if(dPhiSplittedTrack<0)
 412       {
 413        dPhiSplittedTrack+=TMath::TwoPi(); // to ensure angle is in [0,2Pi>
 414       }
 415       if(dPhiSplittedTrack>=TMath::TwoPi())
 416       {
 417        dPhiSplittedTrack-=TMath::TwoPi(); // to ensure angle is in [0,2Pi>
 418       }
 419      } // end of if(fPhiRange>0.)
 420      // pt:
 421      if(fPtRange>0.)
 422      {
 423       Double_t minPt = dPtOriginalTrack-fPtRange;
 424       Double_t maxPt = dPtOriginalTrack+fPtRange;
 425       if(minPt<0)
 426       {
 427        minPt = 0.; // protection against pt<0 for splitted track
 428       }
 429       dPtSplittedTrack = fMyTRandom3->Uniform(minPt,maxPt);
 430      } // end of if(fPtRange>0.)
 431      // eta:
 432      if(fEtaRange>0.)
 433      {
 434       dEtaSplittedTrack = fMyTRandom3->Uniform(dEtaOriginalTrack-fEtaRange,dEtaOriginalTrack+fEtaRange);
 435      } // end of if(fEtaRange>0.)
 436     } // end of if(d>0 && (dPhiOriginalTrack>=fNonflowSectorMin && dPhiOriginalTrack<fNonflowSectorMax))
 437     Double_t dTmpPhi = -44.;
 438     Double_t dTmpPt = -44.;
 439     Double_t dTmpEta = -44.;
 440     if(d>0)
 441     {
 442      if(dPhiOriginalTrack>=fNonflowSectorMin && dPhiOriginalTrack<fNonflowSectorMax)
 443      {
 444       dTmpPhi = dPhiSplittedTrack;
 445       dTmpPt = dPtSplittedTrack;
 446       dTmpEta = dEtaSplittedTrack;
 447      }
 448     } else
 449       {
 450        dTmpPhi = dPhiOriginalTrack;
 451        dTmpPt = dPtOriginalTrack;
 452        dTmpEta = dEtaOriginalTrack;
 453       }
 454     // make the new track:
 455     AliFlowTrackSimple *pTrack = new AliFlowTrackSimple();
 456     // uniform acceptance:
 457     if(bUniformAcceptance)
 458     {
 459      if(!fCreateJets || (fCreateJets && !(i >= TMath::Nint((1.-fJetTracksFraction)*multiplicityRP))))
 460      {
 461       pTrack->SetPt(dTmpPt);
 462       pTrack->SetEta(dTmpEta);
 463       pTrack->SetPhi(dTmpPhi);
 464       // checking RP cuts:
 465       if(cutsRP->PassesCuts(pTrack))
 466       {
 467        pTrack->SetForRPSelection(kTRUE);
 468        iSelParticlesRP++;
 469       }
 470       // assign particles to subevents:
 471       if(pTrack->Eta()>=fEtaMinA && pTrack->Eta()<=fEtaMaxA)
 472       {
 473        pTrack->SetForSubevent(0);
 474       }
 475       if(pTrack->Eta()>=fEtaMinB && pTrack->Eta()<=fEtaMaxB)
 476       {
 477        pTrack->SetForSubevent(1);
 478       }
 479       // checking POI cuts:
 480       if(cutsPOI->PassesCuts(pTrack))
 481       {
 482        pTrack->SetForPOISelection(kTRUE);
 483        iSelParticlesPOI++;
 484       }
 485       pEvent->AddTrack(pTrack);
 486      } else // to if(!fCreateJets || (fCreateJets && !(i >= (Int_t)(1.-fJetTracksFraction)*multiplicityRP)))
 487        {
 488         for(Int_t j=0;j<TMath::Nint(fJetTracksFraction*multiplicityRP);j++) // fragmenting last sampled particle into the jet
 489         {
 490          AliFlowTrackSimple *pTrackInJet = new AliFlowTrackSimple();
 491          Double_t dTmpPhiWithinJet = fMyTRandom3->Uniform(dTmpPhi-fJetCone*TMath::Pi()/360.,dTmpPhi+fJetCone*TMath::Pi()/360.);
 492          if(dTmpPhiWithinJet<0.)
 493          {
 494           dTmpPhiWithinJet+=TMath::TwoPi(); // to ensure angle is in [0,2Pi>
 495          }
 496          if(dTmpPhiWithinJet>=TMath::TwoPi())
 497          {
 498           dTmpPhiWithinJet-=TMath::TwoPi(); // to ensure angle is in [0,2Pi>
 499          }
 500          pTrackInJet->SetPt(dTmpPt); // to be improved - smear pt a little bit for particles in jet
 501          pTrackInJet->SetEta(dTmpEta); // to be improved - smear eta a little bit for particles in jet
 502          pTrackInJet->SetPhi(dTmpPhiWithinJet);
 503          // checking RP cuts:
 504          if(cutsRP->PassesCuts(pTrackInJet))
 505          {
 506           pTrackInJet->SetForRPSelection(kTRUE);
 507           iSelParticlesRP++;
 508          }
 509          // assign particles to subevents:
 510          if(pTrackInJet->Eta()>=fEtaMinA && pTrackInJet->Eta()<=fEtaMaxA)
 511          {
 512           pTrackInJet->SetForSubevent(0);
 513          }
 514          if(pTrackInJet->Eta()>=fEtaMinB && pTrackInJet->Eta()<=fEtaMaxB)
 515          {
 516           pTrackInJet->SetForSubevent(1);
 517          }
 518          // checking POI cuts:
 519          if(cutsPOI->PassesCuts(pTrackInJet))
 520          {
 521           pTrackInJet->SetForPOISelection(kTRUE);
 522           iSelParticlesPOI++;
 523          }
 524          pEvent->AddTrack(pTrackInJet);
 525          iGoodTracks++;
 526          //delete pTrackInJet;
 527         } // end of for(Int_t j=0;j<Int_t(fJetTracksFraction*multiplicityRP);j++)
 528         i=multiplicityRP;
 529        } // end of else // to if(!fCreateJets || (fCreateJets && !(i >= (Int_t)(1.-fJetTracksFraction)*multiplicityRP)))
 530     } // end of if(bUniformAcceptance)
 531     // non-uniform acceptance, 1st sector:
 532     else if((dTmpPhi > fPhiMin1*Pi/180) && (dTmpPhi < fPhiMax1*Pi/180))
 533     {
 534      if(fMyTRandom3->Uniform(0,1) > 1 - fProbability1)
 535      {
 536       pTrack->SetPt(dTmpPt);
 537       pTrack->SetEta(dTmpEta);
 538       pTrack->SetPhi(dTmpPhi);
 539       // checking RP cuts:
 540       if(cutsRP->PassesCuts(pTrack))
 541       {
 542        pTrack->SetForRPSelection(kTRUE);
 543        iSelParticlesRP++;
 544       }
 545       // assign particles to subevents
 546       if(pTrack->Eta()>=fEtaMinA && pTrack->Eta()<=fEtaMaxA)
 547       {
 548        pTrack->SetForSubevent(0);
 549       }
 550       if(pTrack->Eta()>=fEtaMinB && pTrack->Eta()<=fEtaMaxB)
 551       {
 552        pTrack->SetForSubevent(1);
 553       }
 554       // checking POI cuts:
 555       if(cutsPOI->PassesCuts(pTrack))
 556       {
 557        pTrack->SetForPOISelection(kTRUE);
 558        iSelParticlesPOI++;
 559       }
 560       pEvent->AddTrack(pTrack);
 561       iGoodTracks++;
 562      } // end of if(fMyTRandom3->Uniform(0,1) > 1 - fProbability1)
 563     } // end of else if ((dTmpPhi > fPhiMin1*Pi/180) && (dTmpPhi < fPhiMax1*Pi/180))
 564     // non-uniform acceptance, 2nd sector:
 565     else if((dTmpPhi > fPhiMin2*Pi/180) && (dTmpPhi < fPhiMax2*Pi/180))
 566     {
 567      if(fMyTRandom3->Uniform(0,1) > 1 - fProbability2)
 568      {
 569       pTrack->SetPt(dTmpPt);
 570       pTrack->SetEta(dTmpEta);
 571       pTrack->SetPhi(dTmpPhi);
 572       // checking RP cuts:
 573       if(cutsRP->PassesCuts(pTrack))
 574       {
 575        pTrack->SetForRPSelection(kTRUE);
 576        iSelParticlesRP++;
 577       }
 578       // assign particles to subevents
 579       if(pTrack->Eta()>=fEtaMinA && pTrack->Eta()<=fEtaMaxA)
 580       {
 581        pTrack->SetForSubevent(0);
 582       }
 583       if(pTrack->Eta()>=fEtaMinB && pTrack->Eta()<=fEtaMaxB)
 584       {
 585        pTrack->SetForSubevent(1);
 586       }
 587       // checking POI cuts:
 588       if(cutsPOI->PassesCuts(pTrack))
 589       {
 590        pTrack->SetForPOISelection(kTRUE);
 591        iSelParticlesPOI++;
 592       }
 593       pEvent->AddTrack(pTrack);
 594       iGoodTracks++;
 595      } // end of if(fMyTRandom3->Uniform(0,1) > 1 - fProbability2)
 596     } // end of else if ((dTmpPhi > fPhiMin2*Pi/180) && (dTmpPhi < fPhiMax2*Pi/180))
 597     else
 598     {
 599      pTrack->SetPt(dTmpPt);
 600      pTrack->SetEta(dTmpEta);
 601      pTrack->SetPhi(dTmpPhi);
 602      // checking RP cuts:
 603      if(cutsRP->PassesCuts(pTrack))
 604      {
 605       pTrack->SetForRPSelection(kTRUE);
 606       iSelParticlesRP++;
 607      }
 608      // assign particles to subevents
 609      if (pTrack->Eta()>=fEtaMinA && pTrack->Eta()<=fEtaMaxA)
 610      {
 611       pTrack->SetForSubevent(0);
 612      }
 613      if (pTrack->Eta()>=fEtaMinB && pTrack->Eta()<=fEtaMaxB)
 614      {
 615       pTrack->SetForSubevent(1);
 616      }
 617      // checking POI cuts:
 618      if(cutsPOI->PassesCuts(pTrack))
 619      {
 620       pTrack->SetForPOISelection(kTRUE);
 621       iSelParticlesPOI++;
 622      }
 623      pEvent->AddTrack(pTrack);
 624      iGoodTracks++;
 625     } // end of else
 626     //delete pTrack;
 627    } // end of for(Int_t d=0;d<fNoOfLoops;d++)
 628   } // end of for(Int_t i=0;i<iNewMultiplicityOfRP;i++)
 629
 630   // update the event quantities
 631   pEvent->SetEventNSelTracksRP(iSelParticlesRP);
 632   pEvent->SetMCReactionPlaneAngle(dMCReactionPlaneAngle);
 633
 634   Int_t cycle = 0;
 635   if(fPtDependentHarmonicV1 || fEtaDependentHarmonicV1 ||
 636      fPtDependentHarmonicV2 || fEtaDependentHarmonicV2 ||
 637      fPtDependentHarmonicV4 || fEtaDependentHarmonicV4)
 638   {
 639    cycle = 10;
 640   } else
 641     {
 642      cycle = 100;
 643     }
 644
 645   if ( (++fCount % cycle) == 0) {
 646     if (!dMCReactionPlaneAngle == 0) cout<<" MC Reaction Plane Angle = "<<  dMCReactionPlaneAngle << endl;
 647     else cout<<" MC Reaction Plane Angle = unknown "<< endl;
 648     cout<<" iGoodTracks = "<< iGoodTracks << endl;
 649     cout<<" # of RP selected tracks = "<<iSelParticlesRP<<endl;
 650     cout<<" # of POI selected tracks = "<<iSelParticlesPOI<<endl;
 651     cout << "# " << fCount << " events processed" << endl;
 652   }
 653
 654   return pEvent;
 655
 656 } // end of CreateEventOnTheFly()
 657
 658 //========================================================================
 659
 660 /*
 661 Double_t AliFlowEventSimpleMakerOnTheFly::NonflowBias()
 662 {
 663  // Determine multiplicity and flow harmonics for current event from Glauber moder
 664
 665  Int_t multiplicity = 0;
 666  Double_t v1 = 0.;
 667  Double_t v2 = 0.;
 668  Double_t v4 = 0.;
 669
 670  // Determine multiplicity, v1, v2 and v4 from Glauber model:
 671
 672  // multiplicity = ...
 673  // v1 = ...
 674  // v2 = ...
 675  // v4 = ...
 676
 677  // Set obtained values as parameters in relevant distributions:
 678  fPtSpectra->SetParameter(0,multiplicity);
 679  fPhiDistribution->SetParameter(0,v1);
 680  fPhiDistribution->SetParameter(1,v2);
 681  fPhiDistribution->SetParameter(3,v4);
 682
 683  return multiplicity;
 684
 685 } // end of Double_t AliFlowEventSimpleMakerOnTheFly::NonflowBias()
 686
 687 //========================================================================
 688
 689 */