Merge branch 'feature-movesplit'

[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetModelBaseTask.cxx

diff --git a/PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx b/PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx
index 04a904804b3be0e7df9e9f7612eb67f4608f95e1..8e487b5b3ca7e2d1a34ce8a2ec9fa2079ebac237 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx
+++ b/PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx
@@ -11,6 +11,7 @@
 #include <TLorentzVector.h>
 #include <TRandom3.h>
 #include <TList.h>
+#include <TF2.h>
 
 #include "AliVEvent.h"
 #include "AliAODCaloCluster.h"
@@ -41,6 +42,7 @@ AliJetModelBaseTask::AliJetModelBaseTask() :
   fOutCellsName(),
   fMCParticlesName(),
   fOutMCParticlesName(),
+  fPythiaInfoName(""),
   fIsMC(kFALSE),
   fSuffix(),
   fEtaMin(-1),
@@ -55,8 +57,13 @@ AliJetModelBaseTask::AliJetModelBaseTask() :
   fNTracks(0),
   fMarkMC(99999),
   fPtSpectrum(0),
+  fPtPhiEvPlDistribution(0),
   fDensitySpectrum(0),
+  fDifferentialV2(0),
+  fAddV2(kFALSE),
+  fFlowFluctuations(kFALSE),
   fQAhistos(kFALSE),
+  fPsi(0),
   fIsInit(0),
   fGeom(0),
   fClusters(0),
@@ -72,7 +79,8 @@ AliJetModelBaseTask::AliJetModelBaseTask() :
   fOutMCParticlesMap(0),
   fMCLabelShift(0),
   fEsdMode(kFALSE),
-  fOutput(0)
+  fOutput(0),
+  fPythiaInfo(0x0)
 {
   // Default constructor.
 
@@ -93,6 +101,7 @@ AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) :
   fOutCellsName(""),
   fMCParticlesName(""),
   fOutMCParticlesName(""),
+  fPythiaInfoName(""),
   fIsMC(kFALSE),
   fSuffix("Processed"),
   fEtaMin(-1),
@@ -107,8 +116,13 @@ AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) :
   fNTracks(1),
   fMarkMC(99999),
   fPtSpectrum(0),
+  fPtPhiEvPlDistribution(0),
   fDensitySpectrum(0),
+  fDifferentialV2(0),
+  fAddV2(kFALSE),
+  fFlowFluctuations(kFALSE),
   fQAhistos(drawqa),
+  fPsi(0),
   fIsInit(0),
   fGeom(0),
   fClusters(0),
@@ -124,7 +138,8 @@ AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) :
   fOutMCParticlesMap(0),
   fMCLabelShift(0),
   fEsdMode(kFALSE),
-  fOutput(0)
+  fOutput(0),
+  fPythiaInfo(0x0)
 {
   // Standard constructor.
 
@@ -186,9 +201,9 @@ void AliJetModelBaseTask::UserExec(Option_t *)
   }
 
   if (fDensitySpectrum) {
-    fNTracks = fDensitySpectrum->GetRandom();
-    fNCells = fDensitySpectrum->GetRandom();
-    fNClusters = fDensitySpectrum->GetRandom();
+    fNTracks = TMath::Nint(fDensitySpectrum->GetRandom());
+    fNCells = TMath::Nint(fDensitySpectrum->GetRandom());
+    fNClusters = TMath::Nint(fDensitySpectrum->GetRandom());
   }
   
   // Clear map
@@ -205,8 +220,11 @@ void AliJetModelBaseTask::UserExec(Option_t *)
     }
   }
 
+  if (fPtPhiEvPlDistribution || fAddV2)
+    fPsi = gRandom->Rndm() * TMath::Pi();
+  
   Run();
-
+  
   if (fCaloCells && !fCopyArray) {
     delete fCaloCells;
     fCaloCells = tempCaloCells;
@@ -305,6 +323,10 @@ Bool_t AliJetModelBaseTask::ExecOnce()
     }
   }
 
+  if(fAddV2 && (!fDifferentialV2)) {
+    AliWarning(Form("%s: Cannot add v2 without diffential v2!", GetName()));
+  }
+
   if (!fCaloName.IsNull()) {
     fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
  
@@ -659,23 +681,30 @@ AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId, Int_t labe
 }
 
 //________________________________________________________________________
-AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi, Byte_t type, Double_t etaemc, Double_t phiemc, Double_t ptemc, Bool_t ise, Int_t label, Short_t charge)
+AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi, Byte_t type, Double_t etaemc, Double_t phiemc, Double_t ptemc, Bool_t ise, Int_t label, Short_t charge, Double_t mass)
 {
   // Add a track to the event.
-
-  const Int_t nTracks = fOutTracks->GetEntriesFast();
   
-  if (pt < 0) 
-    pt = GetRandomPt();
-  if (eta < -100) 
-    eta = GetRandomEta();
-  if (phi < 0) 
-    phi = GetRandomPhi();
+  if (pt < 0 && eta < -100 && phi < 0) {
+    GetRandomParticle(pt,eta,phi);
+  }
+  else {
+    if (pt < 0) 
+      pt = GetRandomPt();
+    if (eta < -100) 
+      eta = GetRandomEta();
+    if (phi < 0) 
+      phi = GetRandomPhi(pt);
+  }
 
   if (label >= 0)
     label += fMarkMC+fMCLabelShift;
   else if (label < 0)
     label -= fMarkMC+fMCLabelShift;
+  
+  if(fAddV2) AddV2(phi, pt);
+
+  const Int_t nTracks = fOutTracks->GetEntriesFast();
 
   AliPicoTrack *track = new ((*fOutTracks)[nTracks]) AliPicoTrack(pt, 
                                                                  eta, 
@@ -686,7 +715,8 @@ AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t
                                                                  etaemc, 
                                                                  phiemc,
                                                                  ptemc,
-                                                                 ise);
+                                                                 ise,
+                                                                 mass);
 
   return track;
 }
@@ -708,7 +738,25 @@ AliAODMCParticle* AliJetModelBaseTask::AddMCParticle(AliAODMCParticle *part, Int
   return aodpart;
 }
 
-//________________________________________________________________________
+//_____________________________________________________________________________
+void AliJetModelBaseTask::AddV2(Double_t &phi, Double_t &pt) const
+{
+    // similar to AliFlowTrackSimple::AddV2, except for the flow fluctuations
+    Double_t phi0(phi), v2(0.), f(0.), fp(0.), phiprev(0.);
+    if(fDifferentialV2) v2 = fDifferentialV2->Eval(pt);
+    if(TMath::AreEqualAbs(v2, 0, 1e-5)) return; 
+    // introduce flow fluctuations (gaussian)
+    if(fFlowFluctuations) v2 += TMath::Sqrt(2*(v2*.25)*(v2*.25))*TMath::ErfInverse(2*(gRandom->Uniform(0, fFlowFluctuations))-1);
+    for (Int_t i(0); i < 100; i++) {
+        phiprev=phi; //store last value for comparison
+        f =  phi-phi0+v2*TMath::Sin(2.*(phi-fPsi));
+        fp = 1.0+2.0*v2*TMath::Cos(2.*(phi-fPsi)); //first derivative
+        phi -= f/fp;
+        if (TMath::AreEqualAbs(phiprev, phi, 1e-10)) break; 
+    }
+}
+
+//_____________________________________________________________________________
 void AliJetModelBaseTask::CopyCells()
 {
   if (!fCaloCells)
@@ -893,7 +941,9 @@ Double_t AliJetModelBaseTask::GetRandomPhi(Bool_t emcal)
     if (phimin < EmcalMinPhi) phimin = EmcalMinPhi;
   }
 
-  return gRandom->Rndm() * (phimax - phimin) + phimin;
+  Double_t result = gRandom->Rndm() * (phimax - phimin) + phimin;
+
+  return result;
 }
 
 //________________________________________________________________________
@@ -907,6 +957,46 @@ Double_t AliJetModelBaseTask::GetRandomPt()
     return gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin;
 }
 
+//________________________________________________________________________
+void AliJetModelBaseTask::GetRandomParticle(Double_t &pt, Double_t &eta, Double_t &phi, Bool_t emcal)
+{
+  // Get a random particle.
+  
+  eta = GetRandomEta(emcal);
+
+  if (fPtPhiEvPlDistribution) {
+    Double_t phimax = fPhiMax;
+    Double_t phimin = fPhiMin;
+    
+    if (emcal) {
+      const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
+      const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
+      
+      if (phimax > EmcalMaxPhi) phimax = EmcalMaxPhi;
+      if (phimax < EmcalMinPhi) phimax = EmcalMinPhi;
+      if (phimin > EmcalMaxPhi) phimin = EmcalMaxPhi;
+      if (phimin < EmcalMinPhi) phimin = EmcalMinPhi;
+    }
+    
+    if (fPtPhiEvPlDistribution->GetXmin() > phimax || fPtPhiEvPlDistribution->GetXmax() < phimin) {
+      AliWarning(Form("The hisogram %s does not overlap with the EMCal acceptance limits. It will be ignored.",fPtPhiEvPlDistribution->GetName()));
+      pt = GetRandomPt();
+      phi = GetRandomPhi(emcal);
+    }
+    else {
+      do {
+       fPtPhiEvPlDistribution->GetRandom2(pt,phi);
+       phi += fPsi;
+       if (phi > TMath::Pi() * 2) phi -= TMath::Pi() * 2;
+      } while (phi > phimax || phi < phimin);
+    }
+  }
+  else {
+    pt = GetRandomPt();
+    phi = GetRandomPhi(emcal);
+  }
+}
+
 //________________________________________________________________________
 void AliJetModelBaseTask::Run() 
 {