preparations for glauber inclusion
authorsnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 4 May 2010 15:14:26 +0000 (15:14 +0000)
committersnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 4 May 2010 15:14:26 +0000 (15:14 +0000)
PWG2/FLOW/AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.cxx
PWG2/FLOW/AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.h
PWG2/FLOW/macros/runFlowAnalysisOnTheFly.C

index 6d15588..79d2cdc 100644 (file)
@@ -39,6 +39,7 @@ ClassImp(AliFlowEventSimpleMakerOnTheFly)
 
 
 AliFlowEventSimpleMakerOnTheFly::AliFlowEventSimpleMakerOnTheFly(UInt_t iseed):
+  fUseGlauberModel(kFALSE),
   fMultDistrOfRPsIsGauss(kFALSE),
   fMultiplicityOfRP(0),
   fMultiplicitySpreadOfRP(0.),
@@ -133,16 +134,13 @@ void AliFlowEventSimpleMakerOnTheFly::Init()
 
 //========================================================================
 
-AliFlowEventSimple* AliFlowEventSimpleMakerOnTheFly::CreateEventOnTheFly(AliFlowTrackSimpleCuts *cutsRP, AliFlowTrackSimpleCuts *cutsPOI)
+Int_t AliFlowEventSimpleMakerOnTheFly::DetermineMultiplicity()
 {
-  // method to create event on the fly
-  
-  AliFlowEventSimple* pEvent = new AliFlowEventSimple(fMultiplicityOfRP);
-   
-  // sampling the multiplicity:
-  Int_t iNewMultiplicityOfRP = fMultiplicityOfRP;
+ // Determine multiplicity for current event.
+ Int_t iNewMultiplicityOfRP = fMultiplicityOfRP;
   
-  if(fMultDistrOfRPsIsGauss) {
+ if(fMultDistrOfRPsIsGauss) {
     if (fMultiplicitySpreadOfRP>0.0) iNewMultiplicityOfRP = (Int_t)fMyTRandom3->Gaus(fMultiplicityOfRP,fMultiplicitySpreadOfRP);
     fPtSpectra->SetParameter(0,iNewMultiplicityOfRP);
   } else {
@@ -154,50 +152,116 @@ AliFlowEventSimple* AliFlowEventSimpleMakerOnTheFly::CreateEventOnTheFly(AliFlow
     }
   }
   
-  // set the 'temperature' of RPs
-  fPtSpectra->SetParameter(1,fTemperatureOfRP);  
+ return iNewMultiplicityOfRP;
+} // end of AliFlowEventSimpleMakerOnTheFly::DetermineMultiplicity()
+
+//========================================================================
+
+void AliFlowEventSimpleMakerOnTheFly::DetermineV1()
+{
+ // Determine flow harmonics v1 for current event (if v1 is not pt or eta dependent).
+
+ Double_t dNewV1RP=fV1RP;
+ if(fV1SpreadRP>0.0) {dNewV1RP = fMyTRandom3->Gaus(fV1RP,fV1SpreadRP);}
+ fPhiDistribution->SetParameter(0,dNewV1RP);
   
-  // sampling the reaction plane
-  Double_t dMCReactionPlaneAngle = fMyTRandom3->Uniform(0.,TMath::TwoPi());
-  fPhiDistribution->SetParameter(2,dMCReactionPlaneAngle);
+} // end of void AliFlowEventSimpleMakerOnTheFly::DetermineV1()
   
-  // sampling the V1:
-  if(!(fPtDependentHarmonicV1 || fEtaDependentHarmonicV1))
-  {
-   Double_t dNewV1RP=fV1RP;
-   if(fV1SpreadRP>0.0) {dNewV1RP = fMyTRandom3->Gaus(fV1RP,fV1SpreadRP);}
-   fPhiDistribution->SetParameter(0,dNewV1RP);
-  }
-  
-  // sampling the V2:
-  if(!(fPtDependentHarmonicV2 || fEtaDependentHarmonicV2)) 
+//========================================================================
+
+void AliFlowEventSimpleMakerOnTheFly::DetermineV4()
+{
+ // Determine flow harmonics v4 for current event (if v4 is not pt or eta dependent).
+ Double_t dNewV4RP = fV4RP;
+ if(fV4SpreadRP>0.0) dNewV4RP = fMyTRandom3->Gaus(fV4RP,fV4SpreadRP);
+ fPhiDistribution->SetParameter(3,dNewV4RP);
+
+} // end of void AliFlowEventSimpleMakerOnTheFly::DetermineV4()
+
+//========================================================================
+
+void AliFlowEventSimpleMakerOnTheFly::DetermineV2()
+{
+ // Determine flow harmonics v2 for current event (if v2 is not pt or eta dependent).
+
+ Double_t dNewV2RP = fV2RP;
+ if(fConstantV2IsSampledFromGauss) 
+ {
+  if(fV2SpreadRP>0.0) 
   {
-   Double_t dNewV2RP = fV2RP;
-   if(fConstantV2IsSampledFromGauss) 
+   dNewV2RP = fMyTRandom3->Gaus(fV2RP,fV2SpreadRP);}
+  fPhiDistribution->SetParameter(1,dNewV2RP);
+ } else if(fMinV2RP < fMaxV2RP) 
    {
-    if(fV2SpreadRP>0.0) 
-    {
-     dNewV2RP = fMyTRandom3->Gaus(fV2RP,fV2SpreadRP);
-    }
-    fPhiDistribution->SetParameter(1,dNewV2RP);
-   } else if(fMinV2RP < fMaxV2RP) 
+    dNewV2RP = fMyTRandom3->Uniform(fMinV2RP,fMaxV2RP);   
+    fPhiDistribution->SetParameter(1,dNewV2RP);  
+   } else if(fMinV2RP == fMaxV2RP)
      {
-      dNewV2RP = fMyTRandom3->Uniform(fMinV2RP,fMaxV2RP);    
-          fPhiDistribution->SetParameter(1,dNewV2RP);  
-     } else if(fMinV2RP == fMaxV2RP)
-        {
-             dNewV2RP = fMinV2RP;
-             fPhiDistribution->SetParameter(1,dNewV2RP);          
-        }
-  } // end of if(!(bPtDependentHarmonicV2 || bEtaDependentHarmonicV2))
+      dNewV2RP = fMinV2RP;
+      fPhiDistribution->SetParameter(1,dNewV2RP);          
+     }
+
+} // end of void AliFlowEventSimpleMakerOnTheFly::DetermineV2()
+
+//========================================================================
+
+Int_t AliFlowEventSimpleMakerOnTheFly::GlauberModel()
+{
+ // Determine multiplicity and flow harmonics for current event from Glauber moder
+ Int_t multiplicity = 0;
+ Double_t v1 = 0.;
+ Double_t v2 = 0.;
+ Double_t v4 = 0.;
+ // Determine multiplicity, v1, v2 and v4 from Glauber model:
+  
+ // multiplicity = ... 
+ // v1 = ...
+ // v2 = ...
+ // v4 = ...
+ // Set obtained values as parameters in relevant distributions:
+ fPtSpectra->SetParameter(0,multiplicity);
+ fPhiDistribution->SetParameter(0,v1);
+ fPhiDistribution->SetParameter(1,v2);
+ fPhiDistribution->SetParameter(3,v4);
+ return multiplicity;
+} // end of Int_t AliFlowEventSimpleMakerOnTheFly::GlauberModel()
+
+//========================================================================
+
+AliFlowEventSimple* AliFlowEventSimpleMakerOnTheFly::CreateEventOnTheFly(AliFlowTrackSimpleCuts *cutsRP, AliFlowTrackSimpleCuts *cutsPOI)
+{
+  // Method to create event on the fly.
   
-  // sampling the V4:
-  if(!(fPtDependentHarmonicV4 || fEtaDependentHarmonicV4))
+  // Determine multiplicity and flow harmonics (if not pt or eta dependent) for current event:
+  Int_t multiplicityRP = 0;
+  if(!fUseGlauberModel)
   {
-   Double_t dNewV4RP = fV4RP;
-   if(fV4SpreadRP>0.0) dNewV4RP = fMyTRandom3->Gaus(fV4RP,fV4SpreadRP);
-   fPhiDistribution->SetParameter(3,dNewV4RP);
-  }
+   multiplicityRP = DetermineMultiplicity(); 
+   if(!(fPtDependentHarmonicV1||fEtaDependentHarmonicV1)) {DetermineV1();}
+   if(!(fPtDependentHarmonicV2||fEtaDependentHarmonicV2)) {DetermineV2();}
+   if(!(fPtDependentHarmonicV4||fEtaDependentHarmonicV4)) {DetermineV4();}
+  } else
+    {
+     // Determine multipliciy and flow harmonics from Glauber model:
+     multiplicityRP = GlauberModel();
+    }  
+  
+  AliFlowEventSimple* pEvent = new AliFlowEventSimple(multiplicityRP);
+    
+  // set the 'temperature' of RPs
+  fPtSpectra->SetParameter(1,fTemperatureOfRP);  
+  
+  // sampling the reaction plane
+  Double_t dMCReactionPlaneAngle = fMyTRandom3->Uniform(0.,TMath::TwoPi());
+  fPhiDistribution->SetParameter(2,dMCReactionPlaneAngle);
   
   // eta:
   Double_t dEtaMin = -1.; // to be improved 
@@ -225,7 +289,7 @@ AliFlowEventSimple* AliFlowEventSimpleMakerOnTheFly::CreateEventOnTheFly(AliFlow
     bUniformAcceptance = kFALSE;
   }
   // loop over original tracks:
-  for(Int_t i=0;i<iNewMultiplicityOfRP;i++) 
+  for(Int_t i=0;i<multiplicityRP;i++) 
   {
    // sample the pt and eta for original track: 
    dPtOriginalTrack = fPtSpectra->GetRandom();
index c8c3655..ecd4dec 100644 (file)
@@ -32,6 +32,11 @@ class AliFlowEventSimpleMakerOnTheFly {
 
   virtual void Init(); 
   
+  Int_t DetermineMultiplicity(); // determine multiplicity for current event
+  virtual void DetermineV1(); // determine flow harmonics v1 for current event (if v1 is not pt or eta dependent)
+  virtual void DetermineV2(); // determine flow harmonics v2 for current event (if v2 is not pt or eta dependent)
+  virtual void DetermineV4(); // determine flow harmonics v4 for current event (if v4 is not pt or eta dependent)
+  Int_t GlauberModel(); // determine multiplicity and flow harmonics for current event from Glauber moder
   AliFlowEventSimple* CreateEventOnTheFly(AliFlowTrackSimpleCuts *cutsRP, AliFlowTrackSimpleCuts *cutsPOI);  // create an event on the fly
  
     
@@ -40,6 +45,9 @@ class AliFlowEventSimpleMakerOnTheFly {
   //                        *****************************
   //................................................................................................
   // setters and getters for global parameters:
+  void SetUseGlauberModel(Bool_t const ugm) {this->fUseGlauberModel = ugm;};
+  Bool_t GetUseGlauberModel() const {return this->fUseGlauberModel;};
+  
   void SetMultDistrOfRPsIsGauss(Bool_t const mdorig) {this->fMultDistrOfRPsIsGauss = mdorig;};
   Bool_t GetMultDistrOfRPsIsGauss() const {return this->fMultDistrOfRPsIsGauss;};
   
@@ -160,6 +168,7 @@ class AliFlowEventSimpleMakerOnTheFly {
   
   //................................................................................................
   // global parameters:
+  Bool_t    fUseGlauberModel;        // if kTRUE multiplicity and flow harmonics are determined e-b-e from Glauber model
   Bool_t    fMultDistrOfRPsIsGauss;  // 1.) if kTRUE  = multiplicitiy of RPs is sampled e-b-e from Gaussian distribution with
                                      //                 mean = fMultiplicityOfRP and spread = fMultiplicitySpreadOfRP
                                      // 2.) if kFALSE = multiplicitiy of RPs is sampled e-b-e uniformly from 
index ca49710..fb920b0 100644 (file)
@@ -83,6 +83,15 @@ Double_t etaRange = 0.0; // If the original track with pseudorapidity eta is spl
 Double_t nonflowSectorMin = 0.0; // detector's sector in which tracks will be splitted starts at this angle (in degrees)                         
 Double_t nonflowSectorMax = 360.0; // detector's sector in which tracks will be splitted ends at this angle (in degrees)                        
 
+//===GLAUBER MODEL================================================================
+Bool_t bUseGlauberModel = kFALSE; // 1.) if kTRUE = multiplicity and constant flow harmonics are 
+                                  //                determined event-by-event from Glauber model
+                                  //                (pt and/or eta dependence of flow harmonics not supported;
+                                  //                settings for multiplicity and flow harmonics bellow are irrelevant) 
+                                  // 2.) if kFALSE = multiplicity and flow harmonics are determined
+                                  //                 independently event-by-event with bellow settings
+                                  //                 (pt and/or eta dependence of flow harmonics supported)
+
 //===FLOW HARMONICS===============================================================
 // harmonics V1, V2, V4... are constants or functions of pt and eta:         
 Bool_t bPtDependentHarmonicV1 = kFALSE; 
@@ -135,7 +144,7 @@ Bool_t bMultDistrOfRPsIsGauss = kTRUE;
                     
 Int_t iMultiplicityOfRP = 500;        // mean multiplicity of RPs (if sampled from Gaussian)
 Double_t dMultiplicitySpreadOfRP = 0; // multiplicity spread of RPs (if sampled from Gaussian)
-Int_t iMinMultOfRP = 400;              // minimal multiplicity of RPs (if sampled uniformly)
+Int_t iMinMultOfRP = 400;              // minimal multiplicity of RPs(if sampled uniformly)
 Int_t iMaxMultOfRP = 600;             // maximal multiplicity of RPs (if sampled uniformly)
                     
 //===DETECTOR ACCEPTANCE===============================================================
@@ -410,17 +419,22 @@ int runFlowAnalysisOnTheFly(Int_t nEvts=1000, Int_t mode=mLocal)
  eventMakerOnTheFly->SetEtaRange(etaRange);
  eventMakerOnTheFly->SetNonflowSectorMin(nonflowSectorMin*TMath::Pi()/180.);
  eventMakerOnTheFly->SetNonflowSectorMax(nonflowSectorMax*TMath::Pi()/180.);
- if(bMultDistrOfRPsIsGauss)
+ eventMakerOnTheFly->SetUseGlauberModel(bUseGlauberModel);
+ if(!bUseGlauberModel)
  {
-  eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
-  eventMakerOnTheFly->SetMultiplicityOfRP(iMultiplicityOfRP);
-  eventMakerOnTheFly->SetMultiplicitySpreadOfRP(dMultiplicitySpreadOfRP);
- } else
-   {
-    eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
-    eventMakerOnTheFly->SetMinMultOfRP(iMinMultOfRP);
-    eventMakerOnTheFly->SetMaxMultOfRP(iMaxMultOfRP); 
-   }
+  if(bMultDistrOfRPsIsGauss)
+  {
+   eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
+   eventMakerOnTheFly->SetMultiplicityOfRP(iMultiplicityOfRP);
+   eventMakerOnTheFly->SetMultiplicitySpreadOfRP(dMultiplicitySpreadOfRP);
+  } else
+    {
+     eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
+     eventMakerOnTheFly->SetMinMultOfRP(iMinMultOfRP);
+     eventMakerOnTheFly->SetMaxMultOfRP(iMaxMultOfRP); 
+    }
+ } // end of if(!bUseGlauberModel)
   
  eventMakerOnTheFly->SetTemperatureOfRP(dTemperatureOfRP);
  eventMakerOnTheFly->SetPtDependentHarmonicV1(bPtDependentHarmonicV1);
@@ -432,8 +446,11 @@ int runFlowAnalysisOnTheFly(Int_t nEvts=1000, Int_t mode=mLocal)
  // V1:
  if(!(bPtDependentHarmonicV1 || bEtaDependentHarmonicV1))
  {
-  eventMakerOnTheFly->SetV1RP(dV1RP);
-  eventMakerOnTheFly->SetV1SpreadRP(dV1SpreadRP);  
+  if(!bUseGlauberModel)
+  {
+   eventMakerOnTheFly->SetV1RP(dV1RP);
+   eventMakerOnTheFly->SetV1SpreadRP(dV1SpreadRP);  
+  } 
  } else // (pt,eta) dependent V1
    {
     eventMakerOnTheFly->SetV1vsPtEtaMax(dV1vsPtEtaMax);
@@ -445,16 +462,19 @@ int runFlowAnalysisOnTheFly(Int_t nEvts=1000, Int_t mode=mLocal)
  // V2:
  if(!(bPtDependentHarmonicV2 || bEtaDependentHarmonicV2)) // constant V2
  { 
-  eventMakerOnTheFly->SetConstantV2IsSampledFromGauss(bConstantV2IsSampledFromGauss);
-  if(bConstantV2IsSampledFromGauss) // Gauss
+  if(!bUseGlauberModel)
   {
-   eventMakerOnTheFly->SetV2RP(dV2RP);
-   eventMakerOnTheFly->SetV2SpreadRP(dV2SpreadRP);  
-  } else // uniform
-    {
-     eventMakerOnTheFly->SetMinV2RP(dMinV2RP);
-     eventMakerOnTheFly->SetMaxV2RP(dMaxV2RP);  
-    }
+   eventMakerOnTheFly->SetConstantV2IsSampledFromGauss(bConstantV2IsSampledFromGauss);
+   if(bConstantV2IsSampledFromGauss) // Gauss
+   {
+    eventMakerOnTheFly->SetV2RP(dV2RP);
+    eventMakerOnTheFly->SetV2SpreadRP(dV2SpreadRP);  
+   } else // uniform
+     {
+      eventMakerOnTheFly->SetMinV2RP(dMinV2RP);
+      eventMakerOnTheFly->SetMaxV2RP(dMaxV2RP);  
+     }
+  } // end of if(!bUseGlauberModel)  
  } else // (pt,eta) dependent V2
    {
     eventMakerOnTheFly->SetV2vsPtEtaMax(dV2vsPtEtaMax);
@@ -470,8 +490,11 @@ int runFlowAnalysisOnTheFly(Int_t nEvts=1000, Int_t mode=mLocal)
  // V4:
  if(!(bPtDependentHarmonicV4 || bEtaDependentHarmonicV4))
  {
-  eventMakerOnTheFly->SetV4RP(dV4RP);
-  eventMakerOnTheFly->SetV4SpreadRP(dV4SpreadRP);  
+  if(!bUseGlauberModel)
+  {
+   eventMakerOnTheFly->SetV4RP(dV4RP);
+   eventMakerOnTheFly->SetV4SpreadRP(dV4SpreadRP);  
+  }
  } 
  
  // non-uniform acceptance:
@@ -661,6 +684,18 @@ void CheckUserSettings()
   break;
  }
  
+ if(bUseGlauberModel)
+ {
+  if(bPtDependentHarmonicV1||bEtaDependentHarmonicV1||
+     bPtDependentHarmonicV2||bEtaDependentHarmonicV2||
+     bPtDependentHarmonicV4||bEtaDependentHarmonicV4)
+  {   
+   cout<<"WARNING: When using Glauber model pt and/or eta dependence of flow harmonics is not supported !!!!"<<endl;
+   cout<<"         Set all booleans bPtDependentHarmonic* to kFALSE in the macro."<<endl;
+   exit(0); 
+  }
+ }
+
  if(bPtDependentHarmonicV4 && !bPtDependentHarmonicV2))
  {
   cout<<"WARNING: V4(pt,eta) is determined as pow(V2(pt,eta),2.) !!!!"<<endl;
@@ -688,7 +723,7 @@ void CheckUserSettings()
   cout<<"         You must also have bPtDependentHarmonicV2 = bPtDependentHarmonicV4 in the macro."<<endl;
   exit(0); 
  }
+  
 } // end of void CheckUserSettings()
 
 void LoadLibraries(const anaModes mode) {