Update from Salvatore for full jet embedding
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 26 Jan 2013 23:11:01 +0000 (23:11 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 26 Jan 2013 23:11:01 +0000 (23:11 +0000)
14 files changed:
PWGJE/CMakelibPWGJEEMCALJetTasks.pkg
PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx [new file with mode: 0644]
PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.h [new file with mode: 0644]
PWGJE/EMCALJetTasks/AliJetEmbeddingTask.cxx
PWGJE/EMCALJetTasks/AliJetEmbeddingTask.h
PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx
PWGJE/EMCALJetTasks/AliJetModelBaseTask.h
PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx
PWGJE/EMCALJetTasks/AliJetResponseMaker.h
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.cxx
PWGJE/EMCALJetTasks/macros/AddTaskJetEmbeddingFromAOD.C [new file with mode: 0644]
PWGJE/EMCALJetTasks/macros/AddTaskJetRespPtHard.C
PWGJE/EMCALJetTasks/macros/AddTaskJetResponseMaker.C
PWGJE/PWGJEEMCALJetTasksLinkDef.h

index bc4516d..3da6d6c 100644 (file)
@@ -39,6 +39,7 @@ set ( SRCS
  EMCALJetTasks/AliHadCorrTask.cxx
  EMCALJetTasks/AliJetEmbeddingTask.cxx
  EMCALJetTasks/AliJetEmbeddingFromGenTask.cxx
+ EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx
  EMCALJetTasks/AliJetModelBaseTask.cxx
  EMCALJetTasks/AliJetRandomizerTask.cxx
  EMCALJetTasks/AliJetResponseMaker.cxx
diff --git a/PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx b/PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx
new file mode 100644 (file)
index 0000000..c1c86ba
--- /dev/null
@@ -0,0 +1,380 @@
+// $Id: AliJetEmbeddingFromAODTask.cxx $
+//
+// Jet embedding from AOD task.
+//
+// Author: S.Aiola, C.Loizides
+
+#include "AliJetEmbeddingFromAODTask.h"
+
+//#include <iostream>
+
+#include <TFile.h>
+#include <TTree.h>
+#include <TClonesArray.h>
+#include <TObjArray.h>
+#include <TObjString.h>
+#include <TGrid.h>
+#include <TH2C.h>
+
+#include "AliVEvent.h"
+#include "AliAODTrack.h"
+#include "AliESDtrack.h"
+#include "AliPicoTrack.h"
+#include "AliVTrack.h"
+#include "AliVCluster.h"
+#include "AliVCaloCells.h"
+#include "AliCentrality.h"
+#include "AliVHeader.h"
+#include "AliVVertex.h"
+#include "AliAODHeader.h"
+#include "AliLog.h"
+#include "AliInputEventHandler.h"
+
+ClassImp(AliJetEmbeddingFromAODTask)
+
+//________________________________________________________________________
+AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() : 
+  AliJetModelBaseTask("AliJetEmbeddingFromAODTask"),
+  fFileList(0),
+  fAODTreeName(),
+  fAODHeaderName(),
+  fAODVertexName(),
+  fAODTrackName(),
+  fAODClusName(),
+  fAODCellsName(),
+  fMinCentrality(0),
+  fMaxCentrality(10),
+  fTriggerMask(AliVEvent::kAny),
+  fZVertexCut(10),
+  fIncludeNoITS(kTRUE),
+  fTotalFiles(2200),
+  fEsdTreeMode(kFALSE),
+  fCurrentFileID(0),
+  fCurrentAODFileID(-1),
+  fCurrentAODFile(0),
+  fAODHeader(0),
+  fAODVertex(0),
+  fAODTracks(0),
+  fAODClusters(0),
+  fAODCaloCells(0),
+  fHistFileIDs(0)
+{
+  // Default constructor.
+  SetSuffix("AODEmbedding");
+  fMarkMC = kFALSE;
+  fAODfilterBits[0] = -1;
+  fAODfilterBits[1] = -1;
+}
+
+//________________________________________________________________________
+AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t drawqa) : 
+  AliJetModelBaseTask(name, drawqa),
+  fFileList(0),
+  fAODTreeName("aodTree"),
+  fAODHeaderName("header"),
+  fAODVertexName("vertices"),
+  fAODTrackName("tracks"),
+  fAODClusName("caloClusters"),
+  fAODCellsName("emcalCells"),
+  fMinCentrality(0),
+  fMaxCentrality(10),
+  fTriggerMask(AliVEvent::kAny),
+  fZVertexCut(10),
+  fIncludeNoITS(kTRUE),
+  fTotalFiles(2200),
+  fEsdTreeMode(kFALSE),
+  fCurrentFileID(0),
+  fCurrentAODFileID(-1),
+  fCurrentAODFile(0),
+  fAODHeader(0),
+  fAODVertex(0),
+  fAODTracks(0),
+  fAODClusters(0),
+  fAODCaloCells(0),
+  fHistFileIDs(0)
+{
+  // Standard constructor.
+  SetSuffix("AODEmbedding");
+  fMarkMC = kFALSE;
+  fAODfilterBits[0] = -1;
+  fAODfilterBits[1] = -1;
+}
+
+//________________________________________________________________________
+AliJetEmbeddingFromAODTask::~AliJetEmbeddingFromAODTask()
+{
+  // Destructor
+
+  if (fCurrentAODFile) {
+    fCurrentAODFile->Close();
+    delete fCurrentAODFile;
+  }
+}
+
+//________________________________________________________________________
+void AliJetEmbeddingFromAODTask::UserCreateOutputObjects()
+{
+  if (!fQAhistos)
+    return;
+
+  AliJetModelBaseTask::UserCreateOutputObjects();
+  
+  fHistFileIDs = new TH2C("fHistFileIDs", "fHistFileIDs", fTotalFiles, -0.5, fTotalFiles - 0.5, fFileList->GetEntriesFast(), -0.5, fFileList->GetEntriesFast() -0.5);
+  fHistFileIDs->GetXaxis()->SetTitle("File no. (PYTHIA)");
+  fHistFileIDs->GetYaxis()->SetTitle("File no. (Embedded AOD)");
+  fOutput->Add(fHistFileIDs);
+
+  PostData(1, fOutput);
+}
+
+//________________________________________________________________________
+Bool_t AliJetEmbeddingFromAODTask::UserNotify()
+{
+  TString path(fInputHandler->GetTree()->GetCurrentFile()->GetPath());
+  if (path.EndsWith("/"))
+    path.Remove(path.Length()-1);
+  path.Remove(path.Last('/'));
+  path.Remove(0, path.Last('/')+1);
+  if (!path.IsDec()) {
+    AliWarning(Form("Could not extract file number from path %s", path.Data()));
+    return kTRUE;
+  }
+  fCurrentFileID = path.Atoi();
+  fCurrentAODFileID = fFileList->GetEntriesFast() * fCurrentFileID / fTotalFiles;
+  AliInfo(Form("Start embedding from file ID %d", fCurrentAODFileID));
+  return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliJetEmbeddingFromAODTask::ExecOnce() 
+{
+  if (fAODTreeName.Contains("aod"))
+    fEsdTreeMode = kFALSE;
+  else
+    fEsdTreeMode = kTRUE;
+
+  return AliJetModelBaseTask::ExecOnce();
+}
+
+//________________________________________________________________________
+Bool_t AliJetEmbeddingFromAODTask::OpenNextFile() 
+{
+  if (fCurrentAODFile) {
+    fCurrentAODFile->Close();
+    delete fCurrentAODFile;
+    fCurrentAODFile = 0;
+  }
+
+  if (fCurrentAODFileID >= fFileList->GetEntriesFast())
+    return kFALSE;
+
+  TObjString *objFileName = static_cast<TObjString*>(fFileList->At(fCurrentAODFileID));
+  TString fileName = objFileName->GetString();
+  
+  if (fileName.BeginsWith("alien://") && !gGrid) {
+      AliInfo("Trying to connect to AliEn ...");
+      TGrid::Connect("alien://");
+  }
+
+  fCurrentAODFile = TFile::Open(fileName);
+
+  if (!fCurrentAODFile || fCurrentAODFile->IsZombie()) {
+    return kFALSE;
+  }
+
+  if (fQAhistos)
+    fHistFileIDs->Fill(fCurrentFileID, fCurrentAODFileID);
+  
+  fCurrentAODFileID++;
+  return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliJetEmbeddingFromAODTask::GetNextEntry() 
+{
+  static TTree *tree = 0;
+  static Int_t entry = 0;
+
+  Int_t attempts = 0;
+
+  do {
+    
+    if (!fCurrentAODFile || !tree || entry >= tree->GetEntries()) {
+      if (!OpenNextFile())
+       return kFALSE;
+      
+      tree = static_cast<TTree*>(fCurrentAODFile->Get(fAODTreeName));
+      if (!tree)
+       return kFALSE;
+
+      if (!fAODHeaderName.IsNull()) 
+       tree->SetBranchAddress(fAODHeaderName, &fAODHeader);
+
+      if (!fAODVertexName.IsNull()) 
+       tree->SetBranchAddress(fAODVertexName, &fAODVertex);
+      
+      if (!fAODTrackName.IsNull()) 
+       tree->SetBranchAddress(fAODTrackName, &fAODTracks);
+      
+      if (!fAODClusName.IsNull()) 
+       tree->SetBranchAddress(fAODClusName, &fAODClusters);
+      
+      if (!fAODCellsName.IsNull()) 
+       tree->SetBranchAddress(fAODCellsName, &fAODCaloCells);
+      
+      entry = 0;
+    }
+    
+    tree->GetEntry(entry);
+    entry++;
+    attempts++;
+
+    if (attempts == 1000) 
+      AliWarning("After 1000 attempts no event has been accepted by the event selection (trigger, centrality...)!");
+
+  } while (!IsAODEventSelected() && tree);
+
+  return (tree!=0);
+}
+
+//________________________________________________________________________
+Bool_t AliJetEmbeddingFromAODTask::IsAODEventSelected()
+{
+  if (!fEsdTreeMode && fAODHeader) {
+    AliAODHeader *aodHeader = static_cast<AliAODHeader*>(fAODHeader);
+
+    UInt_t offlineTrigger = aodHeader->GetOfflineTrigger();
+  
+    if ((offlineTrigger & fTriggerMask) == 0) {
+      AliDebug(2, Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.", offlineTrigger, fTriggerMask));
+      return kFALSE;
+    }
+    
+    AliCentrality *cent = aodHeader->GetCentralityP();
+    Float_t centVal = cent->GetCentralityPercentile("V0M");
+    if (centVal < fMinCentrality || centVal >= fMaxCentrality) {
+      AliDebug(2, Form("Event rejected due to centrality selection. Event centrality: %f, centrality range selection: %f to %f", centVal, fMinCentrality, fMaxCentrality));
+      return kFALSE;
+    }
+  }
+
+  if (fAODVertex) {
+    Double_t vert[3]={0};
+    ((AliVVertex*)fAODVertex->At(0))->GetXYZ(vert);
+    if (TMath::Abs(vert[2]) > fZVertexCut) {
+      AliDebug(2,Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f", vert[2], fZVertexCut));
+      return kFALSE;
+    }
+  }
+
+  return kTRUE;
+}
+
+//________________________________________________________________________
+void AliJetEmbeddingFromAODTask::Run() 
+{
+  if (!GetNextEntry()) {
+    AliError("Unable to get AOD event to embed. Nothing will be embedded.");
+    return;
+  }
+
+  if (fTracks && fOutTracks) {
+
+    if (fCopyArray)
+      CopyTracks();
+
+    if (fAODTracks) {
+      AliDebug(2, Form("%d tracks will be processed for embedding.", fAODTracks->GetEntriesFast()));
+      for (Int_t i = 0; i < fAODTracks->GetEntriesFast(); i++) {
+       AliVTrack *track = static_cast<AliVTrack*>(fAODTracks->At(i));
+       if (!track) {
+         AliError(Form("Could not find track %d in branch %s of tree %s!", i, fAODTrackName.Data(), fAODTreeName.Data()));
+         continue;
+       }
+       
+       Int_t type = 0;
+       Bool_t isEmc = kFALSE;
+       if (!fEsdTreeMode) {
+         AliAODTrack *aodtrack = static_cast<AliAODTrack*>(track);
+         if (aodtrack->TestFilterBit(fAODfilterBits[0])) {
+           type = 0;
+         }
+         else if (aodtrack->TestFilterBit(fAODfilterBits[1])) {
+           if ((aodtrack->GetStatus()&AliESDtrack::kITSrefit)==0) {
+             if (fIncludeNoITS)
+               type = 2;
+             else
+               continue;
+           }
+           else {
+             type = 1;
+           }
+         }
+         else { /*not a good track*/
+           continue;
+         }
+
+         if (TMath::Abs(aodtrack->GetTrackEtaOnEMCal()) < 0.75 && 
+             aodtrack->GetTrackPhiOnEMCal() > 70 * TMath::DegToRad() && 
+             aodtrack->GetTrackPhiOnEMCal() < 190 * TMath::DegToRad())
+           isEmc = kTRUE;
+       }
+       else { /*not AOD mode, let's see if it is PicoTrack*/
+         AliPicoTrack *ptrack = dynamic_cast<AliPicoTrack*>(track);
+         if (ptrack) {
+           type = ptrack->GetTrackType();
+           isEmc = ptrack->IsEMCAL();
+         }
+       }
+       
+       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);
+      }
+    }
+  }
+
+  if (fClusters && fOutClusters) {
+
+    if (fCopyArray)
+      CopyClusters();
+
+    if (fAODClusters) {
+      for (Int_t i = 0; i < fAODClusters->GetEntriesFast(); i++) {
+       AliVCluster *clus = static_cast<AliVCluster*>(fAODClusters->At(i));
+       if (!clus) {
+         AliError(Form("Could not find cluster %d in branch %s of tree %s!", i, fAODClusName.Data(), fAODTreeName.Data()));
+         continue;
+       }
+       TLorentzVector vect;
+       Double_t vert[3] = {0,0,0};
+       clus->GetMomentum(vect,vert);
+       AddCluster(clus->E(), vect.Eta(), vect.Phi());
+      }
+    }
+  }
+
+  if (fCaloCells && fOutCaloCells) {
+
+    Int_t totalCells = fCaloCells->GetNumberOfCells();
+    if (fAODCaloCells)
+      totalCells += fAODCaloCells->GetNumberOfCells();
+
+    SetNumberOfOutCells(totalCells);
+
+    CopyCells();
+
+    if (fAODCaloCells) {
+      for (Short_t i = 0; i < fAODCaloCells->GetNumberOfCells(); i++) {
+       Short_t mclabel = -1;
+       Double_t efrac = 0.;
+       Double_t time = -1;
+       Short_t cellNum = -1;
+       Double_t amp = -1;
+       
+       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);
+      }
+    }
+  }
+}
diff --git a/PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.h b/PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.h
new file mode 100644 (file)
index 0000000..b29db54
--- /dev/null
@@ -0,0 +1,75 @@
+#ifndef ALIJETEMBEDDINGFROMAODTASK_H
+#define ALIJETEMBEDDINGFROMAODTASK_H
+
+// $Id: AliJetEmbeddingFromAODTask.h  $
+
+class TFile;
+class TObjArray;
+class TClonesArray;
+class TString;
+class AliVCaloCells;
+class AliVHeader;
+class TH2;
+
+#include "AliJetModelBaseTask.h"
+
+class AliJetEmbeddingFromAODTask : public AliJetModelBaseTask {
+ public:
+  AliJetEmbeddingFromAODTask();
+  AliJetEmbeddingFromAODTask(const char *name, Bool_t drawqa=kFALSE); 
+  virtual ~AliJetEmbeddingFromAODTask();
+
+  void           UserCreateOutputObjects();
+  Bool_t         UserNotify();
+
+  void           SetFileList(TObjArray *list)                      { fFileList           = list  ; }
+  void           SetAODTreeName(const char *t)                     { fAODTreeName        = t     ; }
+  void           SetAODHeaderName(const char *t)                   { fAODHeaderName      = t     ; }
+  void           SetAODTracksName(const char *n)                   { fAODTrackName       = n     ; }
+  void           SetAODClusName(const char *n)                     { fAODClusName        = n     ; }
+  void           SetAODCellsName(const char *n)                    { fAODCellsName       = 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  ; }
+  void           SetIncludeNoITS(Bool_t f)                         { fIncludeNoITS       = f     ; }
+  void           SetTotalFiles(Int_t n)                            { fTotalFiles         = n     ; }
+
+ protected:
+  Bool_t         ExecOnce()            ;// intialize task
+  void           Run()                 ;// do jet model action
+  Bool_t         OpenNextFile()        ;// open next file in fFileList
+  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
+  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)
+
+ private:
+  AliJetEmbeddingFromAODTask(const AliJetEmbeddingFromAODTask&);            // not implemented
+  AliJetEmbeddingFromAODTask &operator=(const AliJetEmbeddingFromAODTask&); // not implemented
+
+  ClassDef(AliJetEmbeddingFromAODTask, 1) // Jet embedding from AOD task
+};
+#endif
index 5fa78e2..59dc705 100644 (file)
@@ -6,19 +6,6 @@
 
 #include "AliJetEmbeddingTask.h"
 
-#include <TClonesArray.h>
-#include <TLorentzVector.h>
-#include <TRandom3.h>
-
-#include "AliAnalysisManager.h"
-#include "AliVEvent.h"
-#include "AliVCluster.h"
-#include "AliEMCALDigit.h"
-#include "AliEMCALRecPoint.h"
-#include "AliPicoTrack.h"
-#include "AliEMCALGeometry.h"
-#include "AliLog.h"
-
 ClassImp(AliJetEmbeddingTask)
 
 //________________________________________________________________________
index ab4f08d..b065d62 100644 (file)
@@ -3,9 +3,6 @@
 
 // $Id$
 
-class TClonesArray;
-class AliEMCALGeometry;
-
 #include "AliJetModelBaseTask.h"
 
 class AliJetEmbeddingTask : public AliJetModelBaseTask {
index af91d1b..d8b2897 100644 (file)
 #include <TH1.h>
 #include <TLorentzVector.h>
 #include <TRandom3.h>
+#include <TList.h>
 
+#include "AliVEvent.h"
 #include "AliAODCaloCluster.h"
-#include "AliAnalysisManager.h"
+#include "AliESDCaloCluster.h"
+#include "AliVCluster.h"
 #include "AliEMCALDigit.h"
-#include "AliEMCALGeometry.h"
 #include "AliEMCALRecPoint.h"
-#include "AliESDCaloCluster.h"
-#include "AliLog.h"
+#include "AliESDCaloCells.h"
+#include "AliAODCaloCells.h"
+#include "AliVCaloCells.h"
 #include "AliPicoTrack.h"
-#include "AliVCluster.h"
-#include "AliVEvent.h"
+#include "AliEMCALGeometry.h"
+#include "AliLog.h"
 
 ClassImp(AliJetModelBaseTask)
 
@@ -32,6 +35,8 @@ AliJetModelBaseTask::AliJetModelBaseTask() :
   fOutTracksName(),
   fCaloName(),
   fOutCaloName(),
+  fCellsName(),
+  fOutCellsName(),
   fSuffix(),
   fEtaMin(-1),
   fEtaMax(1),
@@ -41,27 +46,36 @@ AliJetModelBaseTask::AliJetModelBaseTask() :
   fPtMax(0),
   fCopyArray(kTRUE),
   fNClusters(0),
+  fNCells(0),
   fNTracks(0),
   fMarkMC(kTRUE),
   fPtSpectrum(0),
+  fQAhistos(kFALSE),
   fIsInit(0),
   fGeom(0),
   fClusters(0),
   fOutClusters(0),
   fTracks(0),
-  fOutTracks(0)
+  fOutTracks(0),
+  fCaloCells(0),
+  fOutCaloCells(0),
+  fAddedCells(0),
+  fEsdMode(kFALSE),
+  fOutput(0)
 {
   // Default constructor.
 }
 
 //________________________________________________________________________
-AliJetModelBaseTask::AliJetModelBaseTask(const char *name) : 
+AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) : 
   AliAnalysisTaskSE(name),
   fGeomName(""),
   fTracksName("PicoTracks"),
   fOutTracksName("PicoTracksEmbedded"),
   fCaloName("CaloClustersCorr"),
   fOutCaloName("CaloClustersCorrEmbedded"),
+  fCellsName(""),
+  fOutCellsName(""),
   fSuffix("Processed"),
   fEtaMin(-1),
   fEtaMax(1),
@@ -71,17 +85,28 @@ AliJetModelBaseTask::AliJetModelBaseTask(const char *name) :
   fPtMax(60),
   fCopyArray(kTRUE),
   fNClusters(0),
+  fNCells(0),
   fNTracks(1),
   fMarkMC(kTRUE),
   fPtSpectrum(0),
+  fQAhistos(drawqa),
   fIsInit(0),
   fGeom(0),
   fClusters(0),
   fOutClusters(0),
   fTracks(0),
-  fOutTracks(0)
+  fOutTracks(0),
+  fCaloCells(0),
+  fOutCaloCells(0),
+  fAddedCells(0),
+  fEsdMode(kFALSE),
+  fOutput(0)
 {
   // Standard constructor.
+
+  if (fQAhistos) {
+    DefineOutput(1, TList::Class()); 
+  }
 }
 
 //________________________________________________________________________
@@ -91,6 +116,20 @@ AliJetModelBaseTask::~AliJetModelBaseTask()
 }
 
 //________________________________________________________________________
+void AliJetModelBaseTask::UserCreateOutputObjects()
+{
+  // Create user output.
+  if (!fQAhistos)
+    return;
+
+  OpenFile(1);
+  fOutput = new TList();
+  fOutput->SetOwner();
+
+  PostData(1, fOutput);
+}
+
+//________________________________________________________________________
 void AliJetModelBaseTask::UserExec(Option_t *) 
 {
   // Execute per event.
@@ -108,10 +147,257 @@ void AliJetModelBaseTask::UserExec(Option_t *)
       fOutClusters->Delete();
   }
 
+  if (fCaloCells) {
+    fAddedCells = 0;
+    if (!fCopyArray)
+      fCaloCells = static_cast<AliVCaloCells*>(fCaloCells->Clone(Form("%s_old",fCaloCells->GetName())));
+  }
+
   Run();
 }
 
 //________________________________________________________________________
+Bool_t AliJetModelBaseTask::ExecOnce()
+{
+  // Init task.
+
+  delete gRandom;
+  gRandom = new TRandom3(0);
+
+  fEsdMode = InputEvent()->InheritsFrom("AliESDEvent");
+
+  if (fPtMax < fPtMin) {
+    AliWarning (Form("PtMax (%f) < PtMin (%f), setting PtMax = PtMin = %f", fPtMax, fPtMin, fPtMin));
+    fPtMax = fPtMin;
+  }
+
+  if (fEtaMax < fEtaMin) {
+    AliWarning (Form("EtaMax (%f) < EtaMin (%f), setting EtaMax = EtaMin = %f", fEtaMax, fEtaMin, fEtaMin));
+    fEtaMax = fEtaMin;
+  }
+
+  if (fPhiMax < fPhiMin) {
+    AliWarning (Form("PhiMax (%f) < PhiMin (%f), setting PhiMax = PhiMin = %f", fPhiMax, fPhiMin, fPhiMin));
+    fPhiMax = fPhiMin;
+  }
+
+  if (!fCaloCells && !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;
+    }
+    
+    if (!fCaloCells->InheritsFrom("AliVCaloCells")) {
+      AliError(Form("%s: Collection %s does not contain a AliVCaloCells object!", GetName(), fCellsName.Data())); 
+      fCaloCells = 0;
+      return kFALSE;
+    }
+
+    if (!fOutCaloCells) {
+      fOutCellsName = fCellsName;
+      if (fCopyArray) {
+       fOutCellsName += fSuffix;
+
+       if (fEsdMode) 
+         fOutCaloCells = new AliESDCaloCells(fOutCellsName,fOutCellsName);
+       else
+         fOutCaloCells = new AliAODCaloCells(fOutCellsName,fOutCellsName);
+
+       if (InputEvent()->FindListObject(fOutCellsName)) {
+         AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCellsName.Data())); 
+         return kFALSE;
+       }
+       else {
+         InputEvent()->AddObject(fOutCaloCells);
+       }
+      }
+      else {
+       fOutCaloCells = fCaloCells;
+      }
+    }
+  }
+
+  if (fNTracks > 0 && !fTracks && !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;
+    }
+    
+    if (!fTracks->GetClass()->GetBaseClass("AliPicoTrack")) {
+      AliError(Form("%s: Collection %s does not contain AliPicoTrack objects!", GetName(), fTracksName.Data())); 
+      fTracks = 0;
+      return kFALSE;
+    }
+
+    if (!fOutTracks) {
+      fOutTracksName = fTracksName;
+      if (fCopyArray) {
+       fOutTracksName += fSuffix;
+       fOutTracks = new TClonesArray("AliPicoTrack", fTracks->GetSize());
+       fOutTracks->SetName(fOutTracksName);
+       if (InputEvent()->FindListObject(fOutTracksName)) {
+         AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutTracksName.Data())); 
+         return kFALSE;
+       }
+       else {
+         InputEvent()->AddObject(fOutTracks);
+       }
+      }
+      else {
+       fOutTracks = fTracks;
+      }
+    }
+  }
+
+  if (fNClusters > 0 && !fClusters && !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;
+    }
+
+    if (!fClusters->GetClass()->GetBaseClass("AliVCluster")) {
+      AliError(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloName.Data())); 
+      fClusters = 0;
+      return kFALSE;
+    }
+
+    if (!fOutClusters) {
+      fOutCaloName = fCaloName;
+      if (fCopyArray) {
+       fOutCaloName += fSuffix;
+       fOutClusters = new TClonesArray(fClusters->GetClass()->GetName(), fClusters->GetSize());
+       fOutClusters->SetName(fOutCaloName);
+       if (InputEvent()->FindListObject(fOutCaloName)) {
+         AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCaloName.Data())); 
+         return kFALSE;
+       }
+       else {
+         InputEvent()->AddObject(fOutClusters);
+       }
+      }
+      else {
+       fOutClusters = fClusters;
+      }
+    }
+  }
+
+if (!fGeom && (fClusters || fCaloCells)) {
+    if (fGeomName.Length() > 0) {
+      fGeom = AliEMCALGeometry::GetInstance(fGeomName);
+      if (!fGeom)
+       AliError(Form("Could not get geometry with name %s!", fGeomName.Data()));
+    } else {
+      fGeom = AliEMCALGeometry::GetInstance();
+      if (!fGeom) 
+       AliError("Could not get default geometry!");
+    }
+  }
+  
+  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;
+}
+
+//________________________________________________________________________
+Int_t AliJetModelBaseTask::SetNumberOfOutCells(Int_t n)
+{
+  if (fOutCaloCells->GetNumberOfCells() < n) {
+    fOutCaloCells->DeleteContainer();
+    fOutCaloCells->CreateContainer(n);
+  }
+  else {
+    fOutCaloCells->SetNumberOfCells(n);
+  }
+
+  fAddedCells = 0;
+
+  return n;
+}
+
+//________________________________________________________________________
+void AliJetModelBaseTask::CopyCells()
+{
+  for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) {
+    Short_t mclabel = -1;
+    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();
+}
+
+//________________________________________________________________________
+Int_t AliJetModelBaseTask::AddCell(Double_t e, Double_t eta, Double_t phi)
+{
+  // Add a cell to the event.
+
+  Int_t absId = 0;
+  if (eta < -100 || phi < 0) {
+    GetRandomCell(eta, phi, absId);
+  }
+  else {
+    fGeom->EtaPhiFromIndex(absId, eta, phi);
+  }
+
+  if (absId == -1) {
+    AliWarning(Form("Unable to embed cell in eta = %f, phi = %f!"
+                   " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])", 
+                   eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax));
+    return 0;
+  } 
+
+  if (e < 0) {
+    Double_t pt = GetRandomPt();
+    TLorentzVector nPart;
+    nPart.SetPtEtaPhiM(pt, eta, phi, 0);
+    e = nPart.E();
+  }
+
+  return AddCell(e, absId);
+}
+
+//________________________________________________________________________
+Int_t AliJetModelBaseTask::AddCell(Double_t e, Int_t absId, Double_t time)
+{
+  // Add a cell to the event.
+
+  Int_t label = fMarkMC ? 100 : -1;
+
+  Bool_t r = fOutCaloCells->SetCell(fAddedCells, absId, e, time, label, 0);
+
+  if (r) {
+    fAddedCells++;
+    return fAddedCells;
+  }
+  else {
+    return 0;
+  }
+}
+
+
+//________________________________________________________________________
 AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Double_t eta, Double_t phi)
 {
   // Add a cluster to the event.
@@ -183,7 +469,7 @@ AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId)
 }
 
 //________________________________________________________________________
-AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi)
+AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi, Byte_t type, Double_t etaemc, Double_t phiemc, Bool_t ise)
 {
   // Add a track to the event.
 
@@ -196,16 +482,18 @@ AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t
   if (phi < 0) 
     phi = GetRandomPhi();
 
-  Int_t label = fMarkMC ? 100 : 0;
+  Int_t label = fMarkMC ? 100 : -1;
 
   AliPicoTrack *track = new ((*fOutTracks)[nTracks]) AliPicoTrack(pt, 
                                                                  eta, 
                                                                  phi, 
-                                                                 1, 
-                                                                 label,    // MC flag!      
-                                                                 0, 
-                                                                 0, 
-                                                                 kFALSE);
+                                                                 1,
+                                                                 label,    // MC flag!
+                                                                 type, 
+                                                                 etaemc, 
+                                                                 phiemc, 
+                                                                 ise);
+
   return track;
 }
 
@@ -214,11 +502,10 @@ void AliJetModelBaseTask::CopyClusters()
 {
   // Copy all the clusters in the new collection
 
-  Bool_t esdMode = (Bool_t)(fClusters->GetClass()->GetBaseClass("AliESDCaloCluster") != 0);
   const Int_t nClusters = fClusters->GetEntriesFast();
   Int_t nCopiedClusters = 0;
   
-  if (esdMode) {
+  if (fEsdMode) {
     for (Int_t i = 0; i < nClusters; ++i) {
       AliESDCaloCluster *esdcluster = static_cast<AliESDCaloCluster*>(fClusters->At(i));
       if (!esdcluster || !esdcluster->IsEMCAL())
@@ -253,126 +540,6 @@ void AliJetModelBaseTask::CopyTracks()
 }
 
 //________________________________________________________________________
-Bool_t AliJetModelBaseTask::ExecOnce()
-{
-  // Init task.
-
-  delete gRandom;
-  gRandom = new TRandom3(0);
-
-  if (fPtMax < fPtMin) {
-    AliWarning (Form("PtMax (%f) < PtMin (%f), setting PtMax = PtMin = %f", fPtMax, fPtMin, fPtMin));
-    fPtMax = fPtMin;
-  }
-
-  if (fEtaMax < fEtaMin) {
-    AliWarning (Form("EtaMax (%f) < EtaMin (%f), setting EtaMax = EtaMin = %f", fEtaMax, fEtaMin, fEtaMin));
-    fEtaMax = fEtaMin;
-  }
-
-  if (fPhiMax < fPhiMin) {
-    AliWarning (Form("PhiMax (%f) < PhiMin (%f), setting PhiMax = PhiMin = %f", fPhiMax, fPhiMin, fPhiMin));
-    fPhiMax = fPhiMin;
-  }
-
-  if (fNTracks > 0 && !fTracks && !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;
-    }
-    
-    if (!fTracks->GetClass()->GetBaseClass("AliPicoTrack")) {
-      AliError(Form("%s: Collection %s does not contain AliPicoTrack objects!", GetName(), fTracksName.Data())); 
-      fTracks = 0;
-      return kFALSE;
-    }
-
-    if (!fOutTracks) {
-      fOutTracksName = fTracksName;
-      if (fCopyArray) {
-       fOutTracksName += fSuffix;
-       fOutTracks = new TClonesArray("AliPicoTrack", fTracks->GetSize());
-       fOutTracks->SetName(fOutTracksName);
-       if (InputEvent()->FindListObject(fOutTracksName)) {
-         AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutTracksName.Data())); 
-         return kFALSE;
-       }
-       else {
-         InputEvent()->AddObject(fOutTracks);
-       }
-      }
-      else {
-       fOutTracks = fTracks;
-      }
-    }
-  }
-
-  if (fNClusters > 0 && !fClusters && !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;
-    }
-
-    if (!fClusters->GetClass()->GetBaseClass("AliVCluster")) {
-      AliError(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloName.Data())); 
-      fClusters = 0;
-      return kFALSE;
-    }
-
-    if (!fOutClusters) {
-      fOutCaloName = fCaloName;
-      if (fCopyArray) {
-       fOutCaloName += fSuffix;
-       fOutClusters = new TClonesArray(fClusters->GetClass()->GetName(), fClusters->GetSize());
-       fOutClusters->SetName(fOutCaloName);
-       if (InputEvent()->FindListObject(fOutCaloName)) {
-         AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCaloName.Data())); 
-         return kFALSE;
-       }
-       else {
-         InputEvent()->AddObject(fOutClusters);
-       }
-      }
-      else {
-       fOutClusters = fClusters;
-      }
-    }
-
-    if (!fGeom) {
-      if (fGeomName.Length() > 0) {
-       fGeom = AliEMCALGeometry::GetInstance(fGeomName);
-       if (!fGeom)
-         AliError(Form("Could not get geometry with name %s!", fGeomName.Data()));
-      } else {
-       fGeom = AliEMCALGeometry::GetInstance();
-       if (!fGeom) 
-         AliError("Could not get default geometry!");
-      }
-    }
-    
-    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;
-}
-
-//________________________________________________________________________
 void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId)
 {
   // Get random cell.
index 523db00..f589e89 100644 (file)
@@ -7,6 +7,7 @@ class TClonesArray;
 class AliEMCALGeometry;
 class AliVCluster;
 class AliPicoTrack;
+class AliVCaloCells;
 
 #include <TH1F.h>
 #include <TF1.h>
@@ -16,37 +17,48 @@ class AliPicoTrack;
 class AliJetModelBaseTask : public AliAnalysisTaskSE {
  public:
   AliJetModelBaseTask();
-  AliJetModelBaseTask(const char *name); 
+  AliJetModelBaseTask(const char *name, Bool_t drawqa=kFALSE); 
   virtual ~AliJetModelBaseTask();
 
   void                   UserExec(Option_t* /*option*/);
+  void                   UserCreateOutputObjects();
 
-  void                   SetClusName(const char *n)            { fCaloName     = n;    }
-  void                   SetCopyArray(Bool_t copy)             { fCopyArray    = copy; }
   void                   SetEtaRange(Float_t min, Float_t max) { fEtaMin       = min;  fEtaMax = max; }
-  void                   SetGeometryName(const char *n)        { fGeomName     = n;    }
-  void                   SetMarkMC(Bool_t m)                   { fMarkMC       = m;    }
-  void                   SetNClusters(Int_t n)                 { fNClusters    = n;    }
-  void                   SetNTracks(Int_t n)                   { fNTracks      = n;    }
   void                   SetPhiRange(Float_t min, Float_t max) { fPhiMin       = min;  fPhiMax = max; }
   void                   SetPtRange(Float_t min, Float_t max)  { fPtMin        = min;  fPtMax  = max;  }
   void                   SetPtSpectrum(TH1 *f)                 { fPtSpectrum   = f;    }
   void                   SetPtSpectrum(TF1 *f)                 { fPtSpectrum   = new TH1F("ptSpectrum","ptSpectrum",250,f->GetXmin(),f->GetXmax()); 
                                                                  fPtSpectrum->Add(f); }
-  void                   SetSuffix(const char *s)              { fSuffix       = s;    }
+
+  void                   SetCopyArray(Bool_t copy)             { fCopyArray    = copy; }
   void                   SetTracksName(const char *n)          { fTracksName   = n;    }
+  void                   SetClusName(const char *n)            { fCaloName     = n;    }
+  void                   SetCellsName(const char *n)           { fCellsName    = n;    }
+  void                   SetSuffix(const char *s)              { fSuffix       = s;    }
+
+  void                   SetGeometryName(const char *n)        { fGeomName     = n;    }
+  void                   SetMarkMC(Bool_t m)                   { fMarkMC       = m;    }
+  virtual void           SetNClusters(Int_t n)                 { fNClusters    = n;    }
+  virtual void           SetNCells(Int_t n)                    { fNCells       = n;    }
+  virtual void           SetNTracks(Int_t n)                   { fNTracks      = n;    }
+
 
  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
-  AliPicoTrack          *AddTrack(Double_t pt = -1, Double_t eta = -999, Double_t phi = -1);    // add a track; if values are -1 generate random parameters
+  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
+  void                   CopyCells();
   void                   CopyClusters();
   void                   CopyTracks();
-  virtual Bool_t         ExecOnce();
   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
   Double_t               GetRandomPt();                                                         // generate a random pt value in the given range
+  virtual Bool_t         ExecOnce();                                                            // intialize task
   virtual void           Run();                                                                 // do jet model action
 
   TString                fGeomName;               // EMCal geometry name
@@ -54,6 +66,8 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE {
   TString                fOutTracksName;          // name of output track collection
   TString                fCaloName;               // name of calo cluster collection
   TString                fOutCaloName;            // name of output cluster collection
+  TString                fCellsName;              // name of calo cells collection
+  TString                fOutCellsName;           // name of output cells collection
   TString                fSuffix;                 // suffix to add in the name of new collections
   Float_t                fEtaMin;                 // eta minimum value
   Float_t                fEtaMax;                 // eta maximum value
@@ -63,20 +77,27 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE {
   Float_t                fPtMax;                  // pt maximum value
   Bool_t                 fCopyArray;              // whether or not the array will be copied to a new one before modelling
   Int_t                  fNClusters;              // how many clusters are being processed
+  Int_t                  fNCells;                 // how many cells are being processed
   Int_t                  fNTracks;                // how many tracks are being processed
   Bool_t                 fMarkMC;                 // whether or not mark new tracks/cluster as MC
   TH1                   *fPtSpectrum;             // pt spectrum parametrization to extract random pt values
+  Bool_t                 fQAhistos;               // draw QA histograms
   Bool_t                 fIsInit;                 //=true if initialized
   AliEMCALGeometry      *fGeom;                   //!pointer to EMCal geometry
   TClonesArray          *fClusters;               //!cluster collection
   TClonesArray          *fOutClusters;            //!output cluster collection
   TClonesArray          *fTracks;                 //!track collection
   TClonesArray          *fOutTracks;              //!output track collection
+  AliVCaloCells         *fCaloCells;              //!cells collection
+  AliVCaloCells         *fOutCaloCells;           //!output cells collection
+  Int_t                  fAddedCells;             //!number of added cells
+  Bool_t                 fEsdMode;                //!ESD/AOD mode
+  TList                 *fOutput;                 //!output list for QA histograms
 
  private:
   AliJetModelBaseTask(const AliJetModelBaseTask&);            // not implemented
   AliJetModelBaseTask &operator=(const AliJetModelBaseTask&); // not implemented
 
-  ClassDef(AliJetModelBaseTask, 4) // Jet modelling task
+  ClassDef(AliJetModelBaseTask, 5) // Jet modelling task
 };
 #endif
index e7955c5..aafdd5e 100644 (file)
@@ -6,20 +6,17 @@
 
 #include "AliJetResponseMaker.h"
 
-#include <TChain.h>
 #include <TClonesArray.h>
 #include <TH1F.h>
 #include <TH2F.h>
-#include <TList.h>
+#include <TH3F.h>
 #include <TLorentzVector.h>
 
-#include "AliAnalysisManager.h"
-#include "AliCentrality.h"
-#include "AliFJWrapper.h"
 #include "AliVCluster.h"
 #include "AliVTrack.h"
 #include "AliEmcalJet.h"
 #include "AliGenPythiaEventHeader.h"
+#include "AliLog.h"
 #include "AliMCEvent.h"
 
 ClassImp(AliJetResponseMaker)
@@ -47,24 +44,27 @@ AliJetResponseMaker::AliJetResponseMaker() :
   fMCJets(0),
   fHistNTrials(0),
   fHistEvents(0),
-  fHistMCJetPhiEta(0),
-  fHistMCJetsPt(0),
-  fHistMCJetPhiEtaFiducial(0),
-  fHistMCJetsPtFiducial(0),
+  fHistMCJetsPhiEta(0),
+  fHistMCJetsPtArea(0),
+  fHistMCJetsPhiEtaFiducial(0),
+  fHistMCJetsPtAreaFiducial(0),
   fHistMCJetsNEFvsPt(0),
   fHistMCJetsZvsPt(0),
-  fHistJetPhiEta(0),
-  fHistJetsPt(0),
+  fHistJetsPhiEta(0),
+  fHistJetsPtArea(0),
+  fHistJetsCorrPtArea(0),
   fHistJetsNEFvsPt(0),
   fHistJetsZvsPt(0),
-  fHistClosestDistance(0),
-  fHistClosestDeltaPhi(0),
-  fHistClosestDeltaEta(0),
-  fHistClosestDeltaPt(0),
-  fHistNonMatchedMCJetPt(0),
-  fHistNonMatchedJetPt(0),
+  fHistMatchingLevelMCPt(0),
+  fHistClosestDeltaEtaPhiMCPt(0),
+  fHistClosestDeltaPtMCPt(0),
+  fHistClosestDeltaCorrPtMCPt(0),
+  fHistNonMatchedMCJetsPtArea(0),
+  fHistNonMatchedJetsPtArea(0),
+  fHistNonMatchedJetsCorrPtArea(0),
   fHistPartvsDetecPt(0),
-  fHistMissedMCJets(0)
+  fHistPartvsDetecCorrPt(0),
+  fHistMissedMCJetsPtArea(0)
 {
   // Default constructor.
 
@@ -98,24 +98,27 @@ AliJetResponseMaker::AliJetResponseMaker(const char *name) :
   fMCJets(0),
   fHistNTrials(0),
   fHistEvents(0),
-  fHistMCJetPhiEta(0),
-  fHistMCJetsPt(0),
-  fHistMCJetPhiEtaFiducial(0),
-  fHistMCJetsPtFiducial(0),
+  fHistMCJetsPhiEta(0),
+  fHistMCJetsPtArea(0),
+  fHistMCJetsPhiEtaFiducial(0),
+  fHistMCJetsPtAreaFiducial(0),
   fHistMCJetsNEFvsPt(0),
   fHistMCJetsZvsPt(0),
-  fHistJetPhiEta(0),
-  fHistJetsPt(0),
+  fHistJetsPhiEta(0),
+  fHistJetsPtArea(0),
+  fHistJetsCorrPtArea(0),
   fHistJetsNEFvsPt(0),
   fHistJetsZvsPt(0),
-  fHistClosestDistance(0),
-  fHistClosestDeltaPhi(0),
-  fHistClosestDeltaEta(0),
-  fHistClosestDeltaPt(0),
-  fHistNonMatchedMCJetPt(0),
-  fHistNonMatchedJetPt(0),
+  fHistMatchingLevelMCPt(0),
+  fHistClosestDeltaEtaPhiMCPt(0),
+  fHistClosestDeltaPtMCPt(0),
+  fHistClosestDeltaCorrPtMCPt(0),
+  fHistNonMatchedMCJetsPtArea(0),
+  fHistNonMatchedJetsPtArea(0),
+  fHistNonMatchedJetsCorrPtArea(0),
   fHistPartvsDetecPt(0),
-  fHistMissedMCJets(0)
+  fHistPartvsDetecCorrPt(0),
+  fHistMissedMCJetsPtArea(0)
 {
   // Standard constructor.
 
@@ -165,95 +168,121 @@ void AliJetResponseMaker::UserCreateOutputObjects()
     fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
   }
 
-  fHistJetPhiEta = new TH2F("fHistJetPhiEta", "fHistJetPhiEta", 20, -2, 2, 32, 0, 6.4);
-  fHistJetPhiEta->GetXaxis()->SetTitle("#eta");
-  fHistJetPhiEta->GetYaxis()->SetTitle("#phi");
-  fOutput->Add(fHistJetPhiEta);
+  fHistJetsPhiEta = new TH2F("fHistJetsPhiEta", "fHistJetsPhiEta", 20, -2, 2, 32, 0, 6.4);
+  fHistJetsPhiEta->GetXaxis()->SetTitle("#eta");
+  fHistJetsPhiEta->GetYaxis()->SetTitle("#phi");
+  fOutput->Add(fHistJetsPhiEta);
   
-  fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinBinPt, fMaxBinPt);
-  fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-  fHistJetsPt->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistJetsPt);
+  fHistJetsPtArea = new TH2F("fHistJetsPtArea", "fHistJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
+  fHistJetsPtArea->GetXaxis()->SetTitle("area");
+  fHistJetsPtArea->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)");
+  fHistJetsPtArea->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistJetsPtArea);
+
+  fHistJetsCorrPtArea = new TH2F("fHistJetsCorrPtArea", "fHistJetsCorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, (Int_t)(1.5*fNbins), -fMaxBinPt/2, fMaxBinPt);
+  fHistJetsCorrPtArea->GetXaxis()->SetTitle("area");
+  fHistJetsCorrPtArea->GetYaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
+  fHistJetsCorrPtArea->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistJetsCorrPtArea);
   
   fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
   fHistJetsZvsPt->GetXaxis()->SetTitle("Z");
-  fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
+  fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)");
   fOutput->Add(fHistJetsZvsPt);
   
   fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
   fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
-  fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
+  fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)");
   fOutput->Add(fHistJetsNEFvsPt);
 
-  fHistMCJetPhiEta = new TH2F("fHistMCJetPhiEta", "fHistMCJetPhiEta", 20, -2, 2, 32, 0, 6.4);
-  fHistMCJetPhiEta->GetXaxis()->SetTitle("#eta");
-  fHistMCJetPhiEta->GetYaxis()->SetTitle("#phi");
-  fOutput->Add(fHistMCJetPhiEta);
+  fHistMCJetsPhiEta = new TH2F("fHistMCJetsPhiEta", "fHistMCJetsPhiEta", 20, -2, 2, 32, 0, 6.4);
+  fHistMCJetsPhiEta->GetXaxis()->SetTitle("#eta");
+  fHistMCJetsPhiEta->GetYaxis()->SetTitle("#phi");
+  fOutput->Add(fHistMCJetsPhiEta);
   
-  fHistMCJetsPt = new TH1F("fHistMCJetsPt", "fHistMCJetsPt", fNbins, fMinBinPt, fMaxBinPt);
-  fHistMCJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-  fHistMCJetsPt->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistMCJetsPt);
-
-  fHistMCJetPhiEtaFiducial = new TH2F("fHistMCJetPhiEtaFiducial", "fHistMCJetPhiEtaFiducial", 20, -2, 2, 32, 0, 6.4);
-  fHistMCJetPhiEtaFiducial->GetXaxis()->SetTitle("#eta");
-  fHistMCJetPhiEtaFiducial->GetYaxis()->SetTitle("#phi");
-  fOutput->Add(fHistMCJetPhiEtaFiducial);
+  fHistMCJetsPtArea = new TH2F("fHistMCJetsPtArea", "fHistMCJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
+  fHistMCJetsPtArea->GetXaxis()->SetTitle("area");
+  fHistMCJetsPtArea->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
+  fHistMCJetsPtArea->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistMCJetsPtArea);
+
+  fHistMCJetsPhiEtaFiducial = new TH2F("fHistMCJetsPhiEtaFiducial", "fHistMCJetsPhiEtaFiducial", 20, -2, 2, 32, 0, 6.4);
+  fHistMCJetsPhiEtaFiducial->GetXaxis()->SetTitle("#eta");
+  fHistMCJetsPhiEtaFiducial->GetYaxis()->SetTitle("#phi");
+  fOutput->Add(fHistMCJetsPhiEtaFiducial);
   
-  fHistMCJetsPtFiducial = new TH1F("fHistMCJetsPtFiducial", "fHistMCJetsPtFiducial", fNbins, fMinBinPt, fMaxBinPt);
-  fHistMCJetsPtFiducial->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-  fHistMCJetsPtFiducial->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistMCJetsPtFiducial);
+  fHistMCJetsPtAreaFiducial = new TH2F("fHistMCJetsPtAreaFiducial", "fHistMCJetsPtAreaFiducial", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
+  fHistMCJetsPtAreaFiducial->GetXaxis()->SetTitle("area");
+  fHistMCJetsPtAreaFiducial->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
+  fHistMCJetsPtAreaFiducial->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistMCJetsPtAreaFiducial);
   
   fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
   fHistMCJetsZvsPt->GetXaxis()->SetTitle("Z");
-  fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
+  fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
   fOutput->Add(fHistMCJetsZvsPt);
 
   fHistMCJetsNEFvsPt = new TH2F("fHistMCJetsNEFvsPt", "fHistMCJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
   fHistMCJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
-  fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
+  fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
   fOutput->Add(fHistMCJetsNEFvsPt);
 
-  fHistClosestDistance = new TH1F("fHistClosestDistance", "fHistClosestDistance", 50, 0, fMaxDistance * 1.2);
-  fHistClosestDistance->GetXaxis()->SetTitle("d");
-  fHistClosestDistance->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistClosestDistance);
-
-  fHistClosestDeltaPhi = new TH1F("fHistClosestDeltaPhi", "fHistClosestDeltaPhi", 64, -1.6, 4.8);
-  fHistClosestDeltaPhi->GetXaxis()->SetTitle("#Delta#phi");
-  fHistClosestDeltaPhi->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistClosestDeltaPhi);
-
-  fHistClosestDeltaEta = new TH1F("fHistClosestDeltaEta", "fHistClosestDeltaEta", TMath::CeilNint(fJetMaxEta - fJetMinEta) * 20, fJetMinEta * 2, fJetMaxEta * 2);
-  fHistClosestDeltaEta->GetXaxis()->SetTitle("#Delta#eta");
-  fHistClosestDeltaEta->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistClosestDeltaEta);
-
-  fHistClosestDeltaPt = new TH1F("fHistClosestDeltaPt", "fHistClosestDeltaPt", fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
-  fHistClosestDeltaPt->GetXaxis()->SetTitle("#Delta p_{T}");
-  fHistClosestDeltaPt->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistClosestDeltaPt);
-
-  fHistNonMatchedMCJetPt = new TH1F("fHistNonMatchedMCJetPt", "fHistNonMatchedMCJetPt", fNbins, fMinBinPt, fMaxBinPt);
-  fHistNonMatchedMCJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-  fHistNonMatchedMCJetPt->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistNonMatchedMCJetPt);
-
-  fHistNonMatchedJetPt = new TH1F("fHistNonMatchedJetPt", "fHistNonMatchedJetPt", fNbins, fMinBinPt, fMaxBinPt);
-  fHistNonMatchedJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-  fHistNonMatchedJetPt->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistNonMatchedJetPt);
+  fHistMatchingLevelMCPt = new TH2F("fHistMatchingLevelMCPt", "fHistMatchingLevelMCPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
+  fHistMatchingLevelMCPt->GetXaxis()->SetTitle("Matching level");
+  fHistMatchingLevelMCPt->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
+  fOutput->Add(fHistMatchingLevelMCPt);
+
+  fHistClosestDeltaEtaPhiMCPt = new TH3F("fHistClosestDeltaEtaPhiMCPt", "fHistClosestDeltaEtaPhiMCPt", TMath::CeilNint(fJetMaxEta - fJetMinEta) * 20, fJetMinEta * 2, fJetMaxEta * 2, 64, -1.6, 4.8, fNbins, fMinBinPt, fMaxBinPt);
+  fHistClosestDeltaEtaPhiMCPt->GetXaxis()->SetTitle("#Delta#eta");
+  fHistClosestDeltaEtaPhiMCPt->GetYaxis()->SetTitle("#Delta#phi");
+  fHistClosestDeltaEtaPhiMCPt->GetZaxis()->SetTitle("p_{T}^{gen}");
+  fOutput->Add(fHistClosestDeltaEtaPhiMCPt);
+
+  fHistClosestDeltaPtMCPt = new TH2F("fHistClosestDeltaPtMCPt", "fHistClosestDeltaPtMCPt", fNbins, fMinBinPt, fMaxBinPt, (Int_t)(fNbins*1.5), -fMaxBinPt / 2, fMaxBinPt);
+  fHistClosestDeltaPtMCPt->GetXaxis()->SetTitle("p_{T}^{gen}");  
+  fHistClosestDeltaPtMCPt->GetYaxis()->SetTitle("#Deltap_{T}^{rec} (GeV/c)");
+  fHistClosestDeltaPtMCPt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistClosestDeltaPtMCPt);
+
+  fHistClosestDeltaCorrPtMCPt = new TH2F("fHistClosestDeltaCorrPtMCPt", "fHistClosestDeltaCorrPtMCPt", fNbins, fMinBinPt, fMaxBinPt, (Int_t)(fNbins*1.5), -fMaxBinPt / 2, fMaxBinPt);
+  fHistClosestDeltaCorrPtMCPt->GetXaxis()->SetTitle("p_{T}^{gen}");  
+  fHistClosestDeltaCorrPtMCPt->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
+  fHistClosestDeltaCorrPtMCPt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistClosestDeltaCorrPtMCPt);
+
+  fHistNonMatchedMCJetsPtArea = new TH2F("fHistNonMatchedMCJetsPtArea", "fHistNonMatchedMCJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
+  fHistNonMatchedMCJetsPtArea->GetXaxis()->SetTitle("area");
+  fHistNonMatchedMCJetsPtArea->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
+  fHistNonMatchedMCJetsPtArea->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistNonMatchedMCJetsPtArea);
+
+  fHistNonMatchedJetsPtArea = new TH2F("fHistNonMatchedJetsPtArea", "fHistNonMatchedJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
+  fHistNonMatchedJetsPtArea->GetXaxis()->SetTitle("area");
+  fHistNonMatchedJetsPtArea->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)");
+  fHistNonMatchedJetsPtArea->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistNonMatchedJetsPtArea);
+
+  fHistNonMatchedJetsCorrPtArea = new TH2F("fHistNonMatchedJetsCorrPtArea", "fHistNonMatchedJetsCorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, (Int_t)(1.5*fNbins), -fMaxBinPt/2, fMaxBinPt);
+  fHistNonMatchedJetsCorrPtArea->GetXaxis()->SetTitle("area");
+  fHistNonMatchedJetsCorrPtArea->GetYaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
+  fHistNonMatchedJetsCorrPtArea->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistNonMatchedJetsCorrPtArea);
 
   fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
   fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{rec}");
   fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{gen}");
   fOutput->Add(fHistPartvsDetecPt);
 
-  fHistMissedMCJets = new TH1F("fHistMissedMCJets", "fHistMissedMCJets", fNbins, fMinBinPt, fMaxBinPt);
-  fHistMissedMCJets->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-  fHistMissedMCJets->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistMissedMCJets);
+  fHistPartvsDetecCorrPt = new TH2F("fHistPartvsDetecCorrPt", "fHistPartvsDetecCorrPt", (Int_t)(1.5*fNbins), -fMaxBinPt/2, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
+  fHistPartvsDetecCorrPt->GetXaxis()->SetTitle("p_{T}^{corr}");
+  fHistPartvsDetecCorrPt->GetYaxis()->SetTitle("p_{T}^{gen}");
+  fOutput->Add(fHistPartvsDetecCorrPt);
+
+  fHistMissedMCJetsPtArea = new TH2F("fHistMissedMCJetsPtArea", "fHistMissedMCJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
+  fHistMissedMCJetsPtArea->GetXaxis()->SetTitle("area");  
+  fHistMissedMCJetsPtArea->GetYaxis()->SetTitle("p_{T} (GeV/c)");
+  fHistMissedMCJetsPtArea->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistMissedMCJetsPtArea);
 
   PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
 }
@@ -506,30 +535,30 @@ Bool_t AliJetResponseMaker::FillHistograms()
       if (!AcceptBiasJet(jet->MatchedJet()) || 
          jet->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet->MatchedJet()->MaxClusterPt() > fMaxClusterPt ||
          jet->MatchedJet()->Pt() > fMaxBinPt) {
-       fHistMissedMCJets->Fill(jet->Pt(), fEventWeight);
+       fHistMissedMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight);
       }
       else {
-       fHistClosestDistance->Fill(jet->ClosestJetDistance(), fEventWeight);
+       fHistMatchingLevelMCPt->Fill(jet->ClosestJetDistance(), jet->Pt(), fEventWeight);
 
        Double_t deta = jet->MatchedJet()->Eta() - jet->Eta();
-       fHistClosestDeltaEta->Fill(deta, fEventWeight);
-
        Double_t dphi = jet->MatchedJet()->Phi() - jet->Phi();
-       fHistClosestDeltaPhi->Fill(dphi, fEventWeight);
+       fHistClosestDeltaEtaPhiMCPt->Fill(deta, dphi, jet->Pt(), fEventWeight);
 
        Double_t dpt = jet->MatchedJet()->Pt() - jet->Pt();
-       fHistClosestDeltaPt->Fill(dpt, fEventWeight);
+       fHistClosestDeltaPtMCPt->Fill(jet->Pt(), dpt, fEventWeight);
+       fHistClosestDeltaCorrPtMCPt->Fill(jet->Pt(), dpt - fRhoVal * jet->MatchedJet()->Area(), fEventWeight);
 
        fHistPartvsDetecPt->Fill(jet->MatchedJet()->Pt(), jet->Pt(), fEventWeight);
+       fHistPartvsDetecCorrPt->Fill(jet->MatchedJet()->Pt() - fRhoVal * jet->MatchedJet()->Area(), jet->Pt(), fEventWeight);
       }
     }
     else {
-      fHistNonMatchedMCJetPt->Fill(jet->Pt(), fEventWeight);
-      fHistMissedMCJets->Fill(jet->Pt(), fEventWeight);
+      fHistNonMatchedMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight);
+      fHistMissedMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight);
     }
 
-    fHistMCJetsPt->Fill(jet->Pt(), fEventWeight);
-    fHistMCJetPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
+    fHistMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight);
+    fHistMCJetsPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
 
     fHistMCJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight);
       
@@ -544,8 +573,8 @@ Bool_t AliJetResponseMaker::FillHistograms()
     if (jet->Eta() < fJetMinEta || jet->Eta() > fJetMaxEta || jet->Phi() < fJetMinPhi || jet->Phi() > fJetMaxPhi)
       continue;
     
-    fHistMCJetsPtFiducial->Fill(jet->Pt(), fEventWeight);
-    fHistMCJetPhiEtaFiducial->Fill(jet->Eta(), jet->Phi(), fEventWeight);
+    fHistMCJetsPtAreaFiducial->Fill(jet->Area(), jet->Pt(), fEventWeight);
+    fHistMCJetsPhiEtaFiducial->Fill(jet->Eta(), jet->Phi(), fEventWeight);
   }
 
   const Int_t nJets = fJets->GetEntriesFast();
@@ -569,12 +598,14 @@ Bool_t AliJetResponseMaker::FillHistograms()
       continue;
 
     if (!jet->MatchedJet()) {
-      fHistNonMatchedJetPt->Fill(jet->Pt(), fEventWeight);
+      fHistNonMatchedJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight);
+      fHistNonMatchedJetsCorrPtArea->Fill(jet->Area(), jet->Pt() - fRhoVal * jet->Area(), fEventWeight);
     }
 
-    fHistJetsPt->Fill(jet->Pt(), fEventWeight);
+    fHistJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight);
+    fHistJetsCorrPtArea->Fill(jet->Area(), jet->Pt() - fRhoVal * jet->Area(), fEventWeight);
 
-    fHistJetPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
+    fHistJetsPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
 
     fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight);
 
index b25f0ce..629ef76 100644 (file)
@@ -5,8 +5,9 @@
 
 class AliGenPythiaEventHeader;
 class TClonesArray;
-class TH1F;
-class TH2F;
+class TH1;
+class TH2;
+class TH3;
 
 #include "AliAnalysisTaskEmcalJet.h"
 
@@ -36,55 +37,58 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet {
   Bool_t                      RetrieveEventObjects();
   Bool_t                      Run();
 
-  TString                     fMCTracksName;              // name of MC particle collection
-  TString                     fMCJetsName;                // name of MC jet collection
-  Double_t                    fMaxDistance;               // maximum distance between matched particle and detector level jets
-  Bool_t                      fDoWeighting;               // = true, weight using trials and given x section
-  Bool_t                      fEventWeightHist;           // = true create event weight histogram
-  Bool_t                      fMCFiducial;                // = true MC jets in fiducial acceptance
-  Double_t                    fMCminEta;                  // MC jets minimum eta
-  Double_t                    fMCmaxEta;                  // MC jets maximum eta
-  Double_t                    fMCminPhi;                  // MC jets minimum phi
-  Double_t                    fMCmaxPhi;                  // MC jets maximum phi
-  Int_t                       fSelectPtHardBin;           // select one pt hard bin for analysis
-  Bool_t                      fDoMatching;                // whether or not it should run the matching between MC and rec jets
+  TString                     fMCTracksName;                  // name of MC particle collection
+  TString                     fMCJetsName;                    // name of MC jet collection
+  Double_t                    fMaxDistance;                   // maximum distance between matched particle and detector level jets
+  Bool_t                      fDoWeighting;                   // = true, weight using trials and given x section
+  Bool_t                      fEventWeightHist;               // = true create event weight histogram
+  Bool_t                      fMCFiducial;                    // = true MC jets in fiducial acceptance
+  Double_t                    fMCminEta;                      // MC jets minimum eta
+  Double_t                    fMCmaxEta;                      // MC jets maximum eta
+  Double_t                    fMCminPhi;                      // MC jets minimum phi
+  Double_t                    fMCmaxPhi;                      // MC jets maximum phi
+  Int_t                       fSelectPtHardBin;               // select one pt hard bin for analysis
+  Bool_t                      fDoMatching;                    // whether or not it should run the matching between MC and rec jets
 
-  AliGenPythiaEventHeader    *fPythiaHeader;              //!event Pythia header
-  Double_t                    fEventWeight;               //!event weight
-  Int_t                       fPtHardBin;                 //!event pt hard bin
-  Int_t                       fNTrials;                   //!event trials
-  TClonesArray               *fMCTracks;                  //!MC particles
-  TClonesArray               *fMCJets;                    //!MC jets
+  AliGenPythiaEventHeader    *fPythiaHeader;                  //!event Pythia header
+  Double_t                    fEventWeight;                   //!event weight
+  Int_t                       fPtHardBin;                     //!event pt hard bin
+  Int_t                       fNTrials;                       //!event trials
+  TClonesArray               *fMCTracks;                      //!MC particles
+  TClonesArray               *fMCJets;                        //!MC jets
   // General histograms
-  TH1F                       *fHistNTrials;               //!total number of trials per pt hard bin
-  TH1F                       *fHistEvents;                //!total number of events per pt hard bin
-  TH1F                       *fHistEventWeight[11];       //!event weight
+  TH1                        *fHistNTrials;                   //!total number of trials per pt hard bin
+  TH1                        *fHistEvents;                    //!total number of events per pt hard bin
+  TH1                        *fHistEventWeight[11];           //!event weight
   // Particle level jets
-  TH2F                       *fHistMCJetPhiEta;           //!phi-eta distribution of jets
-  TH1F                       *fHistMCJetsPt;              //!inclusive jet pt spectrum
-  TH2F                       *fHistMCJetPhiEtaFiducial;   //!phi-eta distribution of jets in fiducial acceptance (plus lead hadron bias)
-  TH1F                       *fHistMCJetsPtFiducial;      //!inclusive jet pt spectrum in fiducial acceptance (plus lead hadron bias)
-  TH2F                       *fHistMCJetsNEFvsPt;         //!jet neutral energy fraction vs. jet pt
-  TH2F                       *fHistMCJetsZvsPt;           //!constituent Pt over Jet Pt ratio vs. jet pt
+  TH2                        *fHistMCJetsPhiEta;              //!phi-eta distribution of jets
+  TH2                        *fHistMCJetsPtArea;              //!inclusive jet pt vs area histogram
+  TH2                        *fHistMCJetsPhiEtaFiducial;      //!phi-eta distribution of jets in fiducial acceptance (plus lead hadron bias)
+  TH2                        *fHistMCJetsPtAreaFiducial;      //!inclusive jet pt spectrum in fiducial acceptance (plus lead hadron bias)
+  TH2                        *fHistMCJetsNEFvsPt;             //!jet neutral energy fraction vs. jet pt
+  TH2                        *fHistMCJetsZvsPt;               //!constituent Pt over Jet Pt ratio vs. jet pt
   // Detector level jets
-  TH2F                       *fHistJetPhiEta;             //!phi-eta distribution of jets
-  TH1F                       *fHistJetsPt;                //!inclusive jet pt spectrum
-  TH2F                       *fHistJetsNEFvsPt;           //!jet neutral energy fraction vs. jet pt
-  TH2F                       *fHistJetsZvsPt;             //!constituent Pt over Jet Pt ratio vs. jet pt
+  TH2                        *fHistJetsPhiEta;                //!phi-eta distribution of jets
+  TH2                        *fHistJetsPtArea;                //!inclusive jet pt vs. area histogram
+  TH2                        *fHistJetsCorrPtArea;            //!inclusive jet pt vs. area histogram
+  TH2                        *fHistJetsNEFvsPt;               //!jet neutral energy fraction vs. jet pt
+  TH2                        *fHistJetsZvsPt;                 //!constituent Pt over Jet Pt ratio vs. jet pt
   // Detector-particle level matching
-  TH1F                       *fHistClosestDistance;       //!distance between closest particle to detector level jet
-  TH1F                       *fHistClosestDeltaPhi;       //!delta phi between closest particle to detector level jet
-  TH1F                       *fHistClosestDeltaEta;       //!delta eta between closest particle to detector level jet
-  TH1F                       *fHistClosestDeltaPt;        //!delta pt between closest particle to detector level jet
-  TH1F                       *fHistNonMatchedMCJetPt;     //!non-matched mc jet pt distribution
-  TH1F                       *fHistNonMatchedJetPt;       //!non-matched jet pt distribution
-  TH2F                       *fHistPartvsDetecPt;         //!particle vs detector level jet pt
-  TH1F                       *fHistMissedMCJets;          //!mc jets not measured
+  TH2                        *fHistMatchingLevelMCPt;         //!matching level vs. particle level pt
+  TH3                        *fHistClosestDeltaEtaPhiMCPt;    //!delta eta-phi between closest particle level to detector level jet vs. particle level pt
+  TH2                        *fHistClosestDeltaPtMCPt;        //!delta pt between closest particle level to detector level jet vs. particle level pt
+  TH2                        *fHistClosestDeltaCorrPtMCPt;    //!delta pt between closest particle level to detector level jet vs. particle level pt
+  TH2                        *fHistNonMatchedMCJetsPtArea;    //!non-matched mc jet pt distribution
+  TH2                        *fHistNonMatchedJetsPtArea;      //!non-matched jet pt distribution
+  TH2                        *fHistNonMatchedJetsCorrPtArea;  //!non-matched jet pt distribution
+  TH2                        *fHistPartvsDetecPt;             //!particle vs detector level jet pt
+  TH2                        *fHistPartvsDetecCorrPt;         //!particle vs detector level jet pt
+  TH2                        *fHistMissedMCJetsPtArea;         //!mc jets not measured
 
  private:
   AliJetResponseMaker(const AliJetResponseMaker&);            // not implemented
   AliJetResponseMaker &operator=(const AliJetResponseMaker&); // not implemented
 
-  ClassDef(AliJetResponseMaker, 8) // Jet response matrix producing task
+  ClassDef(AliJetResponseMaker, 9) // Jet response matrix producing task
 };
 #endif
index bfb5db7..bb282c0 100644 (file)
@@ -25,6 +25,7 @@
 #include "AliLog.h"
 #include "AliEMCALGeometry.h"
 #include "AliEMCALGeoParams.h"
+#include "AliPicoTrack.h"
 
 #include "AliAnalysisTaskSAQA.h"
 
@@ -498,12 +499,16 @@ Float_t AliAnalysisTaskSAQA::DoTrackLoop()
     }
     else {
       fHistTrPhiEtaPt[fCentBin][3]->Fill(track->Eta(), track->Phi(), track->Pt());
+      Int_t type = 0;
 
-      Int_t label = track->GetLabel();
-      if (label >= 0 && label < 3)
-       fHistTrPhiEtaPt[fCentBin][label]->Fill(track->Eta(), track->Phi(), track->Pt());
+      AliPicoTrack* ptrack = dynamic_cast<AliPicoTrack*>(track);
+      if (ptrack)
+       type = ptrack->GetTrackType();
+
+      if (type >= 0 && type < 3)
+       fHistTrPhiEtaPt[fCentBin][type]->Fill(track->Eta(), track->Phi(), track->Pt());
       else
-       AliWarning(Form("%s: track label %d not recognized!", GetName(), label));
+       AliWarning(Form("%s: track type %d not recognized!", GetName(), type));
     }
 
     if (!vtrack)
diff --git a/PWGJE/EMCALJetTasks/macros/AddTaskJetEmbeddingFromAOD.C b/PWGJE/EMCALJetTasks/macros/AddTaskJetEmbeddingFromAOD.C
new file mode 100644 (file)
index 0000000..1f0df3c
--- /dev/null
@@ -0,0 +1,140 @@
+// $Id: AddTaskJetEmbeddingFromAOD.C  $
+
+TObjArray* GenerateFileList(const char* list, Int_t nFiles);
+
+AliJetEmbeddingFromAODTask* AddTaskJetEmbeddingFromAOD(
+  const char     *tracksName    = "Tracks",
+  const char     *cellsName     = "EMCALCells",
+  const char     *fileList      = "files.txt",
+  const char     *aodTreeName   = "aodTree",
+  const char     *aodTracksName = "tracks",
+  const char     *aodCellsName  = "emcalCells",
+  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     nCells        = 1234567890,
+  const Bool_t    copyArray     = kTRUE,
+  const Int_t     nFiles        = 1234567890,
+  const Bool_t    makeQA        = kFALSE,
+  const Double_t  minPt         = 0,
+  const Double_t  maxPt         = 1000,
+  const Double_t  minEta        = -0.9,
+  const Double_t  maxEta        = 0.9,
+  const Double_t  minPhi        = 0,
+  const Double_t  maxPhi        = TMath::Pi() * 2,
+  const char     *taskName      = "JetEmbeddingFromAODTask"
+)
+{  
+  // Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr)
+  {
+    ::Error("AddTaskJetEmbeddingFromAOD", "No analysis manager to connect to.");
+    return NULL;
+  }  
+  
+  // Check the analysis type using the event handlers connected to the analysis manager.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler())
+  {
+    ::Error("AddTaskJetEmbeddingFromAOD", "This task requires an input event handler");
+    return NULL;
+  }
+  
+  //-------------------------------------------------------
+  // Init the task and do settings
+  //-------------------------------------------------------
+
+  AliJetEmbeddingFromAODTask *jetEmb = new AliJetEmbeddingFromAODTask(taskName,makeQA);
+  jetEmb->SetTracksName(tracksName);
+  jetEmb->SetCellsName(cellsName);
+  jetEmb->SetFileList(GenerateFileList(fileList, nFiles));
+  jetEmb->SetAODTreeName(aodTreeName);
+  jetEmb->SetAODTracksName(aodTracksName);
+  jetEmb->SetAODCellsName(aodCellsName);
+  jetEmb->SetCentralityRange(minCent, maxCent);
+  jetEmb->SetTriggerMask(mask);
+  jetEmb->SetNCells(nCells);
+  jetEmb->SetNTracks(nTracks);
+  jetEmb->SetCopyArray(copyArray);
+  jetEmb->SetEtaRange(minEta, maxEta);
+  jetEmb->SetPhiRange(minPhi, maxPhi);
+  jetEmb->SetPtRange(minPt, maxPt);
+
+  jetEmb->SetIncludeNoITS(includeNoITS);
+  TString runPeriod(runperiod);
+  runPeriod.ToLower();
+  if (runPeriod == "lhc11h") {
+    jetEmb->SetAODfilterBits(256,512); // hybrid tracks for LHC11h
+  }
+  else if (runPeriod == "lhc11a" || runPeriod == "lhc12a15a" || runPeriod == "lhc12a15e") {
+    jetEmb->SetAODfilterBits(256,16); // hybrid tracks for LHC11a, LHC12a15a and LHC12a15e
+  }
+  else {
+    if (runPeriod.IsNull())
+      ::Warning("Run period %s not known. It will use IsHybridGlobalConstrainedGlobal.");
+  }
+
+  //-------------------------------------------------------
+  // Final settings, pass to manager and set the containers
+  //-------------------------------------------------------
+
+  mgr->AddTask(jetEmb);
+    
+  // Create containers for input/output
+  mgr->ConnectInput(jetEmb, 0, mgr->GetCommonInputContainer());
+
+  if (makeQA) {
+    TString contName = taskName;
+    contName += "_histos";
+    AliAnalysisDataContainer *outc = mgr->CreateContainer(contName,
+                                                         TList::Class(),
+                                                         AliAnalysisManager::kOutputContainer,
+                                                         "AnalysisResults.root");
+    mgr->ConnectOutput(jetEmb, 1, outc);
+  }
+
+  return jetEmb;
+}
+
+TObjArray* GenerateFileList(const char* list, Int_t nFiles)
+{
+  TObjArray *array = new TObjArray(9999);
+
+  TString myList = list;
+  if (myList.Contains("alien:///")) {
+    TFile::Cp(myList,"file:./list.txt");
+    myList = "./list.txt";
+  }
+
+  // Open the input stream
+  ifstream in;
+  in.open(myList.Data());
+
+  Int_t count = 0;
+
+  // Read the input list of files and add them to the chain
+  TString line;
+  while (in.good()) {
+    if (nFiles != 1234567890) {
+      if (count >= nFiles)
+       break;
+    }
+
+    in >> line;
+
+    if (line.Length() == 0)
+      continue;
+
+    TObjString *aodFile = new TObjString(line);
+    array->Add(aodFile);
+    
+    count++;
+  }
+
+  return array;
+}
index 5123657..4adf515 100644 (file)
@@ -3,6 +3,7 @@
 AliJetResponseMaker* AddTaskJetRespPtHard(const char *ntracks            = "Tracks",
                                          const char *nclusters          = "CaloClusters",
                                          const char *njets              = "Jets",
+                                         const char *nrho               = "Rho",
                                          const char *nmctracks          = "MCParticles",
                                          const char *nmcjets            = "MCJets",
                                          Double_t    jetradius          = 0.2,
@@ -25,7 +26,7 @@ AliJetResponseMaker* AddTaskJetRespPtHard(const char *ntracks            = "Trac
   AliJetResponseMaker *jetTask = new AliJetResponseMaker[maxPtBin - minPtBin + 1];
 
   for (Int_t i = minPtBin; i <= maxPtBin; i++) {
-    AddTaskJetResponseMaker(ntracks, nclusters, njets, nmctracks, nmcjets,
+    AddTaskJetResponseMaker(ntracks, nclusters, njets, nrho, nmctracks, nmcjets,
                            jetradius, jetptcut, jetareacut, jetBiasTrack, 
                            jetBiasClus, maxDistance, type, i, taskname, jetTask + i - minPtBin);
     jetTask[i - minPtBin].SetDoMatching(domatch);
index 005f58d..d0724f6 100644 (file)
@@ -4,6 +4,7 @@ AliJetResponseMaker* AddTaskJetResponseMaker(
   const char *ntracks            = "Tracks",
   const char *nclusters          = "CaloClusters",
   const char *njets              = "Jets",
+  const char *nrho               = "Rho",
   const char *nmctracks          = "MCParticles",
   const char *nmcjets            = "MCJets",
   Double_t    jetradius          = 0.2,
@@ -61,6 +62,7 @@ AliJetResponseMaker* AddTaskJetResponseMaker(
   jetTask->SetTracksName(ntracks);
   jetTask->SetClusName(nclusters);
   jetTask->SetJetsName(njets);
+  jetTask->SetRhoName(nrho);
   jetTask->SetMCJetsName(nmcjets);
   jetTask->SetMCTracksName(nmctracks);
   jetTask->SetJetRadius(jetradius);
index 9cbd888..210d81f 100644 (file)
@@ -18,6 +18,7 @@
 #pragma link C++ class AliHadCorrTask+;
 #pragma link C++ class AliJetEmbeddingTask+;
 #pragma link C++ class AliJetEmbeddingFromGenTask+;
+#pragma link C++ class AliJetEmbeddingFromAODTask+;
 #pragma link C++ class AliJetModelBaseTask+;
 #pragma link C++ class AliJetRandomizerTask+;
 #pragma link C++ class AliJetResponseMaker+;