]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/FLOW/AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.cxx
adding subevents onthefly
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowEventSimpleMakerOnTheFly.cxx
index 6892c98fe1faeab50b663c545f04b879e6cf2e35..fcb7575ff80e4ac3ac5e4672f9374ded8ab61b3f 100644 (file)
@@ -39,16 +39,44 @@ ClassImp(AliFlowEventSimpleMakerOnTheFly)
 
 
 AliFlowEventSimpleMakerOnTheFly::AliFlowEventSimpleMakerOnTheFly(UInt_t iseed):
+  fMultDistrOfRPsIsGauss(kFALSE),
   fMultiplicityOfRP(0),
   fMultiplicitySpreadOfRP(0.),
-  fPtFormula(NULL),
-  fPhiFormula(NULL),
+  fMinMultOfRP(0),
+  fMaxMultOfRP(0),  
+  fTemperatureOfRP(0.),  
+  fUseConstantHarmonics(kFALSE),
+  fV1RP(0.), 
+  fV1SpreadRP(0.), 
+  fV2DistrOfRPsIsGauss(kFALSE),
+  fV2RP(0.), 
+  fV2SpreadRP(0.), 
+  fMinV2RP(0.),
+  fMaxV2RP(0.),
+  fV4RP(0.), 
+  fV4SpreadRP(0.), 
+  fV2RPMax(0.), 
+  fPtCutOff(0.),
+  fPhiMin1(0.),              
+  fPhiMax1(0.),             
+  fProbability1(0.),       
+  fPhiMin2(0.),   
+  fPhiMax2(0.),            
+  fProbability2(0.),          
+  fPtSpectra(NULL),
+  fPhiDistribution(NULL),
   fMyTRandom3(NULL),
-  fCount(0)
- {
+  fCount(0),
+  fNoOfLoops(1),
+  fEtaMinA(-1.0),
+  fEtaMaxA(-0.01),
+  fEtaMinB(0.01),
+  fEtaMaxB(1.0)
+{
   // constructor
-   fMyTRandom3 = new TRandom3(iseed); 
- }
+  fMyTRandom3 = new TRandom3(iseed);   
+  gRandom->SetSeed(fMyTRandom3->Integer(65539));
+}
 
 
 //========================================================================
@@ -57,81 +85,239 @@ AliFlowEventSimpleMakerOnTheFly::AliFlowEventSimpleMakerOnTheFly(UInt_t iseed):
 AliFlowEventSimpleMakerOnTheFly::~AliFlowEventSimpleMakerOnTheFly()
 {
  // destructor
 if (fPtFormula) delete fPtFormula;
 if (fPhiFormula) delete fPhiFormula;
 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 of pions (Boltzman):   
+ 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(-pow(0.13957*0.13957+x*x,0.5)/[1])",dPtMin,dPtMax);  
+ fPtSpectra->SetParName(0,"Multiplicity of RPs");  
+ fPtSpectra->SetParName(1,"Temperature 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);
+   
+  // sampling the multiplicity:
+  Int_t iNewMultiplicityOfRP = fMultiplicityOfRP;
   
-  //reaction plane
-  Double_t fMCReactionPlaneAngle = fMyTRandom3->Uniform(0.,TMath::TwoPi());
-  
-  // pt:   
-  Double_t dPtMin = 0.; // to be improved 
-  Double_t dPtMax = 10.; // to be improved 
-  
-  fPtFormula = new TF1("PtFormula","[0]*x*TMath::Exp(-x*x)",dPtMin,dPtMax);
-  
-  fPtFormula->SetParName(0,"Multiplicity of RPs");
-  fPtFormula->SetParameter(0,fMultiplicityOfRP);
+  if(fMultDistrOfRPsIsGauss) {
+    if (fMultiplicitySpreadOfRP>0.0) iNewMultiplicityOfRP = (Int_t)fMyTRandom3->Gaus(fMultiplicityOfRP,fMultiplicitySpreadOfRP);
+    fPtSpectra->SetParameter(0,iNewMultiplicityOfRP);
+  } else {
+    if (fMinMultOfRP != fMaxMultOfRP) {
+      iNewMultiplicityOfRP = (Int_t)fMyTRandom3->Uniform(fMinMultOfRP,fMaxMultOfRP);
+      fPtSpectra->SetParameter(0,iNewMultiplicityOfRP);
+    }
+    else {
+      fPtSpectra->SetParameter(0,fMinMultOfRP);
+    }
+  }
+    
+  // 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);
 
-  // phi:
-  Double_t dPhiMin = 0.; // to be improved 
-  Double_t dPhiMax = TMath::TwoPi(); // to be improved 
-  
-  fPhiFormula = new TF1("phiDistribution","(1)*(1+2.*[0]*TMath::Cos(2*x))",dPhiMin,dPhiMax);
-  Double_t dV2 = 0.044; // to be improved
-  fPhiFormula->SetParName(0,"elliptic flow");
-  fPhiFormula->SetParameter(0,dV2);
+  // sampling the V1:
+  Double_t dNewV1RP=fV1RP;
+  if(fV1SpreadRP>0.0) {dNewV1RP = fMyTRandom3->Gaus(fV1RP,fV1SpreadRP);}
+  fPhiDistribution->SetParameter(0,dNewV1RP);
  
+  // sampling the V2:
+  if(fUseConstantHarmonics)
+  {
+   Double_t dNewV2RP = fV2RP;
+   if(fV2DistrOfRPsIsGauss)
+   {
+    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
+        {
+         dNewV2RP = fMinV2RP;
+         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 
   Double_t dEtaMax = 1.; // to be improved 
-
-
+  
   Int_t iGoodTracks = 0;
   Int_t iSelParticlesRP = 0;
   Int_t iSelParticlesPOI = 0;
-  for(Int_t i=0;i<fMultiplicityOfRP;i++) {
-    AliFlowTrackSimple* pTrack = new AliFlowTrackSimple();
-    pTrack->SetPt(fPtFormula->GetRandom());
-    pTrack->SetEta(fMyTRandom3->Uniform(dEtaMin,dEtaMax));
-    pTrack->SetPhi(fPhiFormula->GetRandom()+fMCReactionPlaneAngle);
-    pTrack->SetForRPSelection(kTRUE);
-    iSelParticlesRP++;
-    pTrack->SetForPOISelection(kTRUE);
-    iSelParticlesPOI++;
-
-    pEvent->TrackCollection()->Add(pTrack);
-    iGoodTracks++;
+  Double_t dTmpPt = 0.;
+  Double_t dTmpEta = 0.;
+  Double_t dTmpV2 = 0.;
+  Double_t dTmpPhi = 0.;
+  Bool_t bUniformAcceptance = kTRUE;
+  Double_t Pi = TMath::Pi();
+  if(!((fPhiMin1==0.) && (fPhiMax1==0.) && (fPhiMin2==0.) && (fPhiMax2==0.)))
+  {
+   bUniformAcceptance = kFALSE;
   }
+  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();
+      
+      // uniform acceptance:
+      if(bUniformAcceptance)
+      {
+       pTrack->SetPt(dTmpPt); 
+       pTrack->SetEta(dTmpEta); 
+       pTrack->SetPhi(dTmpPhi); 
+       pTrack->SetForRPSelection(kTRUE); 
+       iSelParticlesRP++; 
+       // assign particles to subevents
+       if (pTrack->Eta()>=fEtaMinA && pTrack->Eta()<=fEtaMaxA) {
+        pTrack->SetForSubevent(0);
+       }
+       if (pTrack->Eta()>=fEtaMinB && pTrack->Eta()<=fEtaMaxB) {
+        pTrack->SetForSubevent(1);
+       }
+       pTrack->SetForPOISelection(kTRUE); 
+       iSelParticlesPOI++; 
+       pEvent->TrackCollection()->Add(pTrack); 
+       iGoodTracks++; 
+      }
+      // non-uniform acceptance, 1st sector:
+      else if ((dTmpPhi > fPhiMin1*Pi/180) && (dTmpPhi < fPhiMax1*Pi/180))
+      {
+       if(fMyTRandom3->Uniform(0,1) > 1 - fProbability1)
+       {
+        pTrack->SetPt(dTmpPt);
+        pTrack->SetEta(dTmpEta);
+        pTrack->SetPhi(dTmpPhi);
+        pTrack->SetForRPSelection(kTRUE);
+        iSelParticlesRP++;
+       // assign particles to subevents
+       if (pTrack->Eta()>=fEtaMinA && pTrack->Eta()<=fEtaMaxA) {
+         pTrack->SetForSubevent(0);
+       }
+       if (pTrack->Eta()>=fEtaMinB && pTrack->Eta()<=fEtaMaxB) {
+         pTrack->SetForSubevent(1);
+       }
+        pTrack->SetForPOISelection(kTRUE);
+        iSelParticlesPOI++;
+        pEvent->TrackCollection()->Add(pTrack);
+        iGoodTracks++;
+       }
+      } 
+      // non-uniform acceptance, 2nd sector:
+      else if ((dTmpPhi > fPhiMin2*Pi/180) && (dTmpPhi < fPhiMax2*Pi/180))
+      {
+       if(fMyTRandom3->Uniform(0,1) > 1 - fProbability2)
+       {
+        pTrack->SetPt(dTmpPt);
+        pTrack->SetEta(dTmpEta);
+        pTrack->SetPhi(dTmpPhi);
+        pTrack->SetForRPSelection(kTRUE);
+        iSelParticlesRP++;
+       // assign particles to subevents
+       if (pTrack->Eta()>=fEtaMinA && pTrack->Eta()<=fEtaMaxA) {
+         pTrack->SetForSubevent(0);
+       }
+       if (pTrack->Eta()>=fEtaMinB && pTrack->Eta()<=fEtaMaxB) {
+         pTrack->SetForSubevent(1);
+       }
+       
+        pTrack->SetForPOISelection(kTRUE);
+        iSelParticlesPOI++;
+        pEvent->TrackCollection()->Add(pTrack);
+        iGoodTracks++;
+       }
+      } else 
+        {
+         pTrack->SetPt(dTmpPt);
+         pTrack->SetEta(dTmpEta);
+         pTrack->SetPhi(dTmpPhi);
+         pTrack->SetForRPSelection(kTRUE);
+         iSelParticlesRP++;
+        // assign particles to subevents
+        if (pTrack->Eta()>=fEtaMinA && pTrack->Eta()<=fEtaMaxA) {
+          pTrack->SetForSubevent(0);
+        }
+        if (pTrack->Eta()>=fEtaMinB && pTrack->Eta()<=fEtaMaxB) {
+          pTrack->SetForSubevent(1);
+        }
+         pTrack->SetForPOISelection(kTRUE);
+         iSelParticlesPOI++;
+         pEvent->TrackCollection()->Add(pTrack);
+         iGoodTracks++;
+        }
+    }
+  }
+  
+  // update the event quantities
   pEvent->SetEventNSelTracksRP(iSelParticlesRP);  
   pEvent->SetNumberOfTracks(iGoodTracks);//tracks used either for RP or for POI selection
-  pEvent->SetMCReactionPlaneAngle(fMCReactionPlaneAngle);
-
-  if (!fMCReactionPlaneAngle == 0) cout<<" MC Reaction Plane Angle = "<<  fMCReactionPlaneAngle << endl;
-  else cout<<" MC Reaction Plane Angle = unknown "<< endl;
+  pEvent->SetMCReactionPlaneAngle(dMCReactionPlaneAngle);
 
-  cout<<" iGoodTracks = "<< iGoodTracks << endl;
-  cout<<" # of RP selected tracks = "<<iSelParticlesRP<<endl;
-  cout<<" # of POI selected tracks = "<<iSelParticlesPOI<<endl;  
-  cout << "# " << ++fCount << " events processed" << endl;
+  if ( (++fCount % 100) == 0) {
+    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;
+  }
 
- return pEvent;  
 return pEvent;  
  
 } // end of CreateEventOnTheFly()