PWG/FLOW/Base/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 full *
  18  * flow analysis 'on the fly'.      *
  19  *                                  *
  20  * author: Ante Bilandzic           *
  21  *         (abilandzic@gmail.com)   *
  22  ************************************/
  23
  24 #include "Riostream.h"
  25 #include "TMath.h"
  26 #include "TF1.h"
  27 #include "TRandom3.h"
  28 #include "AliFlowEventSimpleMakerOnTheFly.h"
  29 #include "AliFlowEventSimple.h"
  30 #include "AliFlowTrackSimple.h"
  31 #include "AliFlowTrackSimpleCuts.h"
  32
  33 ClassImp(AliFlowEventSimpleMakerOnTheFly)
  34
  35 //========================================================================================================================================
  36
  37 AliFlowEventSimpleMakerOnTheFly::AliFlowEventSimpleMakerOnTheFly(UInt_t uiSeed):
  38 fCount(0),
  39 fMinMult(0),
  40 fMaxMult(0),
  41 fPtSpectra(NULL),
  42 fMass(0.13957),
  43 fTemperature(0.44),
  44 fPhiDistribution(NULL),
  45 fV1(0.),
  46 fV2(0.05),
  47 fV3(0.),
  48 fV4(0.),
  49 fUniformFluctuationsV2(kFALSE),
  50 fMinV2(0.04),
  51 fMaxV2(0.06),
  52 fPtDependentV2(kFALSE),
  53 fV2vsPtCutOff(2.0),
  54 fV2vsPtMax(0.2),
  55 fEtaMinA(-0.8),
  56 fEtaMaxA(-0.5),
  57 fEtaMinB(0.5),
  58 fEtaMaxB(0.8),
  59 fNTimes(1),
  60 fUniformAcceptance(kTRUE),
  61 fPhiMin1(0.),
  62 fPhiMax1(0.),
  63 fProbability1(0.),
  64 fPhiMin2(0.),
  65 fPhiMax2(0.),
  66 fProbability2(0.),
  67 fPi(TMath::Pi())
  68 {
  69  // Constructor.
  70
  71  // Determine seed for gRandom:
  72  delete gRandom;
  73  gRandom = new TRandom3(uiSeed); // if uiSeed is 0, the seed is determined uniquely in space and time via TUUID
  74
  75 } // end of AliFlowEventSimpleMakerOnTheFly::AliFlowEventSimpleMakerOnTheFly(UInt_t uiSeed):
  76
  77 //====================================================================================================================
  78
  79 AliFlowEventSimpleMakerOnTheFly::~AliFlowEventSimpleMakerOnTheFly()
  80 {
  81  // Destructor.
  82
  83  if(fPtSpectra){delete fPtSpectra;}
  84  if(fPhiDistribution){delete fPhiDistribution;}
  85
  86 } // end of AliFlowEventSimpleMakerOnTheFly::~AliFlowEventSimpleMakerOnTheFly()
  87
  88 //====================================================================================================================
  89
  90 void AliFlowEventSimpleMakerOnTheFly::Init()
  91 {
  92  // Book all objects in this method.
  93
  94  // a) Define the pt spectra;
  95  // b) Define the phi distribution.
  96
  97  // a) Define the pt spectra:
  98  Double_t dPtMin = 0.;
  99  Double_t dPtMax = 10.;
 100  fPtSpectra = new TF1("fPtSpectra","x*TMath::Exp(-pow([0]*[0]+x*x,0.5)/[1])",dPtMin,dPtMax); // hardwired is Boltzmann distribution
 101  fPtSpectra->SetParName(0,"Mass");
 102  fPtSpectra->SetParameter(0,fMass);
 103  fPtSpectra->SetParName(1,"Temperature");
 104  fPtSpectra->SetParameter(1,fTemperature);
 105  fPtSpectra->SetTitle("Boltzmann Distribution: f(p_{t}) = p_{t}exp[-(m^{2}+p_{t}^{2})^{1/2}/T];p_{t};f(p_{t})");
 106
 107  // b) Define the phi distribution:
 108  Double_t dPhiMin = 0.;
 109  Double_t dPhiMax = TMath::TwoPi();
 110  fPhiDistribution = new TF1("fPhiDistribution","1+2.*[1]*TMath::Cos(x-[0])+2.*[2]*TMath::Cos(2.*(x-[0]))+2.*[3]*TMath::Cos(3.*(x-[0]))+2.*[4]*TMath::Cos(4.*(x-[0]))",dPhiMin,dPhiMax);
 111  fPhiDistribution->SetParName(0,"Reaction Plane");
 112  fPhiDistribution->SetParameter(0,0.);
 113  fPhiDistribution->SetParName(1,"Directed Flow (v1)");
 114  fPhiDistribution->SetParameter(1,fV1);
 115  fPhiDistribution->SetParName(2,"Elliptic Flow (v2)");
 116  fPhiDistribution->SetParameter(2,fV2);
 117  fPhiDistribution->SetParName(3,"Triangular Flow (v3)");
 118  fPhiDistribution->SetParameter(3,fV3);
 119  fPhiDistribution->SetParName(4,"Quadrangular Flow (v4)");
 120  fPhiDistribution->SetParameter(4,fV4);
 121
 122 } // end of void AliFlowEventSimpleMakerOnTheFly::Init()
 123
 124 //====================================================================================================================
 125
 126 Bool_t AliFlowEventSimpleMakerOnTheFly::AcceptOrNot(AliFlowTrackSimple *pTrack)
 127 {
 128  // For the case of non-uniform acceptance determine in this method if particle is accepted or rejected.
 129
 130  Bool_t bAccept = kTRUE;
 131
 132  if((pTrack->Phi() >= fPhiMin1*fPi/180.) && (pTrack->Phi() < fPhiMax1*fPi/180.) && gRandom->Uniform(0,1) > fProbability1)
 133  {
 134   bAccept = kFALSE; // particle is rejected in the first non-uniform sector
 135  } else if((pTrack->Phi() >= fPhiMin2*fPi/180.) && (pTrack->Phi() < fPhiMax2*fPi/180.) && gRandom->Uniform(0,1) > fProbability2)
 136     {
 137      bAccept = kFALSE; // particle is rejected in the second non-uniform sector
 138     }
 139
 140  return bAccept;
 141
 142 } // end of Bool_t AliFlowEventSimpleMakerOnTheFly::AcceptOrNot(AliFlowTrackSimple *pTrack);
 143
 144 //====================================================================================================================
 145
 146 AliFlowEventSimple* AliFlowEventSimpleMakerOnTheFly::CreateEventOnTheFly(AliFlowTrackSimpleCuts const *cutsRP, AliFlowTrackSimpleCuts const *cutsPOI)
 147 {
 148  // Method to create event 'on the fly'.
 149
 150  // a) Determine the multiplicity of an event;
 151  // b) Determine the reaction plane of an event;
 152  // c) If v2 fluctuates uniformly event-by-event, sample its value from [fMinV2,fMaxV2];
 153  // d) Create event 'on the fly';
 154  // e) Cosmetics for the printout on the screen.
 155
 156  // a) Determine the multiplicity of an event:
 157  Int_t iMult = (Int_t)gRandom->Uniform(fMinMult,fMaxMult);
 158
 159  // b) Determine the reaction plane of an event:
 160  Double_t dReactionPlane = gRandom->Uniform(0.,TMath::TwoPi());
 161  fPhiDistribution->SetParameter(0,dReactionPlane);
 162
 163  // c) If v2 fluctuates uniformly event-by-event, sample its value from [fMinV2,fMaxV2]:
 164  if(fUniformFluctuationsV2)
 165  {
 166   fPhiDistribution->SetParameter(2,gRandom->Uniform(fMinV2,fMaxV2));
 167  }
 168
 169  // d) Create event 'on the fly':
 170  AliFlowEventSimple *pEvent = new AliFlowEventSimple(iMult);
 171  pEvent->SetReferenceMultiplicity(iMult);
 172  pEvent->SetMCReactionPlaneAngle(dReactionPlane);
 173  Int_t nRPs = 0; // number of particles tagged RP in this event
 174  Int_t nPOIs = 0; // number of particles tagged POI in this event
 175  for(Int_t p=0;p<iMult;p++)
 176  {
 177   AliFlowTrackSimple *pTrack = new AliFlowTrackSimple();
 178   pTrack->SetPt(fPtSpectra->GetRandom());
 179   if(fPtDependentV2 && !fUniformFluctuationsV2)
 180   {
 181    // v2(pt): for pt < fV2vsPtCutOff v2 increases linearly, for pt >= fV2vsPtCutOff v2 = fV2vsPtMax
 182    (pTrack->Pt() < fV2vsPtCutOff ?
 183     fPhiDistribution->SetParameter(2,pTrack->Pt()*fV2vsPtMax/fV2vsPtCutOff) :
 184     fPhiDistribution->SetParameter(2,fV2vsPtMax)
 185    );
 186   } // end of if(fPtDependentV2)
 187   pTrack->SetPhi(fPhiDistribution->GetRandom());
 188   pTrack->SetEta(gRandom->Uniform(-1.,1.));
 189   pTrack->SetCharge((gRandom->Integer(2)>0.5 ? 1 : -1));
 190   // Check uniform acceptance:
 191   if(!fUniformAcceptance && !AcceptOrNot(pTrack)){continue;}
 192   // Checking the RP cuts:
 193   if(cutsRP->PassesCuts(pTrack))
 194   {
 195    pTrack->TagRP(kTRUE);
 196    nRPs++;
 197   }
 198   // Checking the POI cuts:
 199   if(cutsPOI->PassesCuts(pTrack))
 200   {
 201    pTrack->TagPOI(kTRUE);
 202    nPOIs++;
 203   }
 204   // Assign particles to eta subevents (needed only for Scalar Product method):
 205   if(pTrack->Eta()>=fEtaMinA && pTrack->Eta()<fEtaMaxA)
 206   {
 207    pTrack->SetForSubevent(0);
 208   }
 209   if(pTrack->Eta()>=fEtaMinB && pTrack->Eta()<fEtaMaxB)
 210   {
 211    pTrack->SetForSubevent(1);
 212   }
 213   pEvent->AddTrack(pTrack);
 214   // Simulating nonflow:
 215   if(fNTimes>1)
 216   {
 217    for(Int_t nt=1;nt<fNTimes;nt++)
 218    {
 219     pEvent->AddTrack(pTrack->Clone());
 220    }
 221   } // end of if(fNTimes>1)
 222  } // end of for(Int_t p=0;p<iMult;p++)
 223  pEvent->SetNumberOfRPs(fNTimes*nRPs);
 224
 225  // e) Cosmetics for the printout on the screen:
 226  Int_t cycle = (fPtDependentV2 ? 10 : 100);
 227  if((++fCount % cycle) == 0)
 228  {
 229   if(TMath::Abs(dReactionPlane)>1.e-44)
 230   {
 231    cout<<" MC Reaction Plane Angle = "<<dReactionPlane<<endl;
 232   } else
 233     {
 234      cout<<" MC Reaction Plane Angle is unknown :'( "<< endl;
 235     }
 236   cout<<" # of simulated tracks  = "<<fNTimes*iMult<<endl;
 237   cout<<" # of RP tagged tracks  = "<<fNTimes*nRPs<<endl;
 238   cout<<" # of POI tagged tracks = "<<fNTimes*nPOIs<<endl;
 239   cout <<"  .... "<<fCount<< " events processed ...."<<endl;
 240  } // end of if((++fCount % cycle) == 0)
 241
 242  return pEvent;
 243
 244 } // end of CreateEventOnTheFly()
 245
 246 //====================================================================================================================
 247