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