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