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