up from salvatore
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Feb 2013 16:58:49 +0000 (16:58 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Feb 2013 16:58:49 +0000 (16:58 +0000)
PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx
PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.h
PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx
PWGJE/EMCALJetTasks/AliJetModelBaseTask.h
PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.cxx
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.h
PWGJE/EMCALJetTasks/macros/AddTaskJetEmbeddingFromAOD.C

index 23c424d..c36ea78 100644 (file)
@@ -23,6 +23,7 @@
 #include "AliVTrack.h"
 #include "AliVCluster.h"
 #include "AliVCaloCells.h"
+#include "AliAODMCParticle.h"
 #include "AliCentrality.h"
 #include "AliVHeader.h"
 #include "AliVVertex.h"
@@ -42,6 +43,7 @@ AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() :
   fAODTrackName(),
   fAODClusName(),
   fAODCellsName(),
+  fAODMCParticlesName(),
   fMinCentrality(0),
   fMaxCentrality(10),
   fTriggerMask(AliVEvent::kAny),
@@ -58,6 +60,7 @@ AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() :
   fAODTracks(0),
   fAODClusters(0),
   fAODCaloCells(0),
+  fAODMCParticles(0),
   fHistFileIDs(0)
 {
   // Default constructor.
@@ -75,8 +78,9 @@ AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t
   fAODHeaderName("header"),
   fAODVertexName("vertices"),
   fAODTrackName("tracks"),
-  fAODClusName("caloClusters"),
+  fAODClusName(""),
   fAODCellsName("emcalCells"),
+  fAODMCParticlesName(AliAODMCParticle::StdBranchName()),
   fMinCentrality(0),
   fMaxCentrality(10),
   fTriggerMask(AliVEvent::kAny),
@@ -92,7 +96,8 @@ AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t
   fAODVertex(0),
   fAODTracks(0),
   fAODClusters(0),
-  fAODCaloCells(0),
+  fAODCaloCells(0),  
+  fAODMCParticles(0),
   fHistFileIDs(0)
 {
   // Standard constructor.
@@ -232,6 +237,9 @@ Bool_t AliJetEmbeddingFromAODTask::GetNextEntry()
       
       if (!fAODCellsName.IsNull()) 
        tree->SetBranchAddress(fAODCellsName, &fAODCaloCells);
+
+      if (!fAODMCParticlesName.IsNull()) 
+       tree->SetBranchAddress(fAODMCParticlesName, &fAODMCParticles);
       
       entry = 0;
     }
@@ -289,9 +297,30 @@ void AliJetEmbeddingFromAODTask::Run()
     return;
   }
 
-  if (fTracks && fOutTracks) {
+  if (fOutMCParticles) {
+
+    if (fCopyArray && fMCParticles)
+      CopyMCParticles();
+
+    if (fAODMCParticles) {
+      for (Int_t i = 0; i < fAODMCParticles->GetEntriesFast(); i++) {
+       AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fAODMCParticles->At(i));
+       if (!part) {
+         AliError(Form("Could not find MC particle %d in branch %s of tree %s!", i, fAODMCParticlesName.Data(), fAODTreeName.Data()));
+         continue;
+       }
+       
+       if (!part->IsPhysicalPrimary()) 
+         continue;
+
+       AddMCParticle(part, i);
+      }
+    }
+  }
+
+  if (fOutTracks) {
 
-    if (fCopyArray)
+    if (fCopyArray && fTracks)
       CopyTracks();
 
     if (fAODTracks) {
@@ -342,14 +371,14 @@ void AliJetEmbeddingFromAODTask::Run()
        }
        
        AliDebug(3, Form("Embedding track with pT = %f, eta = %f, phi = %f", track->Pt(), track->Eta(), track->Phi()));
-       AddTrack(track->Pt(), track->Eta(), track->Phi(), type, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), isEmc);
+       AddTrack(track->Pt(), track->Eta(), track->Phi(), type, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), isEmc, track->GetLabel());
       }
     }
   }
 
-  if (fClusters && fOutClusters) {
+  if (fOutClusters) {
 
-    if (fCopyArray)
+    if (fCopyArray && fClusters)
       CopyClusters();
 
     if (fAODClusters) {
@@ -362,20 +391,24 @@ void AliJetEmbeddingFromAODTask::Run()
        TLorentzVector vect;
        Double_t vert[3] = {0,0,0};
        clus->GetMomentum(vect,vert);
-       AddCluster(clus->E(), vect.Eta(), vect.Phi());
+       AddCluster(clus->E(), vect.Eta(), vect.Phi(), clus->GetLabel());
       }
     }
   }
 
-  if (fCaloCells && fOutCaloCells) {
+  if (fOutCaloCells) {
+
+    Int_t totalCells = 0;
 
-    Int_t totalCells = fCaloCells->GetNumberOfCells();
+    if (fCaloCells)
+      totalCells += fCaloCells->GetNumberOfCells();
     if (fAODCaloCells)
       totalCells += fAODCaloCells->GetNumberOfCells();
 
     SetNumberOfOutCells(totalCells);
 
-    CopyCells();
+    if (fCaloCells)
+      CopyCells();
 
     if (fAODCaloCells) {
       for (Short_t i = 0; i < fAODCaloCells->GetNumberOfCells(); i++) {
@@ -387,7 +420,7 @@ void AliJetEmbeddingFromAODTask::Run()
        
        fAODCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac);
        AliDebug(3,Form("Adding cell with amplitude %f, absolute ID %d, time %f", amp, cellNum, time));
-       AddCell(amp, cellNum, time);
+       AddCell(amp, cellNum, time, mclabel);
       }
     }
 
index 3555413..6b10895 100644 (file)
@@ -28,6 +28,7 @@ class AliJetEmbeddingFromAODTask : public AliJetModelBaseTask {
   void           SetAODTracksName(const char *n)                   { fAODTrackName       = n     ; }
   void           SetAODClusName(const char *n)                     { fAODClusName        = n     ; }
   void           SetAODCellsName(const char *n)                    { fAODCellsName       = n     ; }
+  void           SetAODMCParticlesName(const char *n)              { fAODMCParticlesName = n     ; }
   void           SetCentralityRange(Double_t min, Double_t max)    { fMinCentrality      = min   ; fMaxCentrality    = max ; }
   void           SetTriggerMask(UInt_t mask)                       { fTriggerMask        = mask  ; }
   void           SetAODfilterBits(Int_t b0 = 0, Int_t b1 = 0)      { fAODfilterBits[0]   = b0    ; fAODfilterBits[1] = b1  ; }
@@ -41,36 +42,38 @@ class AliJetEmbeddingFromAODTask : public AliJetModelBaseTask {
   Bool_t         GetNextEntry()        ;// get next entry in current tree
   Bool_t         IsAODEventSelected()  ;// AOD event trigger/centrality selection
 
-  TObjArray     *fFileList         ;//  List of AOD files
-  TString        fAODTreeName      ;//  Name of the tree in the AOD file
-  TString        fAODHeaderName    ;//  Name of the header in the AOD tree
-  TString        fAODVertexName    ;//  Name of the vertex branch in the AOD tree
-  TString        fAODTrackName     ;//  Name of the track collection branch in the AOD tree
-  TString        fAODClusName      ;//  Name of the cluster collection branch in the AOD tree
-  TString        fAODCellsName     ;//  Name of the cell collection branch in the AOD tree
-  Double_t       fMinCentrality    ;//  Minimum centrality
-  Double_t       fMaxCentrality    ;//  Maximum centrality
-  UInt_t         fTriggerMask      ;//  Trigger selection mask
-  Double_t       fZVertexCut       ;//  Z vertex cut
-  Int_t          fAODfilterBits[2] ;//  AOD track filter bit map
-  Bool_t         fIncludeNoITS     ;//  True = includes tracks with failed ITS refit
-  Int_t          fTotalFiles       ;//  Total number of files per pt hard bin
-  Bool_t         fEsdTreeMode      ;//! True = embed from ESD (must be a skimmed ESD!)
-  Int_t          fCurrentFileID    ;//! Current file being processed (trough the event handler)
-  Int_t          fCurrentAODFileID ;//! Current file ID
-  TFile         *fCurrentAODFile   ;//! Current open file
-  Int_t          fPicoTrackVersion ;//! Version of the PicoTrack class (if any) in fCurrentAODFile
-  AliVHeader    *fAODHeader        ;//! AOD header
-  TClonesArray  *fAODVertex        ;//! AOD vertex
-  TClonesArray  *fAODTracks        ;//! AOD track collection
-  TClonesArray  *fAODClusters      ;//! AOD cluster collection
-  AliVCaloCells *fAODCaloCells     ;//! AOD cell collection
-  TH2           *fHistFileIDs      ;//! Current file ID vs. AOD file ID (to be embedded)
+  TObjArray     *fFileList            ;//  List of AOD files
+  TString        fAODTreeName         ;//  Name of the tree in the AOD file
+  TString        fAODHeaderName       ;//  Name of the header in the AOD tree
+  TString        fAODVertexName       ;//  Name of the vertex branch in the AOD tree
+  TString        fAODTrackName        ;//  Name of the track collection branch in the AOD tree
+  TString        fAODClusName         ;//  Name of the cluster collection branch in the AOD tree
+  TString        fAODCellsName        ;//  Name of the cell collection branch in the AOD tree
+  TString        fAODMCParticlesName  ;//  Name of the cell collection branch in the AOD tree
+  Double_t       fMinCentrality       ;//  Minimum centrality
+  Double_t       fMaxCentrality       ;//  Maximum centrality
+  UInt_t         fTriggerMask         ;//  Trigger selection mask
+  Double_t       fZVertexCut          ;//  Z vertex cut
+  Int_t          fAODfilterBits[2]    ;//  AOD track filter bit map
+  Bool_t         fIncludeNoITS        ;//  True = includes tracks with failed ITS refit
+  Int_t          fTotalFiles          ;//  Total number of files per pt hard bin
+  Bool_t         fEsdTreeMode         ;//! True = embed from ESD (must be a skimmed ESD!)
+  Int_t          fCurrentFileID       ;//! Current file being processed (trough the event handler)
+  Int_t          fCurrentAODFileID    ;//! Current file ID
+  TFile         *fCurrentAODFile      ;//! Current open file
+  Int_t          fPicoTrackVersion    ;//! Version of the PicoTrack class (if any) in fCurrentAODFile
+  AliVHeader    *fAODHeader           ;//! AOD header
+  TClonesArray  *fAODVertex           ;//! AOD vertex
+  TClonesArray  *fAODTracks           ;//! AOD track collection
+  TClonesArray  *fAODClusters         ;//! AOD cluster collection
+  AliVCaloCells *fAODCaloCells        ;//! AOD cell collection
+  TClonesArray  *fAODMCParticles      ;//! AOD MC particles collection
+  TH2           *fHistFileIDs         ;//! Current file ID vs. AOD file ID (to be embedded)
 
  private:
   AliJetEmbeddingFromAODTask(const AliJetEmbeddingFromAODTask&);            // not implemented
   AliJetEmbeddingFromAODTask &operator=(const AliJetEmbeddingFromAODTask&); // not implemented
 
-  ClassDef(AliJetEmbeddingFromAODTask, 2) // Jet embedding from AOD task
+  ClassDef(AliJetEmbeddingFromAODTask, 3) // Jet embedding from AOD task
 };
 #endif
index bc201b0..fde3a44 100644 (file)
@@ -7,7 +7,7 @@
 #include "AliJetModelBaseTask.h"
 
 #include <TClonesArray.h>
-#include <TH1.h>
+#include <TH1I.h>
 #include <TLorentzVector.h>
 #include <TRandom3.h>
 #include <TList.h>
@@ -20,6 +20,7 @@
 #include "AliEMCALRecPoint.h"
 #include "AliESDCaloCells.h"
 #include "AliAODCaloCells.h"
+#include "AliAODMCParticle.h"
 #include "AliVCaloCells.h"
 #include "AliPicoTrack.h"
 #include "AliEMCALGeometry.h"
@@ -37,6 +38,8 @@ AliJetModelBaseTask::AliJetModelBaseTask() :
   fOutCaloName(),
   fCellsName(),
   fOutCellsName(),
+  fMCParticlesName(),
+  fOutMCParticlesName(),
   fSuffix(),
   fEtaMin(-1),
   fEtaMax(1),
@@ -60,6 +63,11 @@ AliJetModelBaseTask::AliJetModelBaseTask() :
   fCaloCells(0),
   fOutCaloCells(0),
   fAddedCells(0),
+  fMCParticles(0),
+  fMCParticlesMap(0),
+  fOutMCParticles(0),
+  fOutMCParticlesMap(0),
+  fMCLabelShift(0),
   fEsdMode(kFALSE),
   fOutput(0)
 {
@@ -76,6 +84,8 @@ AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) :
   fOutCaloName("CaloClustersCorrEmbedded"),
   fCellsName(""),
   fOutCellsName(""),
+  fMCParticlesName(""),
+  fOutMCParticlesName(""),
   fSuffix("Processed"),
   fEtaMin(-1),
   fEtaMax(1),
@@ -99,6 +109,11 @@ AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) :
   fCaloCells(0),
   fOutCaloCells(0),
   fAddedCells(0),
+  fMCParticles(0),
+  fMCParticlesMap(0),
+  fOutMCParticles(0),
+  fOutMCParticlesMap(0),
+  fMCLabelShift(0),
   fEsdMode(kFALSE),
   fOutput(0)
 {
@@ -145,6 +160,8 @@ void AliJetModelBaseTask::UserExec(Option_t *)
       fOutTracks->Delete();
     if (fOutClusters)
       fOutClusters->Delete();
+    if (fOutMCParticles)
+      fOutMCParticles->Delete();
   }
 
   AliVCaloCells *tempCaloCells = 0;
@@ -190,14 +207,12 @@ Bool_t AliJetModelBaseTask::ExecOnce()
     fPhiMax = fPhiMin;
   }
 
-  if (!fCaloCells && !fCellsName.IsNull()) {
+  if (!fCellsName.IsNull()) {
     fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCellsName));
     if (!fCaloCells) {
-      AliError(Form("%s: Couldn't retrieve calo cells with name %s!", GetName(), fCellsName.Data()));
-      return kFALSE;
+      AliWarning(Form("%s: Couldn't retrieve calo cells with name %s!", GetName(), fCellsName.Data()));
     }
-    
-    if (!fCaloCells->InheritsFrom("AliVCaloCells")) {
+    else if (!fCaloCells->InheritsFrom("AliVCaloCells")) {
       AliError(Form("%s: Collection %s does not contain a AliVCaloCells object!", GetName(), fCellsName.Data())); 
       fCaloCells = 0;
       return kFALSE;
@@ -205,9 +220,9 @@ Bool_t AliJetModelBaseTask::ExecOnce()
 
     if (!fOutCaloCells) {
       fOutCellsName = fCellsName;
-      if (fCopyArray) {
+      if (fCopyArray) 
        fOutCellsName += fSuffix;
-
+      if (fCopyArray || !fCaloCells) {
        if (fEsdMode) 
          fOutCaloCells = new AliESDCaloCells(fOutCellsName,fOutCellsName);
        else
@@ -227,14 +242,12 @@ Bool_t AliJetModelBaseTask::ExecOnce()
     }
   }
 
-  if (fNTracks > 0 && !fTracks && !fTracksName.IsNull()) {
+  if (fNTracks > 0 && !fTracksName.IsNull()) {
     fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
     if (!fTracks) {
-      AliError(Form("%s: Couldn't retrieve tracks with name %s!", GetName(), fTracksName.Data()));
-      return kFALSE;
+      AliWarning(Form("%s: Couldn't retrieve tracks with name %s!", GetName(), fTracksName.Data()));
     }
-    
-    if (!fTracks->GetClass()->GetBaseClass("AliPicoTrack")) {
+    else if (!fTracks->GetClass()->GetBaseClass("AliPicoTrack")) {
       AliError(Form("%s: Collection %s does not contain AliPicoTrack objects!", GetName(), fTracksName.Data())); 
       fTracks = 0;
       return kFALSE;
@@ -242,9 +255,10 @@ Bool_t AliJetModelBaseTask::ExecOnce()
 
     if (!fOutTracks) {
       fOutTracksName = fTracksName;
-      if (fCopyArray) {
+      if (fCopyArray)
        fOutTracksName += fSuffix;
-       fOutTracks = new TClonesArray("AliPicoTrack", fTracks->GetSize());
+      if (fCopyArray || !fTracks) {
+       fOutTracks = new TClonesArray("AliPicoTrack");
        fOutTracks->SetName(fOutTracksName);
        if (InputEvent()->FindListObject(fOutTracksName)) {
          AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutTracksName.Data())); 
@@ -260,15 +274,13 @@ Bool_t AliJetModelBaseTask::ExecOnce()
     }
   }
 
-  if (fNClusters > 0 && !fClusters && !fCaloName.IsNull()) {
+  if (fNClusters > 0 && !fCaloName.IsNull()) {
     fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
  
     if (!fClusters) {
-      AliError(Form("%s: Couldn't retrieve clusters with name %s!", GetName(), fCaloName.Data()));
-      return kFALSE;
+      AliWarning(Form("%s: Couldn't retrieve clusters with name %s!", GetName(), fCaloName.Data()));
     }
-
-    if (!fClusters->GetClass()->GetBaseClass("AliVCluster")) {
+    else if (!fClusters->GetClass()->GetBaseClass("AliVCluster")) {
       AliError(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloName.Data())); 
       fClusters = 0;
       return kFALSE;
@@ -276,8 +288,9 @@ Bool_t AliJetModelBaseTask::ExecOnce()
 
     if (!fOutClusters) {
       fOutCaloName = fCaloName;
-      if (fCopyArray) {
+      if (fCopyArray) 
        fOutCaloName += fSuffix;
+      if (fCopyArray || !fClusters) {
        fOutClusters = new TClonesArray(fClusters->GetClass()->GetName(), fClusters->GetSize());
        fOutClusters->SetName(fOutCaloName);
        if (InputEvent()->FindListObject(fOutCaloName)) {
@@ -294,6 +307,61 @@ Bool_t AliJetModelBaseTask::ExecOnce()
     }
   }
 
+  if (!fMCParticlesName.IsNull()) {
+    fMCParticles = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCParticlesName));
+    if (!fMCParticles) {
+      AliWarning(Form("%s: Couldn't retrieve MC particles with name %s!", GetName(), fMCParticlesName.Data()));
+    }
+    else {
+      if (!fMCParticles->GetClass()->GetBaseClass("AliAODMCParticle")) {
+       AliError(Form("%s: Collection %s does not contain AliAODMCParticle objects!", GetName(), fMCParticlesName.Data())); 
+       fMCParticles = 0;
+       return kFALSE;
+      }
+      
+      fMCParticlesMap = dynamic_cast<TH1I*>(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);
+       }
+      }
+    }
+
+    if (!fOutMCParticles) {
+      fOutMCParticlesName = fMCParticlesName;
+      if (fCopyArray)
+       fOutMCParticlesName += fSuffix;
+      if (fCopyArray || !fMCParticles) {
+
+       fOutMCParticles = new TClonesArray("AliAODMCParticle");
+       fOutMCParticles->SetName(fOutMCParticlesName);
+       if (InputEvent()->FindListObject(fOutMCParticlesName)) {
+         AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutMCParticlesName.Data())); 
+         return kFALSE;
+       }
+       else {
+         InputEvent()->AddObject(fOutMCParticles);
+       }
+
+       fOutMCParticlesMap = new TH1I(fOutMCParticlesName + "_Map", fOutMCParticlesName + "_Map",9999,0,1);
+       if (InputEvent()->FindListObject(fOutMCParticlesName + "_Map")) {
+         AliFatal(Form("%s: Map %s_Map is already present in the event!", GetName(), fOutMCParticlesName.Data())); 
+         return kFALSE;
+       }
+       else {
+         InputEvent()->AddObject(fOutMCParticlesMap);
+       }
+      }
+      else {
+       fOutMCParticles = fMCParticles;
+       fOutMCParticlesMap = fMCParticlesMap;
+      }
+    }
+  }
+
   if (!fCaloName.IsNull() || !fCellsName.IsNull()) {
     if (!fGeom) {
       if (fGeomName.Length() > 0) {
@@ -349,6 +417,9 @@ Int_t AliJetModelBaseTask::SetNumberOfOutCells(Int_t n)
 //________________________________________________________________________
 void AliJetModelBaseTask::CopyCells()
 {
+  if (!fCaloCells)
+    return;
+
   for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) {
     Short_t mclabel = 0;
     Double_t efrac = 0.;
@@ -396,11 +467,16 @@ Int_t AliJetModelBaseTask::AddCell(Double_t e, Double_t eta, Double_t phi)
 }
 
 //________________________________________________________________________
-Int_t AliJetModelBaseTask::AddCell(Double_t e, Int_t absId, Double_t time)
+Int_t AliJetModelBaseTask::AddCell(Double_t e, Int_t absId, Double_t time, Int_t label)
 {
   // Add a cell to the event.
 
-  Bool_t r = fOutCaloCells->SetCell(fAddedCells, absId, e, time, fMarkMC, 0);
+  if (label == 0)
+    label = fMarkMC;
+  else
+    label += fMCLabelShift;
+
+  Bool_t r = fOutCaloCells->SetCell(fAddedCells, absId, e, time, label, 0);
 
   if (r) {
     fAddedCells++;
@@ -413,7 +489,7 @@ Int_t AliJetModelBaseTask::AddCell(Double_t e, Int_t absId, Double_t time)
 
 
 //________________________________________________________________________
-AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Double_t eta, Double_t phi)
+AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Double_t eta, Double_t phi, Int_t label)
 {
   // Add a cluster to the event.
 
@@ -439,11 +515,11 @@ AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Double_t eta, Double_t
     e = nPart.E();
   }
 
-  return AddCluster(e, absId);
+  return AddCluster(e, absId, label);
 }
       
 //________________________________________________________________________
-AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId)
+AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId, Int_t label)
 {
   // Add a cluster to the event.
 
@@ -479,21 +555,26 @@ AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId)
   cluster->SetEmcCpvDistance(-1);
 
   //MC label
+  if (label == 0)
+    label = fMarkMC;
+  else
+    label += fMCLabelShift;
+
   if (fEsdMode) {
     AliESDCaloCluster *esdClus = static_cast<AliESDCaloCluster*>(cluster);
-    TArrayI parents(1, &fMarkMC);
+    TArrayI parents(1, &label);
     esdClus->AddLabels(parents); 
   }
   else {
     AliAODCaloCluster *aodClus = static_cast<AliAODCaloCluster*>(cluster);
-    aodClus->SetLabel(&fMarkMC, 1); 
+    aodClus->SetLabel(&label, 1); 
   }
   
   return cluster;
 }
 
 //________________________________________________________________________
-AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi, Byte_t type, Double_t etaemc, Double_t phiemc, Bool_t ise)
+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)
 {
   // Add a track to the event.
 
@@ -506,11 +587,16 @@ AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t
   if (phi < 0) 
     phi = GetRandomPhi();
 
+  if (label == 0)
+    label = fMarkMC;
+  else
+    label += fMCLabelShift;
+
   AliPicoTrack *track = new ((*fOutTracks)[nTracks]) AliPicoTrack(pt, 
                                                                  eta, 
                                                                  phi, 
                                                                  1,
-                                                                 fMarkMC,    // MC flag!
+                                                                 label,
                                                                  type, 
                                                                  etaemc, 
                                                                  phiemc, 
@@ -520,9 +606,23 @@ AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t
 }
 
 //________________________________________________________________________
+AliAODMCParticle* AliJetModelBaseTask::AddMCParticle(AliAODMCParticle *part, Int_t origIndex)
+{
+  const Int_t nPart = fOutMCParticles->GetEntriesFast();
+
+  AliAODMCParticle *aodpart = new ((*fOutMCParticles)[nPart]) AliAODMCParticle(*part);
+
+  fOutMCParticlesMap->SetBinContent(origIndex + fMCLabelShift, nPart);
+
+  return aodpart;
+}
+
+//________________________________________________________________________
 void AliJetModelBaseTask::CopyClusters()
 {
   // Copy all the clusters in the new collection
+  if (!fClusters)
+    return;
 
   const Int_t nClusters = fClusters->GetEntriesFast();
   Int_t nCopiedClusters = 0;
@@ -550,6 +650,11 @@ void AliJetModelBaseTask::CopyClusters()
 //________________________________________________________________________
 void AliJetModelBaseTask::CopyTracks()
 {
+  // Copy all the tracks in the new collection
+
+  if (!fTracks)
+    return;
+
   const Int_t nTracks = fTracks->GetEntriesFast();
   Int_t nCopiedTracks = 0;
   for (Int_t i = 0; i < nTracks; ++i) {
@@ -562,6 +667,39 @@ void AliJetModelBaseTask::CopyTracks()
 }
 
 //________________________________________________________________________
+void AliJetModelBaseTask::CopyMCParticles()
+{
+  // Copy all the MC particles in the new collection
+
+  if (!fMCParticles)
+    return;
+
+  const Int_t nPart = fMCParticles->GetEntriesFast();
+  Int_t nCopiedPart = 0;
+  for (Int_t i = 0; i < nPart; ++i) {
+    AliAODMCParticle *part = static_cast<AliAODMCParticle*>(fMCParticles->At(i));
+    if (!part)
+      continue;
+    new ((*fOutMCParticles)[nCopiedPart]) AliAODMCParticle(*part);
+
+    nCopiedPart++;
+  }
+
+  if (!fMCParticlesMap)
+    return;
+
+  Int_t shift = 0;
+
+  for (Int_t i = 0; i < fMCParticlesMap->GetNbinsX()+2; i++) {
+    fOutMCParticlesMap->SetBinContent(i, fMCParticlesMap->GetBinContent(i));
+    if (fMCParticlesMap->GetBinContent(i) != 0)
+      shift = i;
+  }
+
+  fMCLabelShift = shift;
+}
+
+//________________________________________________________________________
 void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId)
 {
   // Get random cell.
index f7e9bb9..3bf7973 100644 (file)
@@ -4,10 +4,12 @@
 // $Id$
 
 class TClonesArray;
+class TH1I;
 class AliEMCALGeometry;
 class AliVCluster;
 class AliPicoTrack;
 class AliVCaloCells;
+class AliAODMCParticle;
 
 #include <TH1F.h>
 #include <TF1.h>
@@ -34,6 +36,7 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE {
   void                   SetTracksName(const char *n)          { fTracksName   = n;    }
   void                   SetClusName(const char *n)            { fCaloName     = n;    }
   void                   SetCellsName(const char *n)           { fCellsName    = n;    }
+  void                   SetMCParticlesName(const char *n)     { fMCParticlesName = n; }
   void                   SetSuffix(const char *s)              { fSuffix       = s;    }
 
   void                   SetGeometryName(const char *n)        { fGeomName     = n;    }
@@ -46,14 +49,16 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE {
  protected:
   Int_t                  SetNumberOfOutCells(Int_t n);                                          // set the number of cells
   Int_t                  AddCell(Double_t e = -1, Double_t eta = -999, Double_t phi = -1);      // add a cell; if values are -1 generate random parameters
-  Int_t                  AddCell(Double_t e, Int_t absId, Double_t time = 0);                   // add a cell with given energy, position and times
-  AliVCluster           *AddCluster(Double_t e = -1, Double_t eta = -999, Double_t phi = -1);   // add a cluster; if values are -1 generate random parameters
-  AliVCluster           *AddCluster(Double_t e, Int_t absId);                                   // add a cluster with given energy and position
+  Int_t                  AddCell(Double_t e, Int_t absId, Double_t time = 0, Int_t label=0);   // add a cell with given energy, position and times
+  AliVCluster           *AddCluster(Double_t e = -1, Double_t eta = -999, Double_t phi = -1, Int_t label=0); // add a cluster; if values are -1 generate random parameters
+  AliVCluster           *AddCluster(Double_t e, Int_t absId, Int_t label=0);                   // add a cluster with given energy and position
   AliPicoTrack          *AddTrack(Double_t pt = -1, Double_t eta = -999, Double_t phi = -1, 
-                                 Byte_t type=0, Double_t etaemc=0, Double_t phiemc=0, Bool_t ise=kFALSE);    // add a track; if values are -1 generate random parameters
+                                 Byte_t type=0, Double_t etaemc=0, Double_t phiemc=0, Bool_t ise=kFALSE, Int_t label=0); // add a track; if values are -1 generate random parameters
+  AliAODMCParticle      *AddMCParticle(AliAODMCParticle *part, Int_t origIndex);                // add a MC particle
   void                   CopyCells();
   void                   CopyClusters();
   void                   CopyTracks();
+  void                   CopyMCParticles();
   void                   GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId);             // generate a random cell in the calorimeter
   Double_t               GetRandomEta();                                                        // generate a random eta value in the given range
   Double_t               GetRandomPhi();                                                        // generate a random phi value in the given range
@@ -68,6 +73,8 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE {
   TString                fOutCaloName;            // name of output cluster collection
   TString                fCellsName;              // name of calo cells collection
   TString                fOutCellsName;           // name of output cells collection
+  TString                fMCParticlesName;        // name of MC particle collection
+  TString                fOutMCParticlesName;     // name of output MC particle collection
   TString                fSuffix;                 // suffix to add in the name of new collections
   Float_t                fEtaMin;                 // eta minimum value
   Float_t                fEtaMax;                 // eta maximum value
@@ -91,6 +98,11 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE {
   AliVCaloCells         *fCaloCells;              //!cells collection
   AliVCaloCells         *fOutCaloCells;           //!output cells collection
   Int_t                  fAddedCells;             //!number of added cells
+  TClonesArray          *fMCParticles;            //!MC particles collection
+  TH1I                  *fMCParticlesMap;          //!MC particles mapping
+  TClonesArray          *fOutMCParticles;         //!output MC particles collection
+  TH1I                  *fOutMCParticlesMap;      //!MC particles mapping
+  Int_t                  fMCLabelShift;           //!MC label shift
   Bool_t                 fEsdMode;                //!ESD/AOD mode
   TList                 *fOutput;                 //!output list for QA histograms
 
@@ -98,6 +110,6 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE {
   AliJetModelBaseTask(const AliJetModelBaseTask&);            // not implemented
   AliJetModelBaseTask &operator=(const AliJetModelBaseTask&); // not implemented
 
-  ClassDef(AliJetModelBaseTask, 6) // Jet modelling task
+  ClassDef(AliJetModelBaseTask, 7) // Jet modelling task
 };
 #endif
index d06d655..4194f08 100644 (file)
@@ -16,8 +16,8 @@
 #include "AliVTrack.h"
 #include "AliEmcalJet.h"
 #include "AliGenPythiaEventHeader.h"
-#include "AliLog.h"
 #include "AliMCEvent.h"
+#include "AliLog.h"
 #include "AliRhoParameter.h"
 
 ClassImp(AliJetResponseMaker)
@@ -267,7 +267,7 @@ void AliJetResponseMaker::UserCreateOutputObjects()
     fHistDeltaCorrPtvsMatchingLevel->GetXaxis()->SetTitle("Matching level");  
     fHistDeltaCorrPtvsMatchingLevel->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
     fHistDeltaCorrPtvsMatchingLevel->GetZaxis()->SetTitle("counts");
-    fOutput->Add(fHistDeltaCorrPtvsJet2Pt);
+    fOutput->Add(fHistDeltaCorrPtvsMatchingLevel);
   }
 
   fHistNonMatchedJets1PtArea = new TH2F("fHistNonMatchedJets1PtArea", "fHistNonMatchedJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
@@ -461,24 +461,24 @@ Bool_t AliJetResponseMaker::RetrieveEventObjects()
 
   if (fRho2)
     fRho2Val = fRho2->GetVal();
-  
-  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
-
-  if (!fPythiaHeader)
-    return kFALSE;
 
   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
+  
+  if (MCEvent())
+    fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 
-  Double_t pthard = fPythiaHeader->GetPtHard();
-
-  for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
-    if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
-      break;
+  if (fPythiaHeader) {
+    Double_t pthard = fPythiaHeader->GetPtHard();
+    
+    for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
+      if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
+       break;
+    }
+    
+    fNTrials = fPythiaHeader->Trials();
   }
 
-  fNTrials = fPythiaHeader->Trials();
-
   return kTRUE;
 }
 
@@ -620,10 +620,10 @@ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet
        AliWarning(Form("Could not find track %d!", iTrack));
        continue;
       }
-      Int_t MClabel = track->GetLabel();
+      Int_t MClabel = TMath::Abs(track->GetLabel());
       Int_t index = -1;
          
-      if (MClabel <= 0) {// this is not a MC particle; remove it completely
+      if (MClabel == 0) {// this is not a MC particle; remove it completely
        AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
        totalPt1 -= track->Pt();
        d1 -= track->Pt();
@@ -662,10 +662,10 @@ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet
          Int_t cellId = clus->GetCellAbsId(iCell);
          Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
 
-         Int_t MClabel = fCaloCells->GetCellMCLabel(cellId);
+         Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
          Int_t index = -1;
          
-         if (MClabel <= 0) {// this is not a MC particle; remove it completely
+         if (MClabel == 0) {// this is not a MC particle; remove it completely
            AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
            totalPt1 -= part.Pt() * cellFrac;
            d1 -= part.Pt() * cellFrac;
@@ -693,10 +693,10 @@ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet
        }
       }
       else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
-       Int_t MClabel = clus->GetLabel();
+       Int_t MClabel = TMath::Abs(clus->GetLabel());
        Int_t index = -1;
            
-       if (MClabel <= 0) {// this is not a MC particle; remove it completely
+       if (MClabel == 0) {// this is not a MC particle; remove it completely
          AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
          totalPt1 -= part.Pt();
          d1 -= part.Pt();
index 6ab421a..156c5fb 100644 (file)
@@ -62,6 +62,7 @@ AliAnalysisTaskSAQA::AliAnalysisTaskSAQA() :
 
   for (Int_t i = 0; i < 4; i++) {
     for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
+    fHistTrPhiEtaPtNegLab[i] = 0;
     fHistClusPhiEtaEnergy[i] = 0;
     fHistJetsPhiEtaPt[i] = 0;
     fHistJetsPtArea[i] = 0;
@@ -101,6 +102,7 @@ AliAnalysisTaskSAQA::AliAnalysisTaskSAQA(const char *name) :
 
   for (Int_t i = 0; i < 4; i++) {
     for (Int_t j = 0; j < 4; j++) fHistTrPhiEtaPt[i][j] = 0;
+    fHistTrPhiEtaPtNegLab[i] = 0;
     fHistClusPhiEtaEnergy[i] = 0;
     fHistJetsPhiEtaPt[i] = 0;
     fHistJetsPtArea[i] = 0;
@@ -143,6 +145,14 @@ void AliAnalysisTaskSAQA::UserCreateOutputObjects()
        fHistTrPhiEtaPt[i][j]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
        fOutput->Add(fHistTrPhiEtaPt[i][j]);
       }
+      if (!fParticleLevel) {
+       histname = Form("fHistTrPhiEtaPtNegLab_%d",i);
+       fHistTrPhiEtaPtNegLab[i] = new TH3F(histname,histname, 100, -1, 1, 201, 0, TMath::Pi() * 2.01, fNbins, fMinBinPt, fMaxBinPt);
+       fHistTrPhiEtaPtNegLab[i]->GetXaxis()->SetTitle("#eta");
+       fHistTrPhiEtaPtNegLab[i]->GetYaxis()->SetTitle("#phi");
+       fHistTrPhiEtaPtNegLab[i]->GetZaxis()->SetTitle("p_{T} (GeV/c)");
+       fOutput->Add(fHistTrPhiEtaPtNegLab[i]);
+      }
     }
 
     fHistTrEmcPhiEta = new TH2F("fHistTrEmcPhiEta","Phi-Eta emcal distribution of tracks", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
@@ -499,6 +509,9 @@ Float_t AliAnalysisTaskSAQA::DoTrackLoop()
     }
     else {
       fHistTrPhiEtaPt[fCentBin][3]->Fill(track->Eta(), track->Phi(), track->Pt());
+      if (fHistTrPhiEtaPtNegLab[fCentBin] && track->GetLabel() < 0)
+       fHistTrPhiEtaPtNegLab[fCentBin]->Fill(track->Eta(), track->Phi(), track->Pt());
+
       Int_t type = 0;
 
       AliPicoTrack* ptrack = dynamic_cast<AliPicoTrack*>(track);
index 54c07a2..72ce435 100644 (file)
@@ -50,6 +50,7 @@ class AliAnalysisTaskSAQA : public AliAnalysisTaskEmcalJet {
 
   // Tracks
   TH3F                       *fHistTrPhiEtaPt[4][4];     //!Phi-Eta-Pt distribution of tracks
+  TH3F                       *fHistTrPhiEtaPtNegLab[4];  //!Phi-Eta-Pt distribution of tracks with negative labels
   TH2F                       *fHistTrEmcPhiEta;          //!Phi-Eta emcal propagated distribution of tracks
   TH2F                       *fHistTrPhiEtaNonProp;      //!Phi-Eta distribution of non emcal propagated tracks
   TH2F                       *fHistDeltaEtaPt;           //!Eta-EtaProp vs. Pt
index 3be0ea9..562aeaa 100644 (file)
@@ -4,17 +4,22 @@ TObjArray* GenerateFileList(const char* list, Int_t nFiles);
 
 AliJetEmbeddingFromAODTask* AddTaskJetEmbeddingFromAOD(
   const char     *tracksName    = "Tracks",
+  const char     *clusName      = "",
   const char     *cellsName     = "EMCALCells",
+  const char     *MCPartName    = "",
   const char     *fileList      = "files.txt",
   const char     *aodTreeName   = "aodTree",
   const char     *aodTracksName = "tracks",
+  const char     *aodClusName   = "",
   const char     *aodCellsName  = "emcalCells",
+  const char     *aodMCPartName = "",
   const char     *runperiod     = "",
   Bool_t          includeNoITS  = kTRUE,
   Double_t        minCent       = 0,
   Double_t        maxCent       = 10,
   UInt_t          mask          = AliVEvent::kAny,
   const Int_t     nTracks       = 1234567890,
+  const Int_t     nClus         = 0,
   const Int_t     nCells        = 1234567890,
   const Bool_t    copyArray     = kTRUE,
   const Int_t     nFiles        = 1234567890,
@@ -51,14 +56,19 @@ AliJetEmbeddingFromAODTask* AddTaskJetEmbeddingFromAOD(
 
   AliJetEmbeddingFromAODTask *jetEmb = new AliJetEmbeddingFromAODTask(taskName,makeQA);
   jetEmb->SetTracksName(tracksName);
+  jetEmb->SetClusName(clusName);
   jetEmb->SetCellsName(cellsName);
+  jetEmb->SetMCParticlesName(MCPartName);
   jetEmb->SetFileList(GenerateFileList(fileList, nFiles));
   jetEmb->SetAODTreeName(aodTreeName);
   jetEmb->SetAODTracksName(aodTracksName);
+  jetEmb->SetAODClusName(aodClusName);
   jetEmb->SetAODCellsName(aodCellsName);
+  jetEmb->SetAODMCParticlesName(aodMCPartName);
   jetEmb->SetCentralityRange(minCent, maxCent);
   jetEmb->SetTriggerMask(mask);
   jetEmb->SetNCells(nCells);
+  jetEmb->SetNClusters(nClus);
   jetEmb->SetNTracks(nTracks);
   jetEmb->SetCopyArray(copyArray);
   jetEmb->SetEtaRange(minEta, maxEta);