]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetModelBaseTask.cxx
index fde3a44c6cf1123427a3b1d62b0a75e8bbe8904c..8e487b5b3ca7e2d1a34ce8a2ec9fa2079ebac237 100644 (file)
@@ -11,6 +11,7 @@
 #include <TLorentzVector.h>
 #include <TRandom3.h>
 #include <TList.h>
+#include <TF2.h>
 
 #include "AliVEvent.h"
 #include "AliAODCaloCluster.h"
@@ -25,6 +26,7 @@
 #include "AliPicoTrack.h"
 #include "AliEMCALGeometry.h"
 #include "AliLog.h"
+#include "AliNamedArrayI.h"
 
 ClassImp(AliJetModelBaseTask)
 
@@ -40,6 +42,8 @@ AliJetModelBaseTask::AliJetModelBaseTask() :
   fOutCellsName(),
   fMCParticlesName(),
   fOutMCParticlesName(),
+  fPythiaInfoName(""),
+  fIsMC(kFALSE),
   fSuffix(),
   fEtaMin(-1),
   fEtaMax(1),
@@ -53,7 +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),
@@ -69,9 +79,14 @@ AliJetModelBaseTask::AliJetModelBaseTask() :
   fOutMCParticlesMap(0),
   fMCLabelShift(0),
   fEsdMode(kFALSE),
-  fOutput(0)
+  fOutput(0),
+  fPythiaInfo(0x0)
 {
   // Default constructor.
+
+  fVertex[0] = 0;
+  fVertex[1] = 0;
+  fVertex[2] = 0;
 }
 
 //________________________________________________________________________
@@ -86,6 +101,8 @@ AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) :
   fOutCellsName(""),
   fMCParticlesName(""),
   fOutMCParticlesName(""),
+  fPythiaInfoName(""),
+  fIsMC(kFALSE),
   fSuffix("Processed"),
   fEtaMin(-1),
   fEtaMax(1),
@@ -99,7 +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),
@@ -115,13 +138,18 @@ AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) :
   fOutMCParticlesMap(0),
   fMCLabelShift(0),
   fEsdMode(kFALSE),
-  fOutput(0)
+  fOutput(0),
+  fPythiaInfo(0x0)
 {
   // Standard constructor.
 
   if (fQAhistos) {
     DefineOutput(1, TList::Class()); 
   }
+
+  fVertex[0] = 0;
+  fVertex[1] = 0;
+  fVertex[2] = 0;
 }
 
 //________________________________________________________________________
@@ -155,6 +183,14 @@ void AliJetModelBaseTask::UserExec(Option_t *)
   if (!fIsInit)
     return;
 
+  fVertex[0] = 0;
+  fVertex[1] = 0;
+  fVertex[2] = 0;
+
+  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
+  if (vert)
+    vert->GetXYZ(fVertex);
+
   if (fCopyArray) {
     if (fOutTracks)
       fOutTracks->Delete();
@@ -164,6 +200,16 @@ void AliJetModelBaseTask::UserExec(Option_t *)
       fOutMCParticles->Delete();
   }
 
+  if (fDensitySpectrum) {
+    fNTracks = TMath::Nint(fDensitySpectrum->GetRandom());
+    fNCells = TMath::Nint(fDensitySpectrum->GetRandom());
+    fNClusters = TMath::Nint(fDensitySpectrum->GetRandom());
+  }
+  
+  // Clear map
+  if (fOutMCParticlesMap)
+    fOutMCParticlesMap->Clear();
+
   AliVCaloCells *tempCaloCells = 0;
 
   if (fCaloCells) {
@@ -174,9 +220,12 @@ void AliJetModelBaseTask::UserExec(Option_t *)
     }
   }
 
+  if (fPtPhiEvPlDistribution || fAddV2)
+    fPsi = gRandom->Rndm() * TMath::Pi();
+  
   Run();
-
-  if (fCaloCells) {
+  
+  if (fCaloCells && !fCopyArray) {
     delete fCaloCells;
     fCaloCells = tempCaloCells;
   }
@@ -242,7 +291,7 @@ Bool_t AliJetModelBaseTask::ExecOnce()
     }
   }
 
-  if (fNTracks > 0 && !fTracksName.IsNull()) {
+  if (!fTracksName.IsNull()) {
     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
     if (!fTracks) {
       AliWarning(Form("%s: Couldn't retrieve tracks with name %s!", GetName(), fTracksName.Data()));
@@ -274,7 +323,11 @@ Bool_t AliJetModelBaseTask::ExecOnce()
     }
   }
 
-  if (fNClusters > 0 && !fCaloName.IsNull()) {
+  if(fAddV2 && (!fDifferentialV2)) {
+    AliWarning(Form("%s: Cannot add v2 without diffential v2!", GetName()));
+  }
+
+  if (!fCaloName.IsNull()) {
     fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
  
     if (!fClusters) {
@@ -290,8 +343,16 @@ Bool_t AliJetModelBaseTask::ExecOnce()
       fOutCaloName = fCaloName;
       if (fCopyArray) 
        fOutCaloName += fSuffix;
+      TString className;
+      if (fClusters)
+       className = fClusters->GetClass()->GetName();
+      else if (fEsdMode)
+       className = "AliESDCaloCluster";
+      else
+       className = "AliAODCaloCluster";
+       
       if (fCopyArray || !fClusters) {
-       fOutClusters = new TClonesArray(fClusters->GetClass()->GetName(), fClusters->GetSize());
+       fOutClusters = new TClonesArray(className.Data());
        fOutClusters->SetName(fOutCaloName);
        if (InputEvent()->FindListObject(fOutCaloName)) {
          AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCaloName.Data())); 
@@ -319,13 +380,13 @@ Bool_t AliJetModelBaseTask::ExecOnce()
        return kFALSE;
       }
       
-      fMCParticlesMap = dynamic_cast<TH1I*>(InputEvent()->FindListObject(fMCParticlesName + "_Map"));
+      fMCParticlesMap = dynamic_cast<AliNamedArrayI*>(InputEvent()->FindListObject(fMCParticlesName + "_Map"));
 
       if (!fMCParticlesMap) {
        AliWarning(Form("%s: Could not retrieve map for MC particles %s! Will assume MC labels consistent with indexes...", GetName(), fMCParticlesName.Data())); 
-       fMCParticlesMap = new TH1I(fMCParticlesName + "_Map", fMCParticlesName + "_Map",9999,0,1);
-       for (Int_t i = 0; i < 9999; i++) {
-         fMCParticlesMap->SetBinContent(i,i);
+       fMCParticlesMap = new AliNamedArrayI(fMCParticlesName + "_Map", 99999);
+       for (Int_t i = 0; i < 99999; i++) {
+         fMCParticlesMap->AddAt(i,i);
        }
       }
     }
@@ -346,7 +407,7 @@ Bool_t AliJetModelBaseTask::ExecOnce()
          InputEvent()->AddObject(fOutMCParticles);
        }
 
-       fOutMCParticlesMap = new TH1I(fOutMCParticlesName + "_Map", fOutMCParticlesName + "_Map",9999,0,1);
+       fOutMCParticlesMap = new AliNamedArrayI(fOutMCParticlesName + "_Map",99999);
        if (InputEvent()->FindListObject(fOutMCParticlesName + "_Map")) {
          AliFatal(Form("%s: Map %s_Map is already present in the event!", GetName(), fOutMCParticlesName.Data())); 
          return kFALSE;
@@ -362,37 +423,20 @@ Bool_t AliJetModelBaseTask::ExecOnce()
     }
   }
 
-  if (!fCaloName.IsNull() || !fCellsName.IsNull()) {
-    if (!fGeom) {
-      if (fGeomName.Length() > 0) {
-       fGeom = AliEMCALGeometry::GetInstance(fGeomName);
-       if (!fGeom) {
-         AliFatal(Form("Could not get geometry with name %s!", fGeomName.Data()));
-         return kFALSE;
-       }
-      } else {
-       fGeom = AliEMCALGeometry::GetInstance();
-       if (!fGeom) {
-         AliFatal("Could not get default geometry!");
-         return kFALSE;
-       }
+  if (!fGeom) {
+    if (fGeomName.Length() > 0) {
+      fGeom = AliEMCALGeometry::GetInstance(fGeomName);
+      if (!fGeom) {
+       AliFatal(Form("Could not get geometry with name %s!", fGeomName.Data()));
+       return kFALSE;
+      }
+    } else {
+      fGeom = AliEMCALGeometry::GetInstance();
+      if (!fGeom) {
+       AliFatal("Could not get default geometry!");
+       return kFALSE;
       }
     }
-  
-    const Double_t EmcalMinEta = fGeom->GetArm1EtaMin();
-    const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax();
-    const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad();
-    const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad();
-    
-    if (fEtaMax > EmcalMaxEta) fEtaMax = EmcalMaxEta;
-    if (fEtaMax < EmcalMinEta) fEtaMax = EmcalMinEta;
-    if (fEtaMin > EmcalMaxEta) fEtaMin = EmcalMaxEta;
-    if (fEtaMin < EmcalMinEta) fEtaMin = EmcalMinEta;
-    
-    if (fPhiMax > EmcalMaxPhi) fPhiMax = EmcalMaxPhi;
-    if (fPhiMax < EmcalMinPhi) fPhiMax = EmcalMinPhi;
-    if (fPhiMin > EmcalMaxPhi) fPhiMin = EmcalMaxPhi;
-    if (fPhiMin < EmcalMinPhi) fPhiMin = EmcalMinPhi;
   }
 
   return kTRUE;
@@ -414,28 +458,6 @@ Int_t AliJetModelBaseTask::SetNumberOfOutCells(Int_t n)
   return n;
 }
 
-//________________________________________________________________________
-void AliJetModelBaseTask::CopyCells()
-{
-  if (!fCaloCells)
-    return;
-
-  for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) {
-    Short_t mclabel = 0;
-    Double_t efrac = 0.;
-    Double_t time = -1;
-    Short_t cellNum = -1;
-    Double_t amp = -1;
-
-    fCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
-    fOutCaloCells->SetCell(i, cellNum, amp, time, mclabel, efrac);
-  }
-
-  fAddedCells = fCaloCells->GetNumberOfCells();
-
-  AliDebug(2, Form("%d cells from the PYTHIA event", fAddedCells));
-}
-
 //________________________________________________________________________
 Int_t AliJetModelBaseTask::AddCell(Double_t e, Double_t eta, Double_t phi)
 {
@@ -471,22 +493,107 @@ Int_t AliJetModelBaseTask::AddCell(Double_t e, Int_t absId, Double_t time, Int_t
 {
   // Add a cell to the event.
 
-  if (label == 0)
-    label = fMarkMC;
-  else
-    label += fMCLabelShift;
+  if (label < 0)
+    label = 0;
+
+  label += fMarkMC + fMCLabelShift;
+
+  Short_t pos = -1;
+  if (fCaloCells)  
+    pos = fCaloCells->GetCellPosition(absId);
+
+  Double_t efrac = 1;
+  Bool_t increaseOnSuccess = kFALSE;
+
+  if (pos < 0) {
+    increaseOnSuccess = kTRUE;
+    pos = fAddedCells;
+  }
+  else {
+    Short_t cellNumber = -1;
+    Double_t old_e = 0;
+    Double_t old_time = 0;
+    Int_t old_label = 0;
+    Double_t old_efrac = 0;
+    fOutCaloCells->GetCell(pos, cellNumber, old_e, old_time, old_label, old_efrac);
+    
+    efrac = e / (old_e + e);
+
+    if (old_label > 0 && e < old_efrac * old_e) {
+      label = old_label;
+      efrac = old_efrac;
+      time = old_time;
+    }
+    
+    e += old_e;
+  }
 
-  Bool_t r = fOutCaloCells->SetCell(fAddedCells, absId, e, time, label, 0);
+  Bool_t r = fOutCaloCells->SetCell(pos, absId, e, time, label, efrac);
 
   if (r) {
-    fAddedCells++;
+    if (increaseOnSuccess)
+      fAddedCells++;
     return fAddedCells;
   }
-  else {
+  else 
     return 0;
-  }
 }
 
+//________________________________________________________________________
+AliVCluster* AliJetModelBaseTask::AddCluster(AliVCluster *oc)
+{
+  // Add a cluster to the event.
+
+  const Int_t nClusters = fOutClusters->GetEntriesFast();
+  AliVCluster *dc = static_cast<AliVCluster*>(fOutClusters->New(nClusters));
+  dc->SetType(AliVCluster::kEMCALClusterv1);
+  dc->SetE(oc->E());
+  Float_t pos[3] = {0};
+  oc->GetPosition(pos);
+  dc->SetPosition(pos);
+  dc->SetNCells(oc->GetNCells());
+  dc->SetCellsAbsId(oc->GetCellsAbsId());
+  dc->SetCellsAmplitudeFraction(oc->GetCellsAmplitudeFraction());
+  dc->SetID(oc->GetID());
+  dc->SetDispersion(oc->GetDispersion());
+  dc->SetEmcCpvDistance(-1);
+  dc->SetChi2(-1);
+  dc->SetTOF(oc->GetTOF());     //time-of-flight
+  dc->SetNExMax(oc->GetNExMax()); //number of local maxima
+  dc->SetM02(oc->GetM02());
+  dc->SetM20(oc->GetM20());
+  dc->SetDistanceToBadChannel(oc->GetDistanceToBadChannel()); 
+
+  //MC
+  UInt_t nlabels = oc->GetNLabels();
+  Int_t *labels = oc->GetLabels();
+  TArrayI parents;
+  if (nlabels != 0 && labels) 
+    parents.Set(nlabels, labels);
+  else {
+    nlabels = 1;
+    parents.Set(1);
+    parents[0] = 0;
+  }
+
+  if (fMarkMC+fMCLabelShift != 0) {
+    for (UInt_t i = 0; i < nlabels; i++) {
+      parents[i] += fMarkMC+fMCLabelShift;
+    }
+  }
+
+  AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(dc);
+  if (esdClus) {
+    esdClus->AddLabels(parents); 
+  }
+  else {
+    AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(dc);
+    if (aodClus) 
+      aodClus->SetLabel(parents.GetArray(), nlabels); 
+  }
+
+  return dc;
+}
 
 //________________________________________________________________________
 AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Double_t eta, Double_t phi, Int_t label)
@@ -555,10 +662,10 @@ AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId, Int_t labe
   cluster->SetEmcCpvDistance(-1);
 
   //MC label
-  if (label == 0)
-    label = fMarkMC;
-  else
-    label += fMCLabelShift;
+  if (label < 0)
+    label = 0;
+  label += fMarkMC+fMCLabelShift;
 
   if (fEsdMode) {
     AliESDCaloCluster *esdClus = static_cast<AliESDCaloCluster*>(cluster);
@@ -574,33 +681,42 @@ 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, Bool_t ise, Int_t label)
+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);
 
-  if (label == 0)
-    label = fMarkMC;
-  else
-    label += fMCLabelShift;
+  const Int_t nTracks = fOutTracks->GetEntriesFast();
 
   AliPicoTrack *track = new ((*fOutTracks)[nTracks]) AliPicoTrack(pt, 
                                                                  eta, 
                                                                  phi, 
-                                                                 1,
+                                                                 charge,
                                                                  label,
                                                                  type, 
                                                                  etaemc, 
-                                                                 phiemc, 
-                                                                 ise);
+                                                                 phiemc,
+                                                                 ptemc,
+                                                                 ise,
+                                                                 mass);
 
   return track;
 }
@@ -612,11 +728,61 @@ AliAODMCParticle* AliJetModelBaseTask::AddMCParticle(AliAODMCParticle *part, Int
 
   AliAODMCParticle *aodpart = new ((*fOutMCParticles)[nPart]) AliAODMCParticle(*part);
 
-  fOutMCParticlesMap->SetBinContent(origIndex + fMCLabelShift, nPart);
+  if (origIndex + fMCLabelShift >= fOutMCParticlesMap->GetSize())
+    fOutMCParticlesMap->Set((origIndex + fMCLabelShift)*2);
+
+  fOutMCParticlesMap->AddAt(nPart, origIndex + fMCLabelShift);
+  AliDebug(2, Form("Setting bin %d to %d (fMCLabelShift=%d, origIndex=%d)", 
+                  origIndex + fMCLabelShift, fOutMCParticlesMap->At(origIndex + fMCLabelShift), fMCLabelShift, origIndex));
 
   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)
+    return;
+
+  fAddedCells = 0;
+  fCaloCells->Sort();
+  for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) {
+    Int_t mclabel = 0;
+    Double_t efrac = 0.;
+    Double_t time = -1;
+    Short_t cellNum = -1;
+    Double_t amp = -1;
+
+    fCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
+
+    if (!fIsMC) 
+      mclabel = 0;
+
+    fOutCaloCells->SetCell(i, cellNum, amp, time, mclabel, efrac);
+    fAddedCells++;
+  }
+
+  AliDebug(2, Form("%d cells from the current event", fAddedCells));
+}
+
 //________________________________________________________________________
 void AliJetModelBaseTask::CopyClusters()
 {
@@ -632,7 +798,12 @@ void AliJetModelBaseTask::CopyClusters()
       AliESDCaloCluster *esdcluster = static_cast<AliESDCaloCluster*>(fClusters->At(i));
       if (!esdcluster || !esdcluster->IsEMCAL())
        continue;
-      new ((*fOutClusters)[nCopiedClusters]) AliESDCaloCluster(*esdcluster);
+      AliESDCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliESDCaloCluster(*esdcluster);
+      if (!fIsMC) {
+       TArrayI *labels = clus->GetLabelsArray();
+       if (labels)
+         labels->Reset();
+      }
       nCopiedClusters++;
     }
   }
@@ -641,7 +812,9 @@ void AliJetModelBaseTask::CopyClusters()
       AliAODCaloCluster *aodcluster = static_cast<AliAODCaloCluster*>(fClusters->At(i));
       if (!aodcluster || !aodcluster->IsEMCAL())
        continue;
-      new ((*fOutClusters)[nCopiedClusters]) AliAODCaloCluster(*aodcluster);
+      AliAODCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliAODCaloCluster(*aodcluster);
+      if (!fIsMC) 
+       clus->SetLabel(0,0);
       nCopiedClusters++;
     }
   }
@@ -658,10 +831,12 @@ void AliJetModelBaseTask::CopyTracks()
   const Int_t nTracks = fTracks->GetEntriesFast();
   Int_t nCopiedTracks = 0;
   for (Int_t i = 0; i < nTracks; ++i) {
-    AliPicoTrack *track = static_cast<AliPicoTrack*>(fTracks->At(i));
-    if (!track)
+    AliPicoTrack *picotrack = static_cast<AliPicoTrack*>(fTracks->At(i));
+    if (!picotrack)
       continue;
-    new ((*fOutTracks)[nCopiedTracks]) AliPicoTrack(*track);
+    AliPicoTrack *track = new ((*fOutTracks)[nCopiedTracks]) AliPicoTrack(*picotrack);
+    if (!fIsMC && track->GetLabel() != 0) 
+      track->SetLabel(0);
     nCopiedTracks++;
   }
 }
@@ -685,18 +860,19 @@ void AliJetModelBaseTask::CopyMCParticles()
     nCopiedPart++;
   }
 
-  if (!fMCParticlesMap)
+  if (!fMCParticlesMap || !fOutMCParticlesMap)
     return;
 
-  Int_t shift = 0;
+  if (fOutMCParticlesMap->GetSize() < fMCParticlesMap->GetSize())
+    fOutMCParticlesMap->Set(fMCParticlesMap->GetSize() * 2);
 
-  for (Int_t i = 0; i < fMCParticlesMap->GetNbinsX()+2; i++) {
-    fOutMCParticlesMap->SetBinContent(i, fMCParticlesMap->GetBinContent(i));
-    if (fMCParticlesMap->GetBinContent(i) != 0)
-      shift = i;
+  for (Int_t i = 0; i < fMCParticlesMap->GetSize(); i++) {
+    fOutMCParticlesMap->AddAt(fMCParticlesMap->At(i), i);
+    if (fMCParticlesMap->At(i) >= 0)
+      fMCLabelShift = i;
   }
 
-  fMCLabelShift = shift;
+  AliDebug(2,Form("MC particles copied. fMCLabelShift=%d",fMCLabelShift));
 }
 
 //________________________________________________________________________
@@ -709,9 +885,9 @@ void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &abs
   Double_t rndPhi = phi;
   do {
     if (eta < -100)
-      rndEta = GetRandomEta();
+      rndEta = GetRandomEta(kTRUE);
     if (phi < 0)
-      rndPhi = GetRandomPhi();
+      rndPhi = GetRandomPhi(kTRUE);
     fGeom->GetAbsCellIdFromEtaPhi(rndEta, rndPhi, absId);  
     repeats++;
   } while (absId == -1 && repeats < 100);
@@ -727,19 +903,47 @@ void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &abs
 }
 
 //________________________________________________________________________
-Double_t AliJetModelBaseTask::GetRandomEta()
+Double_t AliJetModelBaseTask::GetRandomEta(Bool_t emcal)
 {
   // Get random eta.
 
-  return gRandom->Rndm() * (fEtaMax - fEtaMin) + fEtaMin;
+  Double_t etamax = fEtaMax;
+  Double_t etamin = fEtaMin;
+
+  if (emcal) {
+    const Double_t EmcalMinEta = fGeom->GetArm1EtaMin();
+    const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax();
+    
+    if (etamax > EmcalMaxEta) etamax = EmcalMaxEta;
+    if (etamax < EmcalMinEta) etamax = EmcalMinEta;
+    if (etamin > EmcalMaxEta) etamin = EmcalMaxEta;
+    if (etamin < EmcalMinEta) etamin = EmcalMinEta;
+  }
+
+  return gRandom->Rndm() * (etamax - etamin) + etamin;
 }
 
 //________________________________________________________________________
-Double_t AliJetModelBaseTask::GetRandomPhi()
+Double_t AliJetModelBaseTask::GetRandomPhi(Bool_t emcal)
 {
   // Get random phi.
+  
+  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;
+  }
 
-  return gRandom->Rndm() * (fPhiMax - fPhiMin) + fPhiMin;
+  Double_t result = gRandom->Rndm() * (phimax - phimin) + phimin;
+
+  return result;
 }
 
 //________________________________________________________________________
@@ -753,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() 
 {