few extra options for momentum dependence v_n
authorsnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 4 May 2009 17:16:01 +0000 (17:16 +0000)
committersnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 4 May 2009 17:16:01 +0000 (17:16 +0000)
PWG2/FLOW/AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.cxx
PWG2/FLOW/AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.h
PWG2/FLOW/macros/runFlowAnalysisOnTheFly.C

index 0a25f7b..ab045fa 100644 (file)
@@ -27,6 +27,7 @@
 #include "TMath.h"
 #include "TF1.h"
 #include "TRandom3.h"
+//#include "TUnuran.h"
 
 #include "AliFlowEventSimpleMakerOnTheFly.h"
 #include "AliFlowEventSimple.h"
@@ -41,10 +42,15 @@ ClassImp(AliFlowEventSimpleMakerOnTheFly)
 AliFlowEventSimpleMakerOnTheFly::AliFlowEventSimpleMakerOnTheFly(UInt_t iseed):
   fMultiplicityOfRP(0),
   fMultiplicitySpreadOfRP(0.),
+  fUseConstantHarmonics(kFALSE),
   fV1RP(0.), 
   fV1SpreadRP(0.), 
   fV2RP(0.), 
   fV2SpreadRP(0.), 
+  fV4RP(0.), 
+  fV4SpreadRP(0.), 
+  fV2RPMax(0.), 
+  fPtCutOff(0.), 
   fPtSpectra(NULL),
   fPhiDistribution(NULL),
   fMyTRandom3(NULL),
@@ -52,7 +58,11 @@ AliFlowEventSimpleMakerOnTheFly::AliFlowEventSimpleMakerOnTheFly(UInt_t iseed):
   fNoOfLoops(1)
  {
   // constructor
-   fMyTRandom3 = new TRandom3(iseed); 
+  fMyTRandom3 = new TRandom3(iseed);   
+  gRandom->SetSeed(fMyTRandom3->Integer(65539));
+  //fMyUnuran = new TUnuran(); 
+  //fMyUnuran->SetSeed(iseed);  
  }
 
 
@@ -62,54 +72,71 @@ AliFlowEventSimpleMakerOnTheFly::AliFlowEventSimpleMakerOnTheFly(UInt_t iseed):
 AliFlowEventSimpleMakerOnTheFly::~AliFlowEventSimpleMakerOnTheFly()
 {
  // destructor
-  if (fPtSpectra) delete fPtSpectra;
-  if (fPhiDistribution) delete fPhiDistribution;
-  if (fMyTRandom3) delete  fMyTRandom3;
+ if (fPtSpectra) delete fPtSpectra;
+ if (fPhiDistribution) delete fPhiDistribution;
+ if (fMyTRandom3) delete fMyTRandom3;
 }
 
 
 //========================================================================
 
 
+void AliFlowEventSimpleMakerOnTheFly::Init()
+{
+ // define the pt spectra and phi distribution
+ // pt spectra:   
+ Double_t dPtMin = 0.; // to be improved (move this to the body of contstructor?)
+ Double_t dPtMax = 10.; // to be improved (move this to the body of contstructor?) 
+  
+ fPtSpectra = new TF1("fPtSpectra","[0]*x*TMath::Exp(-x*x)",dPtMin,dPtMax);  
+ fPtSpectra->SetParName(0,"Multiplicity of RPs");  
+ // phi distribution:
+ Double_t dPhiMin = 0.; // to be improved (move this to the body of contstructor?)
+ Double_t dPhiMax = TMath::TwoPi(); // to be improved (move this to the body of contstructor?)
+  
+ 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);
+ fPhiDistribution->SetParName(0,"directed flow");
+ fPhiDistribution->SetParName(1,"elliptic flow"); 
+ fPhiDistribution->SetParName(2,"Reaction Plane");
+ fPhiDistribution->SetParName(3,"harmonic 4"); // to be improved (name)
+}
+
+//========================================================================
+
 AliFlowEventSimple* AliFlowEventSimpleMakerOnTheFly::CreateEventOnTheFly()
 {
   // method to create event on the fly
   
   AliFlowEventSimple* pEvent = new AliFlowEventSimple(fMultiplicityOfRP);
-    
-  // pt:   
-  Double_t dPtMin = 0.; // to be improved 
-  Double_t dPtMax = 10.; // to be improved 
-  
-  fPtSpectra = new TF1("fPtSpectra","[0]*x*TMath::Exp(-x*x)",dPtMin,dPtMax);  
-  fPtSpectra->SetParName(0,"Multiplicity of RPs");  
+   
   // sampling the multiplicity:
-  Int_t fNewMultiplicityOfRP = fMyTRandom3->Gaus(fMultiplicityOfRP,fMultiplicitySpreadOfRP);
-  fPtSpectra->SetParameter(0,fNewMultiplicityOfRP);
+  Int_t iNewMultiplicityOfRP = fMultiplicityOfRP;
+  if(fMultiplicitySpreadOfRP>0.0) iNewMultiplicityOfRP = (Int_t)fMyTRandom3->Gaus(fMultiplicityOfRP,fMultiplicitySpreadOfRP);
+  fPtSpectra->SetParameter(0,iNewMultiplicityOfRP);
   
-  
-  // phi:
-  Double_t dPhiMin = 0.; // to be improved 
-  Double_t dPhiMax = TMath::TwoPi(); // to be improved 
-  
-  fPhiDistribution = new TF1("fPhiDistribution","1+2.*[0]*TMath::Cos(x-[2])+2.*[1]*TMath::Cos(2*(x-[2]))",dPhiMin,dPhiMax);
-  
-  // samling the reaction plane
-  Double_t fMCReactionPlaneAngle = fMyTRandom3->Uniform(0.,TMath::TwoPi());
-  fPhiDistribution->SetParName(2,"Reaction Plane");
-  fPhiDistribution->SetParameter(2,fMCReactionPlaneAngle);
+  // sampling the reaction plane
+  Double_t dMCReactionPlaneAngle = fMyTRandom3->Uniform(0.,TMath::TwoPi());
+  fPhiDistribution->SetParameter(2,dMCReactionPlaneAngle);
 
   // sampling the V1:
-  fPhiDistribution->SetParName(0,"directed flow");
-  Double_t fNewV1RP=fV1RP;
-  if(fV1SpreadRP>0.0) {fNewV1RP = fMyTRandom3->Gaus(fV1RP,fV1SpreadRP);}
-  fPhiDistribution->SetParameter(0,fNewV1RP);
+  Double_t dNewV1RP=fV1RP;
+  if(fV1SpreadRP>0.0) {dNewV1RP = fMyTRandom3->Gaus(fV1RP,fV1SpreadRP);}
+  fPhiDistribution->SetParameter(0,dNewV1RP);
  
   // sampling the V2:
-  fPhiDistribution->SetParName(1,"elliptic flow");
-  Double_t fNewV2RP = fV2RP;
-  if(fV2SpreadRP>0.0) fNewV2RP = fMyTRandom3->Gaus(fV2RP,fV2SpreadRP);
-  fPhiDistribution->SetParameter(1,fNewV2RP);
+  if(fUseConstantHarmonics)
+  {
+   Double_t dNewV2RP = fV2RP;
+   if(fV2SpreadRP>0.0) dNewV2RP = fMyTRandom3->Gaus(fV2RP,fV2SpreadRP);
+   fPhiDistribution->SetParameter(1,dNewV2RP);
+  }
+  
+  // sampling the V4:
+  Double_t dNewV4RP = fV4RP;
+  if(fV4SpreadRP>0.0) dNewV4RP = fMyTRandom3->Gaus(fV4RP,fV4SpreadRP);
+  fPhiDistribution->SetParameter(3,dNewV4RP);
    
   // eta:
   Double_t dEtaMin = -1.; // to be improved 
@@ -118,19 +145,29 @@ AliFlowEventSimple* AliFlowEventSimpleMakerOnTheFly::CreateEventOnTheFly()
   Int_t iGoodTracks = 0;
   Int_t iSelParticlesRP = 0;
   Int_t iSelParticlesPOI = 0;
-  Double_t fTmpPt =0;
-  Double_t fTmpEta =0;
-  Double_t fTmpPhi =0;
-  
-  for(Int_t i=0;i<fNewMultiplicityOfRP;i++) {
-    fTmpPt = fPtSpectra->GetRandom();
-    fTmpEta = fMyTRandom3->Uniform(dEtaMin,dEtaMax);
-    fTmpPhi = fPhiDistribution->GetRandom();
+  Double_t dTmpPt = 0.;
+  Double_t dTmpEta = 0.;
+  Double_t dTmpV2 = 0.;
+  Double_t dTmpPhi = 0.;
+  for(Int_t i=0;i<iNewMultiplicityOfRP;i++) {
+    dTmpPt = fPtSpectra->GetRandom();
+    dTmpEta = fMyTRandom3->Uniform(dEtaMin,dEtaMax);
+    // to be improved:
+    if(!fUseConstantHarmonics) {
+      if(dTmpPt >= fPtCutOff) { 
+       dTmpV2 = fV2RPMax;
+      } else {
+       dTmpV2 = fV2RPMax*(dTmpPt/fPtCutOff);
+      }  
+      fPhiDistribution->SetParameter(1,dTmpV2);         
+    }
+    dTmpPhi = fPhiDistribution->GetRandom();
+    // add the track to the event
     for(Int_t d=0;d<fNoOfLoops;d++) {
       AliFlowTrackSimple* pTrack = new AliFlowTrackSimple();
-      pTrack->SetPt(fTmpPt);
-      pTrack->SetEta(fTmpEta);
-      pTrack->SetPhi(fTmpPhi);
+      pTrack->SetPt(dTmpPt);
+      pTrack->SetEta(dTmpEta);
+      pTrack->SetPhi(dTmpPhi);
       pTrack->SetForRPSelection(kTRUE);
       iSelParticlesRP++;
       pTrack->SetForPOISelection(kTRUE);
@@ -139,21 +176,20 @@ AliFlowEventSimple* AliFlowEventSimpleMakerOnTheFly::CreateEventOnTheFly()
       iGoodTracks++;
     }
   }
+  
+  // update the event quantities
   pEvent->SetEventNSelTracksRP(iSelParticlesRP);  
   pEvent->SetNumberOfTracks(iGoodTracks);//tracks used either for RP or for POI selection
-  pEvent->SetMCReactionPlaneAngle(fMCReactionPlaneAngle);
+  pEvent->SetMCReactionPlaneAngle(dMCReactionPlaneAngle);
 
-  if (!fMCReactionPlaneAngle == 0) cout<<" MC Reaction Plane Angle = "<<  fMCReactionPlaneAngle << endl;
+  if (!dMCReactionPlaneAngle == 0) cout<<" MC Reaction Plane Angle = "<<  dMCReactionPlaneAngle << endl;
   else cout<<" MC Reaction Plane Angle = unknown "<< endl;
 
   cout<<" iGoodTracks = "<< iGoodTracks << endl;
   cout<<" # of RP selected tracks = "<<iSelParticlesRP<<endl;
   cout<<" # of POI selected tracks = "<<iSelParticlesPOI<<endl;  
   cout << "# " << ++fCount << " events processed" << endl;
-
-  delete fPhiDistribution;
-  delete fPtSpectra;
+  
   return pEvent;  
  
 } // end of CreateEventOnTheFly()
index 07209e6..99f0fed 100644 (file)
@@ -19,6 +19,7 @@
 
 class TF1;
 class TRandom3;
+//class TUnuran;
 
 #include "AliFlowEventSimple.h"  //needed as include
     
@@ -29,6 +30,8 @@ class AliFlowEventSimpleMakerOnTheFly {
   AliFlowEventSimpleMakerOnTheFly(UInt_t);    // constructor
   virtual ~AliFlowEventSimpleMakerOnTheFly(); // destructor
 
+  virtual void Init(); 
+  
   AliFlowEventSimple* CreateEventOnTheFly();  // create an event on the fly
  
     
@@ -42,7 +45,11 @@ class AliFlowEventSimpleMakerOnTheFly {
   
   void SetMultiplicitySpreadOfRP(Double_t multSpreadRP) {this->fMultiplicitySpreadOfRP = multSpreadRP;}
   Double_t GetMultiplicitySpreadOfRP() const {return this->fMultiplicitySpreadOfRP;} 
-
+  
+  void SetUseConstantHarmonics(Bool_t const uch) {this->fUseConstantHarmonics = uch;};
+  Bool_t GetUseConstantHarmonics() const {return this->fUseConstantHarmonics;};
+  
+  // constant harmonics:  
   void SetV1RP(Double_t dV1RP) {this->fV1RP = dV1RP;}
   Double_t GetV1RP() const {return this->fV1RP;} 
   
@@ -54,6 +61,19 @@ class AliFlowEventSimpleMakerOnTheFly {
   
   void SetV2SpreadRP(Double_t dV2SpreadRP) {this->fV2SpreadRP = dV2SpreadRP;}
   Double_t GetV2SpreadRP() const {return this->fV2SpreadRP;} 
+  
+  void SetV4RP(Double_t dV4RP) {this->fV4RP = dV4RP;}
+  Double_t GetV4RP() const {return this->fV4RP;} 
+  
+  void SetV4SpreadRP(Double_t dV4SpreadRP) {this->fV4SpreadRP = dV4SpreadRP;}
+  Double_t GetV4SpreadRP() const {return this->fV4SpreadRP;} 
+  
+  // (pt,eta) dependent harmonics:
+  void SetV2RPMax(Double_t dV2RPMax) {this->fV2RPMax = dV2RPMax;}
+  Double_t GetV2RPMax() const {return this->fV2RPMax;} 
+  
+  void SetPtCutOff(Double_t dPtCutOff) {this->fPtCutOff = dPtCutOff;}
+  Double_t GetPtCutOff() const {return this->fPtCutOff;} 
   //................................................................................................
   
   void SetNoOfLoops(Int_t noofl) {this->fNoOfLoops = noofl;}
@@ -66,12 +86,19 @@ class AliFlowEventSimpleMakerOnTheFly {
   
   //................................................................................................
   // global parameters:
-  Int_t     fMultiplicityOfRP;          // multiplicity of RPs
-  Double_t  fMultiplicitySpreadOfRP; // multiplicity spread of RPs  
+  Int_t     fMultiplicityOfRP;       // multiplicity of RPs
+  Double_t  fMultiplicitySpreadOfRP; // multiplicity spread of RPs 
+  Bool_t    fUseConstantHarmonics;      // harmonics V1, V2, V4... are constant (kTRUE) or functions of pt and eta (kFALSE)     
+  // constant harmonics: 
   Double_t  fV1RP;                   // directed flow of RPs
   Double_t  fV1SpreadRP;             // directed flow spread of RPs
   Double_t  fV2RP;                   // elliptic flow of RPs
   Double_t  fV2SpreadRP;             // elliptic flow spread of RPs
+  Double_t  fV4RP;                   // harmonic V4 of RPs
+  Double_t  fV4SpreadRP;             // harmonic V4's spread of RPs
+  // (pt,eta) dependent harmonics:
+  Double_t  fV2RPMax;                // elliptic flow of RPs
+  Double_t  fPtCutOff;               // elliptic flow spread of RPs
   //................................................................................................
   
   //................................................................................................
@@ -80,10 +107,10 @@ class AliFlowEventSimpleMakerOnTheFly {
   TF1*      fPhiDistribution; // azimuthal distribution
   //................................................................................................
   
-  TRandom3* fMyTRandom3; // our random generator
-  Int_t     fCount;
-  Int_t     fNoOfLoops; // number of times to use the same particle (nonflow)
-
+  TRandom3* fMyTRandom3; // our TRandom3 generator
+  //TUnuran*  fMyUnuran;   // our TUnuran generator
+  Int_t     fCount;      // count number of events 
+  Int_t     fNoOfLoops;  // number of times to use the same particle (nonflow)
 
   ClassDef(AliFlowEventSimpleMakerOnTheFly,0) // macro for rootcint
 };
index 80241f7..3c7c065 100644 (file)
@@ -16,20 +16,34 @@ Bool_t FQD   = kTRUE;
 Bool_t MCEP  = kTRUE; 
 //--------------------------------------------------------------------------------------
 
-//......................................................................................  
+Bool_t bSameSeed = kFALSE; // use always the same seed for random generators
+Bool_t bConstantHarmonics = kTRUE; // harmonics V1, V2, V4... are constant (kTRUE) or functions of pt and eta (kFALSE)
+
 // Set the event parameters:
 Int_t iLoops = 1; // number of times to use each track (to simulate nonflow)
+
 Int_t iMultiplicityOfRP = 500; // multiplicity of RPs
 Double_t dMultiplicitySpreadOfRP = 0; // multiplicity spread of RPs
 
+//......................................................................................  
+// if you use (pt,eta) dependent harmonics (bConstantHarmonics = kFALSE):
+Double_t dPtCutOff = 2.0; // V2(pt) is linear up to pt = 2 GeV and for pt > 2 GeV it is constant: V2(pt) = dVRPMax
+Double_t dV2RPMax = 0.2; // maximum value of V2(pt) for pt >= 2GeV
+//...................................................................................... 
+
+//......................................................................................  
+// if you use constant harmonics (bConstantHarmonics = kTRUE):
 Double_t dV2RP = 0.05; // elliptic flow of RPs
 Double_t dV2SpreadRP = 0.; // elliptic flow spread of RPs
 
 Double_t dV1RP = 0.0; // directed flow of RPs
 Double_t dV1SpreadRP = 0.0; // directed flow spread of RPs
-//...................................................................................... 
 
-enum anaModes {mLocal,mLocalSource,mLocalPAR,};
+Double_t dV4RP = 0.0; // harmonic V4 of RPs (to be improved: name needed)
+Double_t dV4SpreadRP = 0.0; // harmonic V4's spread of RPs (to be improved: name needed)
+//......................................................................................  
+
+enum anaModes {mLocal,mLocalSource,mLocalPAR};
 // mLocal: Analyze data on your computer using aliroot
 // mLocalPAR: Analyze data on your computer using root + PAR files
 // mLocalSource: Analyze data on your computer using root + source files
@@ -55,13 +69,24 @@ int runFlowAnalysisOnTheFly(Int_t mode=mLocal, Int_t nEvts=1000)
  
  LoadLibraries(mode);
 
- // Initialize the random generator
- TTimeStamp dt;
- UInt_t sseed = dt.GetNanoSec()/1000;
+ // Initialize the seed for random generator
+ UInt_t sseed = 0;
+ if(bSameSeed) 
+ {
+  sseed = 44; // the default constant value for seed for random generators
+ } 
 
+ if(!bSameSeed)
+ {
+  TTimeStamp dt;
+  sseed = dt.GetNanoSec()/1000;
+ }
  //---------------------------------------------------------------------------------------
  // Initialize the flowevent maker
  AliFlowEventSimpleMakerOnTheFly* eventMakerOnTheFly = new AliFlowEventSimpleMakerOnTheFly(sseed);
+ eventMakerOnTheFly->Init();
   
  //---------------------------------------------------------------------------------------
  // Initialize all the flow methods:  
@@ -164,11 +189,27 @@ int runFlowAnalysisOnTheFly(Int_t mode=mLocal, Int_t nEvts=1000)
  eventMakerOnTheFly->SetNoOfLoops(iLoops);
  eventMakerOnTheFly->SetMultiplicityOfRP(iMultiplicityOfRP);
  eventMakerOnTheFly->SetMultiplicitySpreadOfRP(dMultiplicitySpreadOfRP);
+
  eventMakerOnTheFly->SetV1RP(dV1RP);
  eventMakerOnTheFly->SetV1SpreadRP(dV1SpreadRP);  
- eventMakerOnTheFly->SetV2RP(dV2RP);
- eventMakerOnTheFly->SetV2SpreadRP(dV2SpreadRP);  
-      
+ eventMakerOnTheFly->SetV4RP(dV4RP);
+ eventMakerOnTheFly->SetV4SpreadRP(dV4SpreadRP);  
+ // constant harmonic V2:
+ if(bConstantHarmonics)
+ { 
+  eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
+  eventMakerOnTheFly->SetV2RP(dV2RP);
+  eventMakerOnTheFly->SetV2SpreadRP(dV2SpreadRP);  
+ }
+ // (pt,eta) dependent harmonic V2:
+ if(!bConstantHarmonics)
+ {
+  eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
+  eventMakerOnTheFly->SetV2RPMax(dV2RPMax);
+  eventMakerOnTheFly->SetPtCutOff(dPtCutOff);  
+ }
+       
  //---------------------------------------------------------------------------------------  
  // create and analyze events 'on the fly':